"""PSIM-RN-05: public tail-risk mathematics with invented inputs only.

Python 3.10+, NumPy and SciPy. Run: python tail_risk_v3.py
Prints JSON. Reads no input files, contacts no services and loads no market data.
All losses are additive shocks in arbitrary units over one unspecified horizon.
These are numerical teaching models, not calibrated security-price processes.
"""
import json
import math
import platform
import numpy as np
import scipy
from scipy import sparse
from scipy.integrate import quad
from scipy.optimize import brentq, linprog
from scipy.stats import norm, t


def empirical_risk(losses, alpha, probabilities=None):
    """Left quantile and exact ES of a finite discrete loss distribution."""
    x = np.asarray(losses, dtype=float)
    if x.ndim != 1 or not x.size or not np.isfinite(x).all() or not 0 < alpha < 1:
        raise ValueError('finite nonempty loss vector and 0 < alpha < 1 required')
    p = np.full(x.size, 1/x.size) if probabilities is None else np.asarray(probabilities, dtype=float)
    if p.shape != x.shape or not np.isfinite(p).all() or (p < 0).any() or not np.isclose(p.sum(), 1, atol=1e-12, rtol=0):
        raise ValueError('probabilities must be nonnegative and sum to one')
    p = p / p.sum()
    order = np.argsort(x)
    index = min(int(np.searchsorted(np.cumsum(p[order]), alpha, side='left')), x.size-1)
    q = float(x[order[index]])
    es = q + float(p @ np.maximum(x-q, 0)) / (1-alpha)
    return q, es


def standard_tail(alpha, nu=5):
    """Analytic VaR and ES for normal and variance-one Student-t losses."""
    if not 0 < alpha < 1 or nu <= 2:
        raise ValueError('require 0 < alpha < 1 and nu > 2')
    normal_q = float(norm.ppf(alpha))
    student_q = float(t.ppf(alpha, nu))
    scale = math.sqrt((nu-2)/nu)
    return {'normal': {'var':normal_q, 'es':float(norm.pdf(normal_q)/(1-alpha))},
        'student_t': {'var':scale*student_q,
            'es':float(scale*(nu+student_q**2)/(nu-1)*t.pdf(student_q,nu)/(1-alpha))}}


def coexceedance(level, rho, nu=None):
    """P(Y>its level-quantile | X>its level-quantile), by quadrature."""
    if not 0 < level < 1 or not -1 < rho < 1:
        raise ValueError('require a probability and a nonsingular correlation')
    if nu is None:
        q = norm.ppf(level)
        def integrand(x):
            return norm.pdf(x)*norm.sf((q-rho*x)/math.sqrt(1-rho*rho))
    else:
        if nu <= 0:
            raise ValueError('positive degrees of freedom required')
        q = t.ppf(level,nu)
        def integrand(x):
            scale = math.sqrt((nu+x*x)*(1-rho*rho)/(nu+1))
            return t.pdf(x,nu)*t.sf((q-rho*x)/scale,nu+1)
    value, error = quad(integrand, q, np.inf, epsabs=1e-12, epsrel=1e-9)
    return float(value/(1-level)), float(error/(1-level))


def mixture_parameters():
    # Constructed law shared by all scenarios; no parameter is fitted.
    return 0.04, 0.20, np.array([0.35,0.95,0.65]), np.array([4.0,0.0,1.5])


def mixture_draws(count, seed):
    p, common_scale, residual_scale, jump = mixture_parameters()
    rng = np.random.Generator(np.random.PCG64(seed))
    common = rng.standard_normal((count,1))
    residual = rng.standard_normal((count,3))
    event = (rng.random((count,1)) < p).astype(float)
    return common_scale*common + residual*residual_scale + (event-p)*jump


def mixture_risk(w, alpha):
    """Population risk of a fixed portfolio under the known two-normal mixture."""
    p, common_scale, residual_scale, jump = mixture_parameters()
    exposure = float(w @ jump)
    means = np.array([-p*exposure, (1-p)*exposure])
    probs = np.array([1-p,p])
    sd = math.sqrt((common_scale*w.sum())**2 + float((w*residual_scale) @ (w*residual_scale)))
    def cdf(x):
        return float(probs @ norm.cdf((x-means)/sd))
    q = brentq(lambda x: cdf(x)-alpha, float(means.min()-12*sd), float(means.max()+12*sd), xtol=1e-12)
    z = (q-means)/sd
    es = float(probs @ (means*norm.sf(z) + sd*norm.pdf(z))/(1-alpha))
    variance = sd*sd + p*(1-p)*exposure*exposure
    return {'var':float(q),'es':es,'standard_deviation':math.sqrt(variance),
            'cdf_at_var_error':abs(cdf(q)-alpha)}


