#!/usr/bin/env python3
"""
Standardised 0-25% gold allocation backtest. FROZEN METHODOLOGY (approved 2026-09-02).

  period          1968-05 .. 2024-06 monthly (674 return months), common to all series
  allocations     0, 5, 10, 15, 20, 25% gold
  non-gold split  60/40 PRESERVED WITHIN the non-gold portion (10% gold -> 54/36/10)
  rebalancing     annual, calendar year-end (primary). No-rebalancing = sensitivity only.
  equity          Shiller S&P composite price + dividend/12 reinvested monthly
  bonds           independently constructed FRED GS10 constant-maturity total return
  gold            World Bank Pink Sheet monthly average, price return only
  inflation       Shiller CPI (= FRED CPIAUCNS, validated exact)
  costs           NONE modelled: no fees, taxes, transaction costs, storage or custody
  drawdown        measured from MONTHLY observations

Methodology frozen before results were seen. No parameter is tuned post hoc.
"""
import csv, os, math, json
from collections import OrderedDict
HERE=os.path.dirname(os.path.abspath(__file__))
ALLOC=[0.0,0.05,0.10,0.15,0.20,0.25]

rows=list(csv.reader(open(os.path.join(HERE,'normalised','monthly_series.csv'))))[1:]
M=[r[0] for r in rows]
GOLD={r[0]:float(r[1]) for r in rows}; P={r[0]:float(r[2]) for r in rows}
DIV={r[0]:float(r[3]) for r in rows}; CPI={r[0]:float(r[4]) for r in rows}
BOND={r[0]:float(r[1]) for r in list(csv.reader(open(os.path.join(HERE,'normalised','bond_monthly_returns.csv'))))[1:]}

R=OrderedDict()
for i in range(1,len(M)):
    a,b=M[i-1],M[i]
    R[b]=(( P[b]+DIV[b]/12.0)/P[a]-1.0, BOND[b], GOLD[b]/GOLD[a]-1.0, CPI[b]/CPI[a]-1.0)
K=list(R)

def run(gw, rebal='annual'):
    tgt=[0.60*(1-gw),0.40*(1-gw),gw]; w=list(tgt)
    nom=[]; real=[]; rv=1.0
    for k in K:
        eq,bd,gd,inf=R[k]
        st=sum(w)
        w=[w[0]*(1+eq), w[1]*(1+bd), w[2]*(1+gd)]
        tot=sum(w); nom.append(tot)
        rv*= (tot/st)/(1+inf); real.append(rv)
        if rebal=='annual' and k.endswith('-12'): w=[tot*t for t in tgt]
    return nom, real, tgt

def mets(series):
    n=len(series); yrs=n/12.0
    cagr=series[-1]**(1/yrs)-1
    mr=[series[0]-1]+[series[i]/series[i-1]-1 for i in range(1,n)]
    mu=sum(mr)/n; vol=math.sqrt(sum((x-mu)**2 for x in mr)/(n-1)*12)
    peak=-1e9; mdd=0.0
    for v in series:
        peak=max(peak,v); mdd=min(mdd, v/peak-1)
    return cagr,vol,mdd,series[-1],mr

def cal_years(series):
    out=OrderedDict(); prev=1.0; buf=OrderedDict()
    for i,k in enumerate(K): buf.setdefault(int(k[:4]),[]).append(series[i])
    prev=1.0
    for y,v in buf.items():
        if len(v)==12:
            out[y]=v[-1]/prev-1; prev=v[-1]
        else:
            prev=v[-1]
    return out

# risk-free proxy: 3m T-bill unavailable in our validated set; use realised inflation-free
# construction -> report Sharpe vs the 10y bond sleeve's own return is NOT standard.
# We use the GS10 yield as the risk-free proxy and STATE IT.
RF=OrderedDict()
for r in list(csv.reader(open(os.path.join(HERE,'data','GS10.csv'))))[1:]:
    if len(r)<2 or r[1] in ('.',''): continue
    y,m,_=r[0].split('-'); RF['%04d-%02d'%(int(y),int(m))]=float(r[1])/100.0

def sharpe(mr):
    ex=[mr[i]-RF[K[i]]/12.0 for i in range(len(K)) if K[i] in RF]
    n=len(ex); mu=sum(ex)/n
    sd=math.sqrt(sum((x-mu)**2 for x in ex)/(n-1))
    return (mu*12)/(sd*math.sqrt(12))

