"""Reproduce the frozen 2026-09-30 research package. No credentials or network needed.
Python 3.9+, numpy, pandas, xlrd, matplotlib, Pillow.
"""
from pathlib import Path
import json, hashlib
import numpy as np
import pandas as pd
ROOT=Path(__file__).resolve().parent
D=ROOT/'data'
Y=365.2425
END=pd.Timestamp('2026-08-31'); NOW=pd.Timestamp('2026-09-30'); ORIGIN=pd.Timestamp('2009-01-03')
m=pd.read_csv(D/'m2.csv',parse_dates=['observation_date']).set_index('observation_date').M2SL
btc=pd.read_csv(D/'btc-price.csv',parse_dates=['time']).set_index('time').PriceUSD
spot=json.loads((D/'btc-spot.json').read_text())
P=float(spot['data']['amount'])
r=pd.read_csv(D/'nyu-annual.csv').set_index('year')
sp=(1+r.sp_return).cumprod()*100
gold=r.gold_usd
tr=pd.read_csv(D/'sptr.csv',parse_dates=['date']).set_index('date')['close']
sp_ytd=tr.loc[END]/tr.loc['2025-12-31']
sp_end=sp.loc[2025]*sp_ytd
gold_end=4563.0 # WGC August 2026 commentary, Table 1, rounded USD/troy oz.
def years(a,b):return (pd.Timestamp(b)-pd.Timestamp(a)).days/Y
def m2(date):return float(m.loc[pd.Timestamp(date).to_period('M').to_timestamp()])
def row(name,start,end,a0,a1,m0,m1,note=''):
 n=years(start,end);g=(a1/a0)**(1/n)-1;mg=(m1/m0)**(1/n)-1;adj=(1+g)/(1+mg)-1;wealth=(a1/a0)/(m1/m0)
 assert abs(adj-(wealth**(1/n)-1))<1e-12
 return dict(asset=name,start=str(pd.Timestamp(start).date()),end=str(pd.Timestamp(end).date()),years=n,start_value=float(a0),end_value=float(a1),m2_start_bn=float(m0),m2_end_bn=float(m1),nominal_cagr=g,m2_cagr=mg,adjusted_cagr=adj,nominal_multiple=a1/a0,adjusted_multiple=wealth,note=note)
rows=[row('Gold','1959-12-31',END,gold[1959],gold_end,m2('1959-12-31'),m2(END),'1959 annual-average gold proxy; December M2 monthly average; August 2026 WGC rounded price.'),row('S&P 500 total return','1959-12-31',END,sp[1959],sp_end,m2('1959-12-31'),m2(END),'Annual dividend reinvestment through 2025; Yahoo S&P total-return index splice in 2026.'),row('Bitcoin inception quote','2009-10-05',NOW,1/1309.03,P,m2('2009-10-05'),m2(END),'October 2009 monthly M2 proxy; latest August M2 carried to September 30; not observed September M2.'),row('Bitcoin inception aligned endpoint','2009-10-05',END,1/1309.03,btc[END],m2('2009-10-05'),m2(END),'October monthly M2 proxy; August endpoint month aligned.'),row('Bitcoin daily-series start','2010-07-18',END,btc.iloc[0],btc[END],m2('2010-07-18'),m2(END),'Daily price series; endpoint-month M2 averages.')]
pd.DataFrame(rows).to_csv(D/'historical-calculations.csv',index=False)
# Annual observations for common-start paths; points connect annual observations, not daily paths.
path=[]
for year in range(1959,2026):
 date=pd.Timestamp(f'{year}-12-31');path.append(dict(date=str(date.date()),gold=gold[year],sp_total=sp[year],m2_bn=m2(date),btc=btc.get(date,np.nan)))
path.append(dict(date=str(END.date()),gold=gold_end,sp_total=sp_end,m2_bn=m2(END),btc=btc[END]))
paths=pd.DataFrame(path);paths.to_csv(D/'annual-asset-levels.csv',index=False)
# OLS equal daily observation weights. Excludes the isolated 2009 quote.
x=np.log10((btc.index-ORIGIN).days.to_numpy());y=np.log10(btc.to_numpy());b,a=np.polyfit(x,y,1)
fitted=a+b*x;resid=y-fitted
assert len(btc)==len(pd.date_range(btc.index.min(),btc.index.max()))
def trend(date,bb=b,aa=a):return 10**aa*float((pd.Timestamp(date)-ORIGIN).days)**bb
m2_n=(len(m)-1)/12;mg=(m.iloc[-1]/m.iloc[0])**(1/m2_n)-1
forward=[]
for year in [2030,2035,2040,2050]:
 target=pd.Timestamp(f'{year}-12-31');n=years(NOW,target);future=trend(target)
 gn=(future/P)**(1/n)-1;tn=(future/trend(NOW))**(1/n)-1
 forward.append(dict(target=str(target.date()),years=n,days_since_origin=(target-ORIGIN).days,median_usd=future,market_nominal=gn,market_adjusted=(1+gn)/(1+mg)-1,trend_nominal=tn,trend_adjusted=(1+tn)/(1+mg)-1,m2_growth_assumed=mg,m2_multiple=(1+mg)**n))
pd.DataFrame(forward).to_csv(D/'forward-calculations.csv',index=False)
sensitivity=[]
for start in ['2010-07-18','2011-01-01','2013-01-01','2015-01-01']:
 ss=btc.loc[start:];xx=np.log10((ss.index-ORIGIN).days);bb,aa=np.polyfit(xx,np.log10(ss),1)
 sensitivity.append(dict(start=start,n=len(ss),a=float(aa),b=float(bb),current_trend=trend(NOW,bb,aa),median_2050=trend('2050-12-31',bb,aa)))
train=btc.loc[:'2019-12-31'];bb,aa=np.polyfit(np.log10((train.index-ORIGIN).days),np.log10(train),1)
common=[]
for key,label in [('gold','Gold'),('sp_total','S&P 500 total return'),('btc','Bitcoin')]:
 base=paths[paths.date=='2010-12-31'].iloc[0];last=paths.iloc[-1]
 common.append(row(label,'2010-12-31',END,base[key],last[key],base.m2_bn,last.m2_bn))
start_sens=[]
for year in [1971,1980,2000,2010]:
 for asset,prices,last in [('Gold',gold,gold_end),('S&P 500 total return',sp,sp_end)]:start_sens.append(row(asset,f'{year}-12-31',END,prices[year],last,m2(f'{year}-12-31'),m2(END)))
pd.DataFrame(common).to_csv(D/'common-period-calculations.csv',index=False)
pd.DataFrame(start_sens).to_csv(D/'start-date-sensitivity.csv',index=False)
pd.DataFrame(sensitivity).to_csv(D/'model-sensitivity.csv',index=False)
summary=dict(asof='2026-09-30',spot=spot,m2_long=dict(start='1959-01',end='2026-08',start_bn=float(m.iloc[0]),end_bn=float(m.iloc[-1]),years=m2_n,cagr=float(mg),multiple=float(m.iloc[-1]/m.iloc[0])),historical=rows,common=common,forward=forward,model=dict(a=float(a),b=float(b),n=len(btc),start=str(btc.index.min().date()),end=str(btc.index.max().date()),origin='2009-01-03',days_now=(NOW-ORIGIN).days,median_now=trend(NOW),actual=P,actual_vs_trend=P/trend(NOW)-1,r_squared=float(1-np.sum(resid**2)/np.sum((y-y.mean())**2)),residual_sd_log10=float(np.std(resid,ddof=2)),residual_lag1=float(np.corrcoef(resid[:-1],resid[1:])[0,1]),sensitivity=sensitivity,holdout_2019_fit_current=trend(NOW,bb,aa),max_daily_drawdown=float((btc/btc.cummax()-1).min())),sp_ytd=float(sp_ytd-1),example=dict(nominal=1.07**30,m2=1.067**30,adjusted=(1.07/1.067)**30,adjusted_cagr=1.07/1.067-1))
(ROOT/'calculations.json').write_text(json.dumps(summary,indent=2))
print(json.dumps(summary,indent=2))