def solve_cvar(losses, alpha, cap):
    """Public scenario CVaR LP for fixed linear losses and long-only weights."""
    m,n = losses.shape
    if not np.isfinite(losses).all() or not 0 < alpha < 1 or cap*n < 1:
        raise ValueError('invalid scenarios, confidence level or infeasible cap')
    objective = np.r_[np.zeros(n), 1.0, np.full(m,1/((1-alpha)*m))]
    aub = sparse.hstack([sparse.csc_matrix(losses), -sparse.csc_matrix(np.ones((m,1))),
                         -sparse.eye(m,format='csc')],format='csc')
    aeq = sparse.csc_matrix((np.ones(n), (np.zeros(n,dtype=int),np.arange(n))),shape=(1,n+1+m))
    result = linprog(objective,A_ub=aub,b_ub=np.zeros(m),A_eq=aeq,b_eq=[1.0],
        bounds=[(0,cap)]*n+[(None,None)]+[(0,None)]*m, method='highs',
        options={'dual_feasibility_tolerance':1e-8,'primal_feasibility_tolerance':1e-8})
    if not result.success or result.status != 0:
        raise ArithmeticError('CVaR LP did not solve successfully')
    w = result.x[:n]
    stationarity = objective - aub.T@result.ineqlin.marginals - aeq.T@result.eqlin.marginals - result.lower.marginals - result.upper.marginals
    primal = max(float(np.max(aub@result.x)),float(abs(w.sum()-1)),float(-w.min()),
                 float(w.max()-cap),float(-result.x[n+1:].min()),0.0)
    dual_objective = float(result.eqlin.marginals[0] + cap*result.upper.marginals[:n].sum())
    complementarity = max(float(np.max(np.abs(result.ineqlin.residual*result.ineqlin.marginals))),
        float(np.max(np.abs(result.x[:n]*result.lower.marginals[:n]))),
        float(np.max(np.abs(result.x[n+1:]*result.lower.marginals[n+1:]))),
        float(np.max(np.abs((cap-w)*result.upper.marginals[:n]))))
    diagnostics = {'status':int(result.status),'objective':float(result.fun),
        'primal_max_violation':primal,'stationarity_max_residual':float(np.max(np.abs(stationarity))),
        'absolute_duality_gap':abs(float(result.fun)-dual_objective),
        'complementarity_max_residual':complementarity,
        'dual_sign_violation':max(float(result.ineqlin.marginals.max()),
            float(-result.lower.marginals.min()),float(result.upper.marginals.max()),0.0),
        'variables':n+1+m,'scenario_constraints':m,'scenario_matrix_nonzeros':int(aub.nnz)}
    return w, diagnostics