print('='*94)
print('  STANDARDISED GOLD ALLOCATION BACKTEST — RAW RESULTS')
print('  %s to %s   %d monthly returns   annual rebalancing   60/40 preserved in non-gold'%(K[0],K[-1],len(K)))
print('='*94)
res={}
print('\nNOMINAL')
print('%-6s %7s %7s %8s %8s %10s %9s %16s'%('gold','wts','CAGR','vol','maxDD','worst yr','Sharpe','$100k ->'))
for gw in ALLOC:
    nom,real,tgt=run(gw); c,v,d,e,mr=mets(nom)
    cy=cal_years(nom); wy=min(cy.values()); wyk=min(cy,key=cy.get)
    sh=sharpe(mr)
    res[gw]={'tgt':tgt,'cagr':c,'vol':v,'mdd':d,'end':e,'worst_year':wy,'worst_year_k':wyk,'sharpe':sh}
    print('%-6s %7s %6.2f%% %7.2f%% %7.1f%% %6.1f%%(%d) %9.3f %16s'
          %('%d%%'%(gw*100),'%d/%d'%(round(tgt[0]*100),round(tgt[1]*100)),c*100,v*100,d*100,wy*100,wyk,sh,'${:,.0f}'.format(e*100000)))

print('\nREAL (inflation-adjusted, Shiller CPI)')
print('%-6s %8s %8s %16s'%('gold','CAGR','maxDD','$100k ->'))
for gw in ALLOC:
    nom,real,tgt=run(gw); c,v,d,e,_=mets(real)
    res[gw].update({'real_cagr':c,'real_mdd':d,'real_end':e})
    print('%-6s %7.2f%% %7.1f%% %16s'%('%d%%'%(gw*100),c*100,d*100,'${:,.0f}'.format(e*100000)))

print('\nSENSITIVITY — NO REBALANCING (reported separately, not part of the primary result)')
print('%-6s %8s %8s %8s'%('gold','CAGR','vol','maxDD'))
for gw in ALLOC:
    nom,_,_=run(gw,rebal='none'); c,v,d,e,_=mets(nom)
    res[gw].update({'norebal_cagr':c,'norebal_vol':v,'norebal_mdd':d})
    print('%-6s %7.2f%% %7.2f%% %7.1f%%'%('%d%%'%(gw*100),c*100,v*100,d*100))

# rolling windows + drawdown-frequency vs 0% baseline
def rolling(series,w):
    out=[]
    for i in range(w,len(series)+1):
        s=series[i-w:i]; base=series[i-w-1] if i-w-1>=0 else 1.0
        n=[x/base for x in s]
        yrs=w/12.0; c=n[-1]**(1/yrs)-1
        peak=-1e9; mdd=0.0
        for v in n:
            peak=max(peak,v); mdd=min(mdd,v/peak-1)
        out.append((K[i-1],c,mdd))
    return out

print('\nROLLING WINDOWS')
base={}
for w,lab in ((120,'10-year'),(240,'20-year')):
    b=rolling(run(0.0)[0],w); base[w]=b
    print('\n%s windows (n=%d)'%(lab,len(b)))
    print('%-6s %9s %9s %9s %11s'%('gold','med CAGR','min','max','med maxDD'))
    for gw in ALLOC:
        r=rolling(run(gw)[0],w)
        cs=sorted(x[1] for x in r); ds=sorted(x[2] for x in r)
        med=lambda a: a[len(a)//2] if len(a)%2 else (a[len(a)//2-1]+a[len(a)//2])/2
        print('%-6s %8.2f%% %8.2f%% %8.2f%% %10.1f%%'%('%d%%'%(gw*100),med(cs)*100,cs[0]*100,cs[-1]*100,med(ds)*100))
        res[gw]['roll%d'%w]={'median_cagr':med(cs),'min':cs[0],'max':cs[-1],'median_mdd':med(ds)}

print('\nHOW OFTEN DID GOLD LOWER MAXIMUM DRAWDOWN vs THE 0% PORTFOLIO?')
for w,lab in ((120,'10-year'),(240,'20-year')):
    print('\n  %s rolling windows:'%lab)
    b=base[w]
    for gw in ALLOC[1:]:
        r=rolling(run(gw)[0],w)
        n=sum(1 for i in range(len(r)) if r[i][2] > b[i][2])   # less negative = lower drawdown
        print('    %3d%% gold: lower max drawdown in %d of %d windows (%.1f%%)'%(gw*100,n,len(r),n/len(r)*100))
        res[gw].setdefault('dd_beat',{})['%d'%w]=[n,len(r)]

json.dump({str(k):v for k,v in res.items()}, open(os.path.join(HERE,'normalised','results.json'),'w'), indent=2, default=float)
print('\nwrote normalised/results.json')