from pathlib import Path
from datetime import datetime, timezone
import json, hashlib
import numpy as np
import pandas as pd
import pandas_market_calendars as mcal

OUT=Path(__file__).resolve().parent
SRC=Path('C:/Users/Administrator/Documents/Codex/2026-11-04/windows-codex-ai-1-2-10/outputs/政策年度基调与资产回测/prices.csv')
def sha(p): return hashlib.sha256(p.read_bytes()).hexdigest()
def save(n,x): (OUT/n).write_text(json.dumps(x,ensure_ascii=False,indent=2,allow_nan=False),encoding='utf-8')
cfg=json.loads((OUT/'风险预算实验设定.json').read_text(encoding='utf-8'))
assert sha(SRC)==cfg['source_sha256']
assets=cfg['assets']; raw=pd.read_csv(SRC); raw=raw[raw.symbol.isin(assets)]
assert not raw.duplicated(['date','symbol']).any()
p=raw.pivot(index='date',columns='symbol',values='adjusted_close')[assets].loc['2007-05-30':'2026-11-04']
assert p.notna().all().all() and (p>0).all().all()
cal=mcal.get_calendar('NYSE').valid_days(start_date=p.index[0],end_date=p.index[-1]).strftime('%Y-%m-%d').tolist()
assert p.index.tolist()==cal
returns=p.pct_change(fill_method=None).to_numpy()
baseline=np.array(cfg['base_weights']);cashunit=np.array([0.,0.,0.,1.])

def erc(cov):
    # Coordinate minimization of 0.5*x'C*x - sum(log(x))/3.
    c=cov/np.mean(np.diag(cov));x=1/np.sqrt(np.diag(c));x/=x.sum()
    for it in range(20000):
        for j in range(3):
            a=c[j,j];q=float(c[j]@x-a*x[j]);x[j]=2/3/(q+np.sqrt(q*q+4*a/3))
        rc=x*(c@x);gap=float(np.max(np.abs(rc-1/3)))
        if gap<1e-11:break
    else: raise AssertionError('ERC convergence')
    w=x/x.sum();shares=w*(cov@w)/(w@cov@w)
    assert np.max(abs(shares-1/3))<1e-9
    return w,shares,it+1

def vol(w,c):return float(np.sqrt(w@c@w*252))
def shrink(w,c,target):
    if abs(vol(w,c)-target)<1e-12:return w.copy(),1.
    # Full four-asset covariance, including covariance with BIL.
    d=w-cashunit;a=float(d@c@d);b=float(2*cashunit@c@d);z=float(cashunit@c@cashunit-target**2/252)
    roots=np.roots([a,b,z]);ks=[float(k.real) for k in roots if abs(k.imag)<1e-10 and -1e-10<=k.real<=1+1e-10]
    assert ks,'No reachable risk target';k=max(ks);nw=cashunit+k*d
    assert abs(vol(nw,c)-target)<1e-9
    return nw,k

def signal(i,rr=returns):
    window=rr[i-251:i+1];assert window.shape==(252,4) and np.isfinite(window).all()
    c=np.cov(window,rowvar=False,ddof=1);assert np.linalg.eigvalsh(c).min()>0
    core,shares,it=erc(c[:3,:3]);w=np.r_[.75*core,.25];capped=False
    if w[0]>.5:
        capped=True;w[1:3]*=(.75-.5)/w[1:3].sum();w[0]=.5
    target=min(vol(w,c),vol(baseline,c));b,kb=shrink(baseline,c,target);e,ke=shrink(w,c,target)
    ws={'nominal_fixed':baseline.copy(),'nominal_erc':w,'matched_fixed':b,'matched_erc':e}
    for z in ws.values():
        assert abs(z.sum()-1)<1e-10 and z.min()>=-1e-10 and z[3]>=.25-1e-10 and z[0]<=.5+1e-10
    rc=w[:3]*(c[:3,:3]@w[:3])/(w[:3]@c[:3,:3]@w[:3])
    return ws,{'signal':p.index[i],'execution':p.index[i+1],'history_start':p.index[i-252],
      'covariance':c.tolist(),'target_annual_volatility':target,'nominal_vols':{n:vol(v,c) for n,v in ws.items()},
      'risky_core_RC_shares':rc.tolist(),'stock_cap_applied':capped,'solver_iterations':it,
      'scales':{'fixed':kb,'erc':ke}}