def large_universe(alpha, checks):
    n,k,m,batch = 10_000,50,100_000,2_000
    rng = np.random.Generator(np.random.PCG64(20260911))
    a = rng.normal(0,0.10/math.sqrt(k),size=(n,k))
    a[:,0] = np.linspace(0.03,0.30,n)
    d = rng.uniform(0.10,0.30,n)**2
    w = np.zeros((n,3))
    w[:,0] = 1/n
    w[:1000,1] = 1/1000
    w[:,2] = np.linspace(1.0,2.0,n)
    w[:,2] /= w[:,2].sum()
    exposures = a.T @ w
    projected = exposures.T@exposures + w.T@(d[:,None]*w)
    root = np.linalg.cholesky(projected)
    direct = w.T@(a@(a.T@w) + d[:,None]*w)
    checks['projected_covariance_matches_operator'] = bool(np.allclose(projected,direct,atol=1e-12,rtol=1e-12))
    dense_small = a[:128]@a[:128].T + np.diag(d[:128])
    small_projection = (a[:128].T@w[:128]).T@(a[:128].T@w[:128]) + w[:128].T@(d[:128,None]*w[:128])
    checks['projected_covariance_matches_dense_128'] = bool(np.allclose(small_projection,w[:128].T@dense_small@w[:128],atol=1e-12,rtol=1e-12))
    checks['all_fixed_portfolios_fully_invested'] = bool(np.allclose(w.sum(axis=0),1,atol=1e-12,rtol=0))
    rng_z = np.random.Generator(np.random.PCG64(20260912))
    rng_u = np.random.Generator(np.random.PCG64(20260913))
    losses = np.empty((m,3))
    for start in range(0,m,batch):
        count = min(batch,m-start)
        normal = rng_z.standard_normal((count,3)) @ root.T
        losses[start:start+count] = normal*np.sqrt(3/rng_u.chisquare(5,count))[:,None]
    constant = standard_tail(alpha)['student_t']['es']
    var_constant = standard_tail(alpha)['student_t']['var']
    rows = []
    for i,name in enumerate(('equal_weight','lower_common_exposure_1000_equal_weight','linear_tilt')):
        sd = math.sqrt(float(projected[i,i]))
        q,es = empirical_risk(losses[:,i],alpha)
        # The exact quantile supplies an independent integration diagnostic.
        score = var_constant*sd + np.maximum(losses[:,i]-var_constant*sd,0)/(1-alpha)
        se = float(score.std(ddof=1)/math.sqrt(m))
        exact_es = constant*sd
        checks['large_es_within_six_sampling_errors_'+name] = abs(es-exact_es) < 6*se
        rows.append({'portfolio':name,'model_standard_deviation':sd,'sample_var':q,'sample_es':es,
                     'analytic_es':exact_es,'es_sampling_error_diagnostic':se})
    checks['large_scenarios_finite'] = bool(np.isfinite(losses).all())
    return {'assets':n,'factors':k,'fixed_portfolios':3,'scenarios':m,'batch_size':batch,
        'seeds':{'constructed_covariance':20260911,'normal_draws':20260912,'common_scale':20260913},
        'storage_bytes':{'hypothetical_dense_asset_scenarios':8*n*m,'factor_and_diagonal':a.nbytes+d.nbytes,
            'fixed_weights':w.nbytes,'retained_portfolio_losses':losses.nbytes,'one_portfolio_loss_block':batch*3*8},
        'projected_covariance':projected.tolist(),'results':rows}


def reproduce():
    alpha,nu,rho = 0.99,5,0.5
    checks = {}
    tails = standard_tail(alpha,nu)
    scale = math.sqrt((nu-2)/nu)
    integrated_es,error = quad(lambda x: x*t.pdf(x/scale,nu)/scale,tails['student_t']['var'],np.inf,epsabs=1e-11)
    checks['student_es_matches_density_integral'] = abs(integrated_es/(1-alpha)-tails['student_t']['es']) < 1e-9
    checks['scaled_student_variance_one'] = math.isclose(float(t.var(nu))*scale*scale,1,rel_tol=1e-12)
    checks['normal_var_tail_probability'] = abs(float(norm.sf(tails['normal']['var']))-(1-alpha)) < 1e-12
    checks['student_var_tail_probability'] = abs(float(t.sf(tails['student_t']['var']/scale,nu))-(1-alpha)) < 1e-10
    dependence = []
    for level in (0.95,0.99,0.999):
        g,ge = coexceedance(level,rho)
        s,se = coexceedance(level,rho,nu)
        checks['coexceedance_probability_bounds_'+str(level)] = 0 <= g <= 1 and 0 <= s <= 1
        checks['coexceedance_quadrature_error_'+str(level)] = max(ge,se) < 1e-7
        dependence.append({'quantile_level':level,'gaussian':g,'student_t':s,'max_quadrature_error':max(ge,se)})
    tail_coefficient = float(2*t.cdf(-math.sqrt((nu+1)*(1-rho)/(1+rho)),nu+1))
    q,es = empirical_risk([0,10],0.975,[0.99,0.01])
    checks['discrete_atom_var_zero_es_four'] = q == 0 and math.isclose(es,4,abs_tol=1e-12)
    base = np.array([-2.0,-0.5,0.2,0.2,1.0,4.0])
    _,base_es = empirical_risk(base,0.7)
    _,shift_es = empirical_risk(base+3,0.7)
    _,scale_es = empirical_risk(2*base,0.7)
    checks['empirical_es_cash_translation'] = math.isclose(shift_es,base_es+3,abs_tol=1e-12)
    checks['empirical_es_positive_homogeneity'] = math.isclose(scale_es,2*base_es,abs_tol=1e-12)
    checks['empirical_es_fractional_tail_mass'] = math.isclose(base_es,(4+0.8)/1.8,abs_tol=1e-12)
    p,common,residual,jump = mixture_parameters()
    covariance = common**2*np.ones((3,3)) + np.diag(residual**2) + p*(1-p)*np.outer(jump,jump)
    minimum = np.linalg.solve(covariance,np.ones(3))
    minimum /= minimum.sum()
    cap = 0.8
    checks['population_minimum_variance_satisfies_caps'] = bool(minimum.min() >= 0 and minimum.max() <= cap)
    train = mixture_draws(20_000,20260908)
    evaluation = mixture_draws(100_000,20260909)
    optimized,diagnostics = solve_cvar(train,alpha,cap)
    checks['lp_primal_feasible'] = diagnostics['primal_max_violation'] < 1e-7
    checks['lp_dual_stationarity'] = diagnostics['stationarity_max_residual'] < 1e-7
    checks['lp_dual_signs'] = diagnostics['dual_sign_violation'] < 1e-7
    checks['lp_duality_gap'] = diagnostics['absolute_duality_gap'] < 1e-7
    checks['lp_complementarity'] = diagnostics['complementarity_max_residual'] < 1e-7
    rows = []
    for name,w in (('equal_weight',np.ones(3)/3),('minimum_variance',minimum),('minimum_sample_es',optimized)):
        _,training_es = empirical_risk(train@w,alpha)
        population = mixture_risk(w,alpha)
        draws = evaluation@w
        evaluation_var,evaluation_es = empirical_risk(draws,alpha)
        score = population['var']+np.maximum(draws-population['var'],0)/(1-alpha)
        se = float(score.std(ddof=1)/math.sqrt(len(draws)))
        checks['mixture_es_within_six_sampling_errors_'+name] = abs(evaluation_es-population['es']) < 6*se
        checks['mixture_quantile_root_'+name] = population['cdf_at_var_error'] < 1e-10
        rows.append({'portfolio':name,'weights':w.tolist(),'training_es':training_es,
            'evaluation_var':evaluation_var,'evaluation_es':evaluation_es,'population':population,
            'es_sampling_error_diagnostic':se})
    checks['lp_objective_matches_discrete_es'] = abs(rows[2]['training_es']-diagnostics['objective']) < 1e-8
    checks['lp_beats_feasible_training_comparators'] = rows[2]['training_es'] <= min(rows[0]['training_es'],rows[1]['training_es'])+1e-8
    large = large_universe(alpha,checks)
    if not all(checks.values()):
        raise ArithmeticError('Numerical checks failed: '+', '.join(k for k,v in checks.items() if not v))
    return {'paper':'PSIM-RN-05','classification':'synthetic_teaching_example',
        'scope':'Public methods, invented constants and seeded draws. No market data, proprietary calibration, signals, holdings or internal performance.',
        'runtime':{'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__},
        'inputs':{'alpha':alpha,'student_degrees_of_freedom':nu,'pairwise_correlation':rho,
            'loss_units':'arbitrary additive shocks; one unspecified horizon; no annualisation'},
        'variance_one_models':tails,'student_es_increase_percent':100*(tails['student_t']['es']/tails['normal']['es']-1),
        'dependence':{'finite_quantile_coexceedance':dependence,'gaussian_tail_coefficient':0.0,'student_tail_coefficient':tail_coefficient},
        'discrete_example':{'losses':[0,10],'probabilities':[0.99,0.01],'alpha':0.975,'var':q,'es':es},
        'mixture':{'event_probability':p,'common_normal_loading':common,'residual_standard_deviations':residual.tolist(),
            'centred_event_loadings':jump.tolist(),'population_covariance':covariance.tolist(),'training_scenarios':len(train),
            'evaluation_scenarios':len(evaluation),'training_seed':20260908,'evaluation_seed':20260909,
            'training_tail_mass_scenarios': int(round((1-alpha)*len(train))), 'evaluation_tail_mass_scenarios': int(round((1-alpha)*len(evaluation))), 'weight_cap_nonbinding': bool(optimized.max() < cap - 1e-7),
            'weight_cap':cap,'optimization':diagnostics,'portfolios':rows},
        'large_universe':large,'checks':{key:bool(value) for key,value in checks.items()}}


if __name__ == '__main__':
    print(json.dumps(reproduce(),indent=2,sort_keys=True,allow_nan=False))