"""Reproduce PSIM's interactive methods from public, invented inputs.

Python 3.10+, NumPy and SciPy. Download figures_data_v1.json from the
paper's figure-reproduction section, then run:
  python interactive_methods_20260915_v1.py --source figures_data_v1.json
This writes interactive_methods_20260915_v1.json beside the command.
"""
from pathlib import Path
import argparse,json,math
import numpy as np
from scipy.stats import norm,t
def holm(raw):
    order=sorted(range(len(raw)),key=lambda i:raw[i])
    adjusted=[0.0]*len(raw)
    running=0.0
    for rank,i in enumerate(order):
        running=max(running,(len(raw)-rank)*raw[i])
        adjusted[i]=min(1.0,running)
    return adjusted

def charges(bps,fixed):
    return [100000*bps/10000+100*fixed,60000*bps/10000+20*fixed]

def drawdowns(linked=True):
    wealth=[1,1.1,1.04,1.04,1.0192,.988]
    peak=np.maximum.accumulate(wealth) if linked else [1,1.1,1.1,1.04,1.04,1.04]
    return [1-w/p for w,p in zip(wealth,peak)]

def portfolio_volatility(rho):
    w=np.array([.5,.5])
    covariance=np.array([[.04,.04*rho],[.04*rho,.04]])
    return float(np.sqrt(max(0,w@covariance@w)))

def apparent_h(noise_ratio,loglag):
    u=10**loglag
    return .5*u/(u+noise_ratio)

parser=argparse.ArgumentParser(description=__doc__)
parser.add_argument('--source',type=Path,default=Path('figures_data_v1.json'))
parser.add_argument('--output',type=Path,default=Path('interactive_methods_20260915_v1.json'))
args=parser.parse_args()
published=json.loads(args.source.read_bytes())
gbm=published['gbm']
times=np.array(gbm['time_years'])
prices=np.array(gbm['paths'])
brownian=(np.log(prices/100)-(.06-.5*.2**2)*times[:,None])/.2
x=np.linspace(0,14,701)
tails={'x':x.tolist(),'normal_survival':norm.sf(x).tolist(),'models':{}}
for df in [3,5,10,30]:
    scale=math.sqrt((df-2)/df)
    metrics={}
    for alpha in [.95,.99,.999]:
        z=norm.ppf(alpha); q=t.ppf(alpha,df)
        es=scale*(df+q*q)/(df-1)*t.pdf(q,df)/(1-alpha)
        metrics[str(alpha)]={'normal_var':float(z),'normal_es':float(norm.pdf(z)/(1-alpha)),'student_var':float(scale*q),'student_es':float(es)}
    tails['models'][str(df)]={'survival':t.sf(x/scale,df).tolist(),'metrics':metrics}
data={'scope':'Public mathematics; invented inputs. No market observations or PSIM performance.','source':'/examples/figures_data_v1.json','gbm':{'time':times.tolist(),'brownian':brownian.tolist(),'seed':20260909,'drift':.06,'initial':100},'tails':tails,'examples':{'holm_raw':[.012,.03,.08,.4],'holm_adjusted':holm([.012,.03,.08,.4]),'costs_at_5bp_zero_fixed':charges(5,0),'drawdown_linked':drawdowns(True),'drawdown_reset':drawdowns(False),'equal_weight_volatility_at_correlation_half':portfolio_volatility(.5),'apparent_h_at_noise_004_lag_001':apparent_h(.04,-2)}}
assert np.allclose(holm([.012,.03,.08,.4]),published['holm']['adjusted'])
assert np.allclose(drawdowns(True),published['drawdown']['linked_drawdown'])
assert np.allclose(100*np.exp((.06-.5*.2**2)*times[:,None]+.2*brownian),prices)
assert abs(tails['models']['5']['metrics']['0.99']['student_es']-published['tails']['student_es'])<1e-10
args.output.write_text(json.dumps(data,separators=(',',':'))+'\n',encoding='utf-8')
print('Prepared fixed GBM shocks and analytic variance-one tail curves for four degrees of freedom and three confidence levels.')