"""Constructed teaching examples. Python 3, NumPy, SciPy. Amounts in crore unless stated."""
import json
import numpy as np
from scipy.stats import norm,t,beta,binom
from scipy.special import ndtr,ndtri
from scipy.integrate import quad
from pathlib import Path
N=1_000_000; alpha=.999; p=.02; rho=.15; ead=100; lgd=.45
rng=np.random.default_rng(20261010); z=rng.standard_normal(N); pc=ndtr((ndtri(p)-np.sqrt(rho)*z)/np.sqrt(1-rho)); count=rng.binomial(1000,pc)
def measures(x):
 x=np.asarray(x); ordered=np.sort(x); k=int(np.ceil(alpha*len(x)))-1;q=float(ordered[k]); es=float((np.mean(x*(x>q))+q*(np.mean(x<=q)-alpha))/(1-alpha)); ranks=[int(binom.ppf(.025,len(x),alpha)),int(binom.ppf(.975,len(x),alpha))+1]
 return {'el':float(x.mean()),'var':q,'es':es,'credit_var':q-float(x.mean()),'ci95':[float(ordered[ranks[0]-1]),float(ordered[ranks[1]-1])],'ci_ranks':ranks}
base=measures(ead/1000*lgd*count)
eta=rng.standard_normal(N)
# Portfolio-common severity: Beta(9,11). Independent severity versus dependence on bad economy.
ind_lgd=beta.ppf(ndtr(eta),9,11); dep_lgd=beta.ppf(ndtr(-np.sqrt(.5)*z+np.sqrt(.5)*eta),9,11)
random_ind=measures(ead/1000*count*ind_lgd); random_dep=measures(ead/1000*count*dep_lgd)
# Common chi-square mixing variable is shared across all obligors: multivariate t copula.
nu=5; mix=rng.chisquare(nu,N); tp=ndtr((t.ppf(p,nu)*np.sqrt(mix/nu)-np.sqrt(rho)*z)/np.sqrt(1-rho)); tn=rng.binomial(1000,tp); tc=measures(ead/1000*lgd*tn); taildep=2*t.cdf(-np.sqrt((nu+1)*(1-rho)/(1+rho)),nu+1)
c=ndtri(p); cutoff=-ndtri(alpha)
biv=quad(lambda zz:norm.pdf(zz)*ndtr((c-np.sqrt(rho)*zz)/np.sqrt(1-rho)),-12,cutoff,epsabs=1e-12)[0]; analytic_es=ead*lgd*biv/(1-alpha)
sectors=[]
for name,ex,pd,cor,severity in [('MSME services',60,.02,.15,.45),('Transition-sensitive manufacturing',40,.04,.25,.60)]:
 rate=severity*ndtr((ndtri(pd)+np.sqrt(cor)*ndtri(alpha))/np.sqrt(1-cor));sectors.append(dict(name=name,ead=ex,pd=pd,rho=cor,lgd=severity,marginal=rate,rc=ex*rate,el=ex*pd*severity))
# Panjer recursion: frequency p_n=(a+b/n)*p_(n-1); loss bands 1,2,5, severity probabilities .5,.3,.2.
f=np.zeros(6); f[[1,2,5]]=[.5,.3,.2]
def panjer(a,b,g0,nmax=500):
 g=np.zeros(nmax+1);g[0]=g0
 for n in range(1,nmax+1):g[n]=sum((a+b*j/n)*f[j]*g[n-j] for j in range(1,min(n,5)+1))
 assert abs(g.sum()-1)<1e-10
 return g
pg=panjer(0,4,np.exp(-4));ng=panjer(.5,1.5,1/16)
# Second independent gamma sector: shape 9, scale 2/9 -> mean frequency 2, same severity bands.
sec=panjer(2/11,8*2/11,(9/11)**9);agg=np.convolve(ng,sec)
def dist_stats(g):
 loss=np.arange(len(g))*.1; cum=np.cumsum(g);q=loss[np.searchsorted(cum,alpha)];return dict(el=float(g@loss),var=float(q),es=float((g[loss>q]@loss[loss>q]+q*(cum[np.searchsorted(cum,alpha)]-alpha))/(1-alpha)))
out=dict(seed=20261010,scenarios=N,gaussian=base,lgd_independent=random_ind,lgd_dependent=random_dep,t_copula=tc,t_tail_dependence=float(taildep),analytic_es=analytic_es,sectors=sectors,sector_var=sum(x['rc'] for x in sectors),creditrisk_poisson=dist_stats(pg),creditrisk_gamma=dist_stats(ng),creditrisk_two_sectors=dist_stats(agg))
Path('credit-var-advanced-results.json').write_text(json.dumps(out,indent=2));print(json.dumps(out,indent=2))

# Additional checks of the original article seed and a tail-event importance sampler.
r=out
a=.999;n=1000000;rng=np.random.default_rng(20261003);z=rng.standard_normal(n);p=ndtr((ndtri(.02)-np.sqrt(.15)*z)/np.sqrt(.85));l=.045*rng.binomial(1000,p);s=np.sort(l);r['original_mc']={'el':float(l.mean()),'var':float(s[998999]),'ci95':[float(s[998937]),float(s[999061])]}
rng=np.random.default_rng(20261011);z=rng.normal(-2,1,100000);p=ndtr((ndtri(.02)-np.sqrt(.15)*z)/np.sqrt(.85));d=rng.binomial(1000,p);w=np.exp(2*z+2);weighted=w*(d>178);r['importance_sampling']={'scenarios':100000,'proposal_mean':-2,'tail_prob':float(weighted.mean()),'standard_error':float(weighted.std(ddof=1)/np.sqrt(len(w))),'ordinary_theoretical_se':float(np.sqrt(.0009945461*(1-.0009945461)/100000))}
Path('credit-var-advanced-results.json').write_text(json.dumps(r,indent=2));

r["creditrisk_grid"]={"loss_unit_crore":0.1,"poisson_probabilities":pg.tolist(),"gamma_probabilities":ng.tolist(),"two_sector_probabilities":agg.tolist()}
Path("credit-var-advanced-results.json").write_text(json.dumps(r,indent=2))

# Deterministic checks for accounting, counterparty and stress illustrations.
def discrete_var(losses, probabilities, confidence):
 order=np.argsort(losses);return float(np.asarray(losses)[order][np.searchsorted(np.cumsum(np.asarray(probabilities)[order]),confidence)])
r['ecl']={'twelve_month':45*.01/1.05,'lifetime':45*(.01/1.05+.0198/1.05**2+.029106/1.05**3)}
r['macro_pd']=float(1/(1+np.exp(-(-3.8918202981106265+1.6))))
r['climate_sector_var']=float(40*.65*ndtr((ndtri(.06)+.5*ndtri(.999))/np.sqrt(.75)))
r['counterparty']={'independent_el':.01*6+.01*18,'wrong_way_el':.005*6+.015*18,'independent_var99':discrete_var([0,6,18],[.98,.01,.01],.99),'wrong_way_var99':discrete_var([0,6,18],[.98,.005,.015],.99)}
assert abs(r['ecl']['lifetime']-2.3681632653061224)<1e-12
assert abs(sum(x['rc'] for x in r['sectors'])-r['sector_var'])<1e-12
assert r['original_mc']['ci95'][0]<=8.01<=r['original_mc']['ci95'][1]
Path('credit-var-advanced-results.json').write_text(json.dumps(r,indent=2))