"""sharpe-sortino-max-drawdown-by-hand.py — Sharpe, Sortino and maximum drawdown computed step by step on FRED's S&P 500 series.

Input : datasets/sharpe-sortino-max-drawdown-by-hand.csv - the derived daily simple return series
        computed from FRED's SP500 daily close (a price index without dividends, on FRED's ten-year
        rolling window) as it stood at the 2026-09-05 retrieval, alongside DGS3MO, the daily 3-month
        Treasury yield used as the risk-free proxy. The index levels themselves are not shipped. S&P Dow
        Jones Indices LLC states in the FRED series notes that "Reproduction of S&P 500 in any form is
        prohibited except with the prior written permission of S&P Dow Jones Indices LLC ("S&P")", so
        this file carries only the returns computed from the series, which are this site's own output.
        Source: S&P Dow Jones Indices LLC, https://fred.stlouisfed.org/series/SP500; DGS3MO from the
        Board of Governors of the Federal Reserve System via FRED.
        Run with --download to pull both series from FRED and rebuild the file end to end. FRED's
        ten-year window rolls forward, so a pull today covers a later sample than the article reports.
Method: r_t = P_t/P_{t-1} - 1 on trading days; rf_t = DGS3MO_t/100/252 (previous close, forward-filled);
        excess e_t = r_t - rf_t. Sharpe = mean(e)/std(e, ddof=1) * sqrt(252). Sortino uses the downside
        deviation sqrt(mean(min(e,0)^2)) * sqrt(252). Max drawdown = min over t of P_t / max_{s<=t} P_s - 1.
        Nothing here is a forecast; the numbers describe one sample of one price index.
Run   : python code/sharpe-sortino-max-drawdown-by-hand.py [--download]
Needs : Python 3.13, pandas 3.0.2, numpy 2.3.
"""
import io
import math
import sys
import urllib.request

import numpy as np
import pandas as pd

CSV = "datasets/sharpe-sortino-max-drawdown-by-hand.csv"
DAYS = 252


def download():
    """Pull both series from FRED and write the derived returns plus DGS3MO; no index levels are stored."""
    frames = []
    for sid in ("SP500", "DGS3MO"):
        req = urllib.request.Request(f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={sid}",
                                     headers={"User-Agent": "prism-data-lab/1.0", "Accept": "*/*"})
        txt = urllib.request.urlopen(req, timeout=60).read().decode()
        frames.append(pd.read_csv(io.StringIO(txt), na_values=".", index_col=0).iloc[:, 0].rename(sid))
    df = pd.concat(frames, axis=1)
    df.index.name = "date"
    df["return_simple"] = df["SP500"].dropna().pct_change()
    df[["return_simple", "DGS3MO"]].to_csv(CSV)


def rebased_level(returns: pd.Series) -> pd.Series:
    """Rebuild a rebased performance path, 1.0 on the first day of the sample, from the shipped returns.

    The base date is the row immediately before the first return, which is the first observation in the
    sample. Every statistic printed below - total and annualised return, volatility, Sharpe, Sortino and
    drawdown - is invariant to the scale of the path, so the rebased path reproduces every published
    figure without storing a single index level.
    """
    r = returns.dropna()
    base_date = returns.index[returns.index.get_loc(r.index[0]) - 1]
    level = (1.0 + r).cumprod()
    return pd.concat([pd.Series([1.0], index=[base_date]), level])


def stats(px: pd.Series, rf_annual: pd.Series) -> dict:
    r = px.pct_change().dropna()
    rf = (rf_annual.shift(1).ffill().reindex(r.index) / 100.0) / DAYS
    e = (r - rf).dropna()
    mean_d, sd_d = e.mean(), e.std(ddof=1)
    downside = math.sqrt((np.minimum(e, 0.0) ** 2).mean())
    peak = px.cummax()
    dd = px / peak - 1.0
    trough = dd.idxmin()
    peak_date = px[:trough].idxmax()
    rec = px[trough:]
    recovered = rec[rec >= px[peak_date]].index
    return {
        "days": len(e), "start": px.index[0], "end": px.index[-1],
        "total_return": px.iloc[-1] / px.iloc[0] - 1.0,
        "ann_return": (px.iloc[-1] / px.iloc[0]) ** (DAYS / len(r)) - 1.0,
        "ann_vol": sd_d * math.sqrt(DAYS),
        "mean_rf": rf.mean() * DAYS,
        "sharpe": mean_d / sd_d * math.sqrt(DAYS),
        "sortino": mean_d / downside * math.sqrt(DAYS),
        "max_dd": dd.min(), "dd_peak": peak_date, "dd_trough": trough,
        "dd_recovered": recovered[0] if len(recovered) else None,
        "worst_day": r.min(), "worst_day_date": r.idxmin(), "best_day": r.max(), "best_day_date": r.idxmax(),
    }


def main():
    if "--download" in sys.argv:
        download()
    df = pd.read_csv(CSV, parse_dates=["date"], index_col="date")
    px = rebased_level(df["return_simple"])
    rf = df["DGS3MO"]
    s = stats(px, rf)
    print(f"sample {s['start']:%Y-%m-%d} to {s['end']:%Y-%m-%d}: {s['days']} daily excess returns")
    print("| statistic | value |")
    print("|---|---|")
    print(f"| total price return | {100 * s['total_return']:.1f}% |")
    print(f"| annualised price return (geometric) | {100 * s['ann_return']:.2f}% |")
    print(f"| annualised volatility (ddof=1) | {100 * s['ann_vol']:.2f}% |")
    print(f"| mean risk-free rate used (DGS3MO) | {100 * s['mean_rf']:.2f}% |")
    print(f"| Sharpe ratio (annualised, daily excess returns) | {s['sharpe']:.2f} |")
    print(f"| Sortino ratio (annualised) | {s['sortino']:.2f} |")
    print(f"| maximum drawdown | {100 * s['max_dd']:.1f}% |")
    print(f"| drawdown peak / trough / recovery | {s['dd_peak']:%Y-%m-%d} / {s['dd_trough']:%Y-%m-%d} / "
          f"{s['dd_recovered']:%Y-%m-%d} |")
    print(f"| worst day | {100 * s['worst_day']:.2f}% on {s['worst_day_date']:%Y-%m-%d} |")
    print(f"| best day | {100 * s['best_day']:.2f}% on {s['best_day_date']:%Y-%m-%d} |")
    print("\n| calendar year | trading days | price return | ann. vol | Sharpe | Sortino | max drawdown within year |")
    print("|---|---|---|---|---|---|---|")
    for y, p in px.groupby(px.index.year):
        if len(p) < 40:
            continue
        # include the prior close so the first day of the year has a return
        prev = px[px.index < p.index[0]]
        pp = pd.concat([prev.tail(1), p]) if len(prev) else p
        t = stats(pp, rf)
        print(f"| {y} | {t['days']} | {100 * t['total_return']:.1f}% | {100 * t['ann_vol']:.1f}% | {t['sharpe']:.2f} | "
              f"{t['sortino']:.2f} | {100 * t['max_dd']:.1f}% |")
    # the same Sharpe from monthly returns, to show the frequency dependence
    m = px.resample("ME").last().pct_change().dropna()
    rf_m = (rf.resample("ME").mean().reindex(m.index) / 100.0) / 12
    em = m - rf_m
    print(f"\nSharpe from monthly returns: {em.mean() / em.std(ddof=1) * math.sqrt(12):.2f} ({len(em)} months); "
          f"skew of daily returns {px.pct_change().dropna().skew():.2f}, excess kurtosis {px.pct_change().dropna().kurt():.1f}")


if __name__ == "__main__":
    main()