# Solver checks do not use market return performance to choose parameters.
cw,_,_=erc(np.diag([.04,.01,.09]));expected=1/np.array([.2,.1,.3]);expected/=expected.sum()
assert np.allclose(cw,expected,atol=1e-9)
cw,_,_=erc(np.full((3,3),.005)+np.eye(3)*.02);assert np.allclose(cw,1/3,atol=1e-9)
signals=[i for i in range(252,len(p)-1) if p.index[i][:7]!=p.index[i+1][:7]]
targets={};logs=[]
for i in signals:
    w,log=signal(i);log['weights']={n:z.tolist() for n,z in w.items()};targets[i+1]=w;logs.append(log)
assert signals
for i in [signals[0],signals[len(signals)//2],signals[-1]]:
    changed=returns.copy();changed[i+1:]+=10
    ws,_=signal(i,changed)
    for n,z in ws.items():assert np.array_equal(z,targets[i+1][n])

def trade(h,cash,w,c):
    nav=float(h.sum()+cash);lo,hi=nav/(1+2*c),nav
    for _ in range(55):
        mid=(lo+hi)/2
        if mid+c*np.abs(mid*w-h).sum()>nav:hi=mid
        else:lo=mid
    after=(lo+hi)/2;notional=float(np.abs(after*w-h).sum());fee=c*notional
    assert abs(after+fee-nav)<max(1e-7,nav*1e-12)
    return after,notional,fee
assert abs(trade(np.zeros(2),100.,np.array([.5,.5]),.0005)[0]-100/1.0005)<1e-10
assert abs(trade(np.array([100.,25.]),0.,np.array([.5,.5]),.001)[0]-124.925)<1e-10

def metrics(values,dates):
    v=np.array(values);dd=1-v/np.maximum.accumulate(v);trough=int(dd.argmax());peak=int(np.argmax(v[:trough+1]))
    recovered=np.where(v[trough:]>=v[peak])[0];recovery=trough+int(recovered[0]) if len(recovered) else None
    years=(pd.Timestamp(dates[-1])-pd.Timestamp(dates[0])).days/365.25
    episodes=[];high=v[0];highi=0;under=False
    for j in range(1,len(v)):
        if v[j]>=high:
            if under:episodes.append({'peak':dates[highi],'recovery':dates[j],'calendar_days':(pd.Timestamp(dates[j])-pd.Timestamp(dates[highi])).days,'censored':False})
            high=float(v[j]);highi=j;under=False
        else:under=True
    if under:episodes.append({'peak':dates[highi],'recovery':None,'calendar_days':(pd.Timestamp(dates[-1])-pd.Timestamp(dates[highi])).days,'censored':True})
    annual=[];prev=float(v[0]);lastyear=dates[0][:4]
    for year in sorted(set(d[:4] for d in dates[1:])):
        ids=[j for j,d in enumerate(dates) if d[:4]==year];end=ids[-1]
        annual.append({'year':int(year),'return':float(v[end]/prev-1),'partial':year in [dates[0][:4],dates[-1][:4]]});prev=float(v[end])
    rolling={}
    idx=pd.DatetimeIndex(dates);series=pd.Series(v,index=idx)
    # Month-end endpoints only. Closest available prior date to the calendar-year boundary.
    ends=[j for j in range(len(dates)-1) if dates[j][:7]!=dates[j+1][:7]]+[len(dates)-1]
    for horizon in [3,5,10,20]:
        rows=[]
        for j in ends:
            cutoff=idx[j]-pd.DateOffset(years=horizon);k=int(idx.searchsorted(cutoff,side='right')-1)
            if k<0:continue
            yy=(idx[j]-idx[k]).days/365.25
            rows.append({'start':dates[k],'end':dates[j],'cagr':float((v[j]/v[k])**(1/yy)-1)})
        rolling[str(horizon)]={'observations':len(rows),'overlapping_not_independent':True,
         'min':min(rows,key=lambda r:r['cagr']) if rows else None,'max':max(rows,key=lambda r:r['cagr']) if rows else None,
         'median':float(np.median([r['cagr'] for r in rows])) if rows else None}
    return {'cagr':float((v[-1]/v[0])**(1/years)-1),'cumulative_return':float(v[-1]/v[0]-1),
      'annual_volatility':float(np.std(v[2:]/v[1:-1]-1,ddof=1)*np.sqrt(252)),
      'max_drawdown':float(dd[trough]),'max_drawdown_peak':dates[peak],'max_drawdown_trough':dates[trough],
      'max_drawdown_recovered':dates[recovery] if recovery is not None else None,
      'peak_to_recovery_calendar_days':(pd.Timestamp(dates[recovery])-pd.Timestamp(dates[peak])).days if recovery is not None else None,
      'max_drawdown_elapsed_to_end_days':(pd.Timestamp(dates[-1])-pd.Timestamp(dates[peak])).days,
      'longest_underwater':max(episodes,key=lambda z:z['calendar_days']) if episodes else None,'years':years,'annual':annual,'rolling_years':rolling}

start=min(targets);dates=[p.index[start]]+p.index[start:].tolist();results=[];paths={};trades=[]
for bps in cfg['costs_one_way_bps']:
    for name in next(iter(targets.values())):
        h=np.zeros(4);cash=100000.;navs=[cash];cost=0.;turnover=0.;orders=[];dailyw=[]
        for i in range(start,len(p)):
            if i>start:h*=1+returns[i]
            nav=float(h.sum()+cash)
            if i in targets:
                w=targets[i][name];after,notional,fee=trade(h,cash,w,bps/10000)
                h=after*w;cash=0.;cost+=fee;turnover+=notional/nav
                if bps==5:orders.append({'signal':p.index[i-1],'execution':p.index[i],'weights':w.tolist(),'fees':fee,'traded_notional':notional})
            navs.append(float(h.sum()+cash));dailyw.append((h/(h.sum()+cash)).tolist())
        m=metrics(navs,dates);m.update({'strategy':name,'one_way_bps':bps,'fees_USD':cost,'traded_notional_over_nav_sum':turnover,
          'annual_two_way_turnover':turnover/m['years'],'annual_one_way_equivalent_turnover':turnover/(2*m['years']),
          'rebalance_events':len(targets),'average_daily_weights':dict(zip(assets,np.mean(dailyw,axis=0).tolist()))})
        results.append(m);paths[name+'_'+str(bps)]={'dates':p.index[start:].tolist(),'nav':navs[1:]}
        if bps==5:trades.append({'strategy':name,'orders':orders})
for name in next(iter(targets.values())):
    r=[x for x in results if x['strategy']==name];assert r[0]['cumulative_return']>=r[1]['cumulative_return']>=r[2]['cumulative_return']
assert sha(SRC)==cfg['source_sha256']
out={'created_at':datetime.now(timezone.utc).isoformat(),'settings':cfg,'settings_sha256':sha(OUT/'风险预算实验设定.json'),
 'source_path':str(SRC),'source_sha256':sha(SRC),'price_start':p.index[0],'first_signal':logs[0]['signal'],'first_trade':p.index[start],
 'last_valuation':p.index[-1],'price_sessions':len(p),'valuation_sessions':len(p)-start,'signals':len(logs),'results':results,
 'verification':{'calendar_exact':True,'solver_closed_forms':True,'fee_closed_forms':True,'positive_definite_covariance_all':True,
 'future_perturbation_3_dates':True,'weight_constraints_all':True,'cost_monotonicity':True,'stock_cap_applied_count':sum(x['stock_cap_applied'] for x in logs),
 'max_matched_vol_gap':max(abs(x['nominal_vols']['matched_fixed']-x['nominal_vols']['matched_erc']) for x in logs),'source_unchanged':True},
 'AI_performance':None,'account_changed':False,'limitations':['Retrospective modern adjusted-close proxy, not PIT data-vendor archive','No personal tax or FX','Estimated equal volatility is not equal crisis loss','No 20-year return window','No forecast of 8% or 10%','C6 untouched; target-weight execution differs from frozen-quantity C6']}
save('风险预算实验结果.json',out);save('风险预算实验路径.json',{'assets':assets,'signals':logs,'paths':paths,'trades_5bps':trades})
print(json.dumps({'first_trade':out['first_trade'],'end':out['last_valuation'],'verification':out['verification'],
 'results':[{k:v for k,v in r.items() if k not in ['annual','rolling_years']} for r in results if r['one_way_bps']==5]},ensure_ascii=False,indent=2))