#!/usr/bin/env python3
"""
Constant-maturity 10-year Treasury monthly total return, built from FRED GS10.

METHOD
  Each month we hold a par 10-year bond whose coupon rate equals the prior
  month-end yield y0 (price = 100 at purchase). One month later the bond has
  9y11m remaining and is repriced at the new yield y1 with semi-annual
  discounting. We then roll into a fresh par 10-year bond (constant maturity).

    dirty_price(y1) = SUM_t c / (1+y1/2)^t  +  100 / (1+y1/2)^n
      c = 100 * y0 / 2            semi-annual coupon
      n = 19.8333 periods         9y11m expressed in semi-annual periods
      t = coupon times in periods, first at n - floor(n) offset

  total_return_month = dirty_price(y1) / 100 - 1

  The discounted cash flows ALREADY contain the coupon the holder is entitled
  to, so no separate accrued-interest term is added.

CORRECTION 2026-09-02: an earlier version added an explicit accrued-coupon term
on top of the dirty price. That double-counted the coupon and inflated returns by
roughly y/12 per month (about +6.3pp a year). Detected by the unchanged-yield
invariant below and by validation against Damodaran. INVARIANT NOW ASSERTED IN
CODE: if y1 == y0 the monthly total return must equal y0/12 to within 5bp.

Conventions: month-end constant-maturity yields, semi-annual compounding, full
repricing (no duration approximation), no transaction costs, no within-month
coupon reinvestment.
"""
import csv, os
from collections import OrderedDict

HERE = os.path.dirname(os.path.abspath(__file__))

def dirty_price(y_annual, years_remaining, coupon_annual_rate, face=100.0):
    y = y_annual / 2.0
    c = face * coupon_annual_rate / 2.0
    n = years_remaining * 2.0
    if y <= -0.99:
        raise ValueError('implausible yield')
    n_full = int(n // 1)
    frac = n - n_full                       # periods until the next coupon
    times = [frac + i for i in range(n_full + 1)] if frac > 1e-12 else \
            [i + 1 for i in range(n_full)]
    times = [t for t in times if t <= n + 1e-9]
    pv = sum(c / ((1 + y) ** t) for t in times)
    pv += face / ((1 + y) ** n)
    return pv

def monthly_total_returns(gs10):
    keys = sorted(gs10)
    out = OrderedDict()
    for i in range(1, len(keys)):
        y0, y1 = gs10[keys[i-1]], gs10[keys[i]]
        out[keys[i]] = dirty_price(y1, 10.0 - 1.0/12.0, y0) / 100.0 - 1.0
    return out

def _self_test():
    """Unchanged yield => monthly total return must be y/12."""
    worst = 0.0
    for y in (0.01, 0.02, 0.05, 0.08, 0.12, 0.15):
        r = dirty_price(y, 10.0 - 1.0/12.0, y) / 100.0 - 1.0
        worst = max(worst, abs(r - y/12.0))
    assert worst < 5e-4, 'INVARIANT FAILED: unchanged-yield return != y/12 (max dev %.6f)' % worst
    return worst

def load_gs10(path):
    out = OrderedDict()
    for row in list(csv.reader(open(path)))[1:]:
        if len(row) < 2 or row[1] in ('.', ''):
            continue
        y, m, _ = row[0].split('-')
        out['%04d-%02d' % (int(y), int(m))] = float(row[1]) / 100.0
    return out

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

if __name__ == '__main__':
    dev = _self_test()
    print('self-test PASSED: unchanged-yield invariant, max deviation %.2e' % dev)
    gs10 = load_gs10(os.path.join(HERE, 'data', 'GS10.csv'))
    mr = monthly_total_returns(gs10)
    ann = to_calendar_year(mr)
    print('monthly obs %d (%s..%s) | complete years %d (%d..%d)'
          % (len(mr), min(mr), max(mr), len(ann), min(ann), max(ann)))
    os.makedirs(os.path.join(HERE, 'normalised'), exist_ok=True)
    with open(os.path.join(HERE,'normalised','bond_monthly_returns.csv'),'w',newline='') as f:
        w=csv.writer(f); w.writerow(['month','ust10_total_return_monthly'])
        for k,v in mr.items(): w.writerow([k,'%.10f'%v])
    with open(os.path.join(HERE,'normalised','bond_annual_returns.csv'),'w',newline='') as f:
        w=csv.writer(f); w.writerow(['year','ust10_total_return_annual'])
        for y,v in ann.items(): w.writerow([y,'%.10f'%v])
    print('wrote bond_monthly_returns.csv, bond_annual_returns.csv')