#!/usr/bin/env python3
"""COMPLEXA C320 — reproduce retrospective macrodata / realcons forecasts.

Uses only Python standard library. No look-ahead: model fitting receives strictly
past observations. Annual series uses four quarterly observations per full year.
Command: python reproduce.py --data macrodata.csv --out output
"""
from __future__ import annotations
import argparse,csv,hashlib,json,math
from collections import defaultdict
from pathlib import Path

def read_data(file):
    rows=list(csv.DictReader(file.open(newline='',encoding='utf-8-sig')))
    years=[int(float(x['year'])) for x in rows]
    quarters=[int(float(x['quarter'])) for x in rows]
    values=[float(x['realcons']) for x in rows]
    assert len(rows)==203 and (years[0],quarters[0])==(1959,1) and (years[-1],quarters[-1])==(2009,3)
    assert all((years[i]*4+quarters[i])==(years[0]*4+quarters[0])+i for i in range(len(rows)))
    assert all(v>0 for v in values)
    return list(zip(years,quarters,values))

def write_csv(file,rows,columns):
    with file.open('w',encoding='utf8',newline='') as f:
        wr=csv.DictWriter(f,fieldnames=columns,lineterminator='\n');wr.writeheader();wr.writerows(rows)

def growths(xs):return [100*math.log(xs[i]/xs[i-1]) for i in range(1,len(xs))]

def fit_ar(g,pairs=None):
    # g contains growths through the last observed level, not future values.
    x,y=g[:-1],g[1:]
    if pairs is not None:x,y=x[-pairs:],y[-pairs:]
    if len(x)<2:raise ValueError('AR fitting requires at least two paired growths')
    xm=sum(x)/len(x);ym=sum(y)/len(y)
    ss=sum((v-xm)**2 for v in x)
    phi=sum((u-xm)*(v-ym) for u,v in zip(x,y))/ss if ss else 0.0
    intercept=ym-phi*xm
    return intercept,phi,len(x)

def predict(values,origin,horizon,rule,window=None):
    # Caller passes actual observations but this function slices them at origin.
    past=values[:origin+1]
    g=growths(past)
    if rule=='persist':rate=0.;a=phi=n_pairs=None
    elif rule=='mean':rate=sum(g)/len(g);a=phi=n_pairs=None
    elif rule=='ar1':a,phi,n_pairs=fit_ar(g,window);rate=None
    else:raise ValueError(rule)
    next_rate=g[-1]
    levels=[];growth_prediction=[];level=past[-1]
    for _ in range(horizon):
        next_rate=(a+phi*next_rate) if rule=='ar1' else rate
        growth_prediction.append(next_rate)
        level*=math.exp(next_rate/100)
        levels.append(level)
    return levels,growth_prediction,a,phi,n_pairs

def metrics(pairs):
    errs=[a-b for a,b in pairs]
    return {'n':len(errs),'mae':sum(abs(e) for e in errs)/len(errs),
            'rmse':math.sqrt(sum(e*e for e in errs)/len(errs)),
            'bias':sum(errs)/len(errs)}

def main(data,out):
    out.mkdir(parents=True,exist_ok=True)
    data=data.resolve();raw=data.read_bytes();sha=hashlib.sha256(raw).hexdigest()
    points=read_data(data)
    xs=[v for _,_,v in points]
    periods=[f'{y}-T{q}' for y,q,v in points]
    assert xs[123]==5261.7 and xs[124]==5303.3 and xs[-1]==9256.0
    dev=[]
    for target in range(84,124): # 1980Q1–1989Q4
        for window in (20,40,80,None):
            lv,gg,a,phi,n=predict(xs,target-1,1,'ar1',window)
            actual_growth=100*math.log(xs[target]/xs[target-1]);err=actual_growth-gg[0]
            dev.append(dict(origin=periods[target-1],target=periods[target],window=window if window is not None else 'expanding',
                            predicted_growth=round(gg[0],12),actual_growth=round(actual_growth,12),error=round(err,12),pairs=n))
    devscores={k:sum(abs(float(r['error'])) for r in dev if r['window']==k)/40 for k in (20,40,80,'expanding')}
    print('Development',devscores)
    assert all(abs(devscores[k]-v)<0.00015 for k,v in {20:.6773,40:.6205,80:.6098,'expanding':.6027}.items())
    one=[];four=[]
    for horizon in (1,4):
        for origin in range(123,203-horizon):
            for rule in ('persist','mean','ar1'):
                lv,gg,a,phi,n=predict(xs,origin,horizon,rule)
                pred=lv[-1];truth=xs[origin+horizon]
                row=dict(origin=periods[origin],target=periods[origin+horizon],rule=rule,horizon=horizon,
                         observed=truth,prediction=round(pred,12),error=round(truth-pred,12),
                         predicted_growth=round(sum(gg),12),actual_growth=round(100*math.log(truth/xs[origin]),12),
                         intercept='' if a is None else round(a,12),phi='' if phi is None else round(phi,12),
                         fitted_pairs='' if n is None else n)
                (one if horizon==1 else four).append(row)
    assert len(one)==79*3 and len(four)==76*3
    for k,v in {'persist':56.86,'mean':30.96,'ar1':28.26}.items():
        z=metrics([(float(r['observed']),float(r['prediction'])) for r in one if r['rule']==k]);print('One',k,z)
        assert abs(z['mae']-v)<0.025,(k,z,v)
    for k,v in {'persist':220.54,'mean':95.93,'ar1':91.11}.items():
        z=metrics([(float(r['observed']),float(r['prediction'])) for r in four if r['rule']==k]);print('Four',k,z)
        assert abs(z['mae']-v)<0.025,(k,z,v)
    # First four-step path from 1989Q4, matches published four quarterly forecasts.
    first_levels,first_g,a,phi,n=predict(xs,123,4,'ar1')
    assert all(abs(x-y)<0.025 for x,y in zip(first_levels,[5304.62,5351.83,5400.35,5449.52]))
    annual={}
    for y in range(1959,2009):
        vs=[v for yy,q,v in points if yy==y]
        assert len(vs)==4
        annual[y]=sum(vs)/4
    annual_predictions=[]
    for y in range(1990,2009):
        origin=next(i for i,(yy,q,v) in enumerate(points) if yy==y-1 and q==4)
        observed=annual[y]
        annual_past=[annual[yy] for yy in range(1959,y)]
        annual_g=growths(annual_past)
        for rule in ('persist','mean'):
            r=(annual_past[-1] if rule=='persist' else annual_past[-1]*math.exp(sum(annual_g)/len(annual_g)/100))
            annual_predictions.append(dict(year=y,source='annual',rule=rule,observed=observed,prediction=r,error=observed-r))
        for rule in ('persist','mean','ar1'):
            lv,g,*_=predict(xs,origin,4,rule)
            r=sum(lv)/4
            annual_predictions.append(dict(year=y,source='quarterly',rule=rule,observed=observed,prediction=r,error=observed-r))
    expected={('annual','persist'):217.22,('annual','mean'):76.66,
              ('quarterly','persist'):137.76,('quarterly','mean'):57.13,('quarterly','ar1'):53.93}
    for k,v in expected.items():
        sel=[r for r in annual_predictions if (r['source'],r['rule'])==k]
        res=metrics([(r['observed'],r['prediction']) for r in sel]);print('Annual',k,res)
        assert abs(res['mae']-v)<0.03,(k,res,v)
    col=['origin','target','rule','horizon','observed','prediction','error','predicted_growth','actual_growth','intercept','phi','fitted_pairs']
    write_csv(out/'development.csv',dev,['origin','target','window','predicted_growth','actual_growth','error','pairs'])
    write_csv(out/'rolling_one_quarter.csv',one,col)
    write_csv(out/'rolling_four_quarters.csv',four,col)
    write_csv(out/'rolling_annual_average.csv',annual_predictions,['year','source','rule','observed','prediction','error'])
    write_csv(out/'realcons_quarterly.csv',[{'period':period,'realcons':v} for period,v in zip(periods,xs)],['period','realcons'])
    summary={'source_sha256':sha,'observations':len(xs),'range':[periods[0],periods[-1]],
            'development_mae':{str(k):v for k,v in devscores.items()},'one_quarter':{k:metrics([(r['observed'],r['prediction']) for r in one if r['rule']==k]) for k in ('persist','mean','ar1')},
            'four_quarters':{k:metrics([(r['observed'],r['prediction']) for r in four if r['rule']==k]) for k in ('persist','mean','ar1')},
            'annual':{f'{s}_{k}':metrics([(r['observed'],r['prediction']) for r in annual_predictions if (r['source'],r['rule'])==(s,k)]) for s,k in expected},
            'first_ar4_levels':first_levels,'last_q2009':[(r['rule'],r['prediction']) for r in one if r['target']=='2009-T3']}
    (out/'validation.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2,sort_keys=True)+'\n',encoding='utf8')
    print('PASS C320: archived 203-quarter realcons and all published summary tables')
if __name__=='__main__':
    p=argparse.ArgumentParser();p.add_argument('--data',type=Path,default=Path('macrodata.csv'));p.add_argument('--out',type=Path,default=Path('output'));a=p.parse_args();main(a.data,a.out)