#!/usr/bin/env python3
"""
Monthly risk-free return series from FRED TB3MS (3-month Treasury bill, secondary market).

CONVENTION — DOCUMENTED, NOT ASSUMED
TB3MS is published on a BANK DISCOUNT basis (percent per annum, 360-day year).
A discount rate d is NOT an investment yield: it is quoted against face value, not
against the price paid. Dividing d by 12 would overstate nothing but would also
mis-state the actual return earned, so we convert explicitly:

  1. Discount rate -> price per 100 face for a 91-day bill:
        P = 100 * (1 - d * 91/360)
  2. Price -> bond-equivalent (investment) yield, actual/365:
        BEY = (100 - P)/P * 365/91
  3. BEY -> monthly return, compounded:
        r_m = (1 + BEY)^(1/12) - 1

Step 3 uses geometric conversion rather than BEY/12 so that twelve monthly returns
compound to the annual yield exactly.

Validation reference: Damodaran histretSP.xls "Returns by year" column D,
"3-month T.Bill".
"""
import csv, os, math, json
from collections import OrderedDict
HERE=os.path.dirname(os.path.abspath(__file__))

def load(fn):
    out=OrderedDict()
    for r in list(csv.reader(open(os.path.join(HERE,'data',fn))))[1:]:
        if len(r)<2 or r[1] in ('.',''): continue
        y,m,_=r[0].split('-'); out['%04d-%02d'%(int(y),int(m))]=float(r[1])
    return out

def discount_to_monthly(d_pct, days=91):
    d=d_pct/100.0
    price=100.0*(1.0 - d*days/360.0)
    bey=(100.0-price)/price*365.0/days
    return (1.0+bey)**(1.0/12.0)-1.0, bey

def build():
    tb=load('TB3MS.csv')
    return OrderedDict((k, discount_to_monthly(v)[0]) for k,v in tb.items()), tb

def annual(monthly):
    yr=OrderedDict()
    for k,v in monthly.items(): yr.setdefault(int(k[:4]),[]).append(v)
    out=OrderedDict()
    for y,v in yr.items():
        if len(v)==12:
            p=1.0
            for r in v: p*=(1+r)
            out[y]=p-1
    return out

if __name__=='__main__':
    # self-test: 12 monthly returns must compound back to the BEY
    for d in (1.0, 5.0, 10.0, 15.0):
        m,bey=discount_to_monthly(d)
        assert abs((1+m)**12-1-bey)<1e-12, 'compounding identity failed'
    print('self-test PASSED: 12 monthly returns compound exactly to BEY')

    mr,raw=build(); ann=annual(mr)
    print('monthly obs %d (%s..%s) | complete years %d (%d..%d)'%(len(mr),min(mr),max(mr),len(ann),min(ann),max(ann)))

    import xlrd, warnings; warnings.filterwarnings('ignore')
    ws=xlrd.open_workbook(os.path.join(HERE,'data','histretSP.xls')).sheet_by_name('Returns by year')
    dam={}
    for r in range(20,ws.nrows):
        y=ws.cell_value(r,0); v=ws.cell_value(r,3)
        if isinstance(y,float) and isinstance(v,float) and 1900<y<2100: dam[int(y)]=float(v)
    com=sorted(set(ann)&set(dam))
    a=[ann[y] for y in com]; b=[dam[y] for y in com]; n=len(com)
    ma,mb=sum(a)/n,sum(b)/n
    sa=math.sqrt(sum((x-ma)**2 for x in a)/(n-1)); sb=math.sqrt(sum((x-mb)**2 for x in b)/(n-1))
    cov=sum((a[i]-ma)*(b[i]-mb) for i in range(n))/(n-1)
    d=[a[i]-b[i] for i in range(n)]; ab=sorted(abs(x) for x in d)
    def cagr(v):
        p=1.0
        for x in v: p*=(1+x)
        return p**(1/len(v))-1
    print()
    print('='*74); print('  T-BILL VALIDATION vs Damodaran 3-month T.Bill'); print('='*74)
    print('common years %d (%d-%d)'%(n,min(com),max(com)))
    print('correlation                : %.4f'%(cov/(sa*sb)))
    print('mean absolute difference   : %.3f pp'%(sum(ab)/n*100))
    print('median absolute difference : %.3f pp'%((ab[n//2] if n%2 else (ab[n//2-1]+ab[n//2])/2)*100))
    print('bias                       : %+.3f pp'%(sum(d)/n*100))
    print('annualised ours/Damodaran  : %.3f%% / %.3f%%  (diff %+.3f pp)'%(cagr(a)*100,cagr(b)*100,(cagr(a)-cagr(b))*100))
    print('\nlargest discrepancies:')
    for y in sorted(com,key=lambda y:-abs(ann[y]-dam[y]))[:6]:
        print('   %d  ours %6.3f%%  Damodaran %6.3f%%  diff %+.3f pp'%(y,ann[y]*100,dam[y]*100,(ann[y]-dam[y])*100))
    ok = cov/(sa*sb)>0.99 and abs(sum(d)/n)<0.005
    print('\nVERDICT: %s (corr>0.99 and |bias|<0.5pp)'%('PASS' if ok else 'FAIL'))
    with open(os.path.join(HERE,'normalised','tbill_monthly_returns.csv'),'w',newline='') as f:
        w=csv.writer(f); w.writerow(['month','tbill3m_monthly_return'])
        for k,v in mr.items(): w.writerow([k,'%.10f'%v])
    json.dump({'corr':cov/(sa*sb),'mae_pp':sum(ab)/n*100,'bias_pp':sum(d)/n*100,
               'cagr_ours':cagr(a),'cagr_dam':cagr(b),'n':n,'pass':bool(ok)},
              open(os.path.join(HERE,'normalised','tbill_validation.json'),'w'),indent=2)