import numpy as np, json
from scipy.stats import norm,binom
from scipy.integrate import quad
from pathlib import Path
p=.02;rho=.15;g=.45;alpha=.999;ead=100
x=norm.ppf(alpha);a=np.sqrt(rho/(1-rho));t=(norm.ppf(p)+np.sqrt(rho)*x)/np.sqrt(1-rho);q=norm.cdf(t);dp=a*norm.pdf(t)
v=g*g*.001*q*(1-q);vp=g*g*.001*dp*(1-2*q);ga=-(vp+v*(-x+a*t))/(2*g*dp)
def finite_cdf(k):
 return quad(lambda z:binom.cdf(k,1000,norm.cdf((norm.ppf(p)+np.sqrt(rho)*z)/np.sqrt(1-rho)))*norm.pdf(z),-10,10,epsabs=1e-11,epsrel=1e-10)[0]
k=int(np.floor(1000*q))
while finite_cdf(k)<alpha:k+=1
while finite_cdf(k-1)>=alpha:k-=1
# Synthetic migration data; row/column ordering AAA,AA,A,BBB,BB,B,CCC,D.
P=np.array([[.97,.02,.006,.002,.001,.0005,.0003,.0002],[.005,.965,.02,.005,.002,.001,.001,.001],[.001,.012,.95,.025,.006,.003,.001,.002],[.001,.003,.04,.89,.045,.014,.005,.002],[0,.001,.005,.04,.88,.05,.014,.01],[0,0,.001,.005,.04,.86,.064,.03],[0,0,0,.002,.008,.08,.81,.10],[0,0,0,0,0,0,0,1]])
assert np.allclose(P.sum(axis=1),1)
yields=np.array([.03,.032,.035,.04,.05,.07,.10]);values=np.r_[6+sum(6/(1+yields)**j for j in range(1,5))+100/(1+yields)**4,40.]
probs=P[3];B=values[3];mu=probs@values
loss=B-values;ix=np.argsort(loss);cum=np.cumsum(probs[ix]);var={str(al):float(loss[ix[np.searchsorted(cum,al)]]) for al in [.99,.999]}
# Two BBB bonds, latent Gaussian correlation .25; worst-to-best thresholds.
pr=probs[::-1];vv=values[::-1];edges=np.r_[-np.inf,norm.ppf(np.cumsum(pr)[:-1]),np.inf];joint=np.zeros((8,8));r=.25
for i in range(8):
 for j in range(8):
  joint[i,j]=quad(lambda z:norm.pdf(z)*(norm.cdf((edges[j+1]-r*z)/np.sqrt(1-r*r))-norm.cdf((edges[j]-r*z)/np.sqrt(1-r*r))),edges[i],edges[i+1],epsabs=1e-12)[0]
assert np.allclose(joint.sum(axis=0),pr,atol=1e-9)
ll=(2*B-(vv[:,None]+vv[None,:])).ravel();pp=joint.ravel();ii=np.argsort(ll);cc=np.cumsum(pp[ii]);two={str(al):float(ll[ii[np.searchsorted(cc,al)]]) for al in [.99,.999]}
rates=np.array([.0025,.005,.0075,.01,.015,.02,.03,.04,.05,.08]);zz=norm.ppf(rates);m=float(zz.mean());s2=float(zz.var());rh=s2/(1+s2);pd=norm.cdf(m/np.sqrt(1+s2))
A=150;D=100;sig=.25;T=1;riskfree=.03;drift=.05
d1=(np.log(A/D)+(riskfree+.5*sig**2)*T)/(sig*np.sqrt(T));d2=d1-sig*np.sqrt(T);E=A*norm.cdf(d1)-D*np.exp(-riskfree*T)*norm.cdf(d2);se=norm.cdf(d1)*A/E*sig;dd=(np.log(A/D)+(drift-.5*sig**2)*T)/(sig*np.sqrt(T))
# Order-statistic MC quantile confidence ranks only, not model parameter uncertainty.
lo=int(binom.ppf(.025,1000000,.999));hi=int(binom.ppf(.975,1000000,.999))+1
out=dict(ga=dict(q=q,asymptotic=ead*g*q,correction=ead*ga,adjusted=ead*(g*q+ga),finite_k=k,finite_var=ead*g*k/1000,cdf_prev=finite_cdf(k-1),cdf_k=finite_cdf(k)),migration=dict(matrix=P.tolist(),values=values.tolist(),mean=mu,benchmark=B,el=B-mu,var=var,credit_var={al:v-(B-mu) for al,v in var.items()},two_var=two,two_credit={al:v-2*(B-mu) for al,v in two.items()},joint_default=joint[0,0]),calibration=dict(mean=m,var=s2,rho=rh,pd=pd,rates=rates.tolist()),merton=dict(equity=E,equityvol=se,dd=dd,pd=norm.cdf(-dd),riskneutral_pd=norm.cdf(-d2)),ci_ranks=[lo,hi])
Path('credit-var-worked-results.json').write_text(json.dumps(out,indent=2));print(json.dumps(out,indent=2))