#!/usr/bin/env python3 """Reproduce a monthly US market risk/return forecasting illustration. Downloads the Kenneth French daily Fama-French factor file, aggregates daily market excess returns and squared daily returns to calendar months, runs two simple expanding-window benchmark comparisons (OOS 2000-2025), and writes CSV, JSON and an SVG chart. Python standard library only. """ from __future__ import annotations import csv, io, json, math, urllib.request, zipfile from collections import defaultdict from datetime import date from pathlib import Path ROOT = next(p for p in Path(__file__).resolve().parents if (p / "src" / "content").is_dir()) OUT = ROOT / "public" / "data" / "research" IMG = ROOT / "public" / "images" / "research" URL = "https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_daily_CSV.zip" START, END, TRAIN_MONTHS = "2000-01", "2025-12", 60 def download(): req = urllib.request.Request(URL, headers={"User-Agent": "Lubasiverse-research-reproducibility/1.0"}) with urllib.request.urlopen(req, timeout=60) as r: payload = r.read() z = zipfile.ZipFile(io.BytesIO(payload)) name = next(n for n in z.namelist() if n.lower().endswith(".csv")) return z.read(name).decode("utf-8-sig", errors="replace") def rows_from_source(text): daily = [] for line in text.splitlines(): line = line.strip() if not line or not line[:8].isdigit(): continue cells = [x.strip() for x in line.split(",")] if len(cells) < 5: continue try: dt = date(int(cells[0][:4]), int(cells[0][4:6]), int(cells[0][6:8])) mktrf, rf = float(cells[1]), float(cells[4]) except (ValueError, IndexError): continue if mktrf < -90 or rf < -90: continue daily.append((dt, mktrf / 100.0, rf / 100.0)) return daily def monthly_aggregate(daily): buckets = defaultdict(list) for dt, mktrf, rf in daily: buckets[dt.strftime("%Y-%m")].append((mktrf, rf)) out = [] for month, vals in sorted(buckets.items()): # Compound daily simple returns separately; RF is the matched risk-free series. mkt_gross = math.prod(1 + x for x, _ in vals) rf_gross = math.prod(1 + x for _, x in vals) excess = mkt_gross - rf_gross # Realised variance proxy: sum of squared daily market excess returns. rv = sum(x * x for x, _ in vals) out.append({"month": month, "n_days": len(vals), "excess_return": excess, "realised_variance": rv, "realised_vol_annualised": math.sqrt(12 * rv)}) return out def mean(xs): return sum(xs) / len(xs) def qlike(y, f): # Patton's QLIKE plus outcome-only terms: nonnegative and scale-free. y = max(y, 1e-12); f = max(f, 1e-12) ratio = y / f return ratio - math.log(ratio) - 1 def hac_tstat(diffs, max_lag=12): """Newey-West mean-loss-difference t statistic (Bartlett weights).""" n = len(diffs); mu = mean(diffs) centered = [x - mu for x in diffs] long_run_var = sum(x*x for x in centered) / n for lag in range(1, min(max_lag, n-1) + 1): cov = sum(centered[t] * centered[t-lag] for t in range(lag, n)) / n long_run_var += 2 * (1 - lag / (max_lag + 1)) * cov se = math.sqrt(max(0, long_run_var) / n) tstat = mu / se if se else float('nan') pvalue = math.erfc(abs(tstat) / math.sqrt(2)) if math.isfinite(tstat) else None return {"mean_loss_difference_alt_minus_benchmark": mu, "newey_west_lags_months": min(max_lag, n-1), "newey_west_t_statistic": tstat, "two_sided_normal_approx_p_value": pvalue} def experiment(monthly): ix = {r["month"]: i for i, r in enumerate(monthly)} first, last = ix[START], ix[END] if first < TRAIN_MONTHS or last < first: raise ValueError(f"Insufficient data for requested window: {START}..{END}") forecasts = [] for i in range(first, last + 1): hist = monthly[:i] if len(hist) < TRAIN_MONTHS: continue yret = monthly[i]["excess_return"] yrv = monthly[i]["realised_variance"] rbase = mean([x["excess_return"] for x in hist]) ralt = mean([x["excess_return"] for x in hist[-12:]]) vbase = mean([x["realised_variance"] for x in hist]) valt = hist[-1]["realised_variance"] forecasts.append({"month": monthly[i]["month"], "actual_excess_return": yret, "return_historical_mean_forecast": rbase, "return_trailing_12m_forecast": ralt, "actual_realised_variance_proxy": yrv, "variance_historical_mean_forecast": vbase, "variance_last_month_forecast": valt, "return_sqerr_benchmark": (yret-rbase)**2, "return_sqerr_alt": (yret-ralt)**2, "variance_qlike_benchmark": qlike(yrv,vbase), "variance_qlike_alt": qlike(yrv,valt)}) ret_diffs = [x["return_sqerr_alt"] - x["return_sqerr_benchmark"] for x in forecasts] vol_diffs = [x["variance_qlike_alt"] - x["variance_qlike_benchmark"] for x in forecasts] rb = mean([x["return_sqerr_benchmark"] for x in forecasts]) ra = mean([x["return_sqerr_alt"] for x in forecasts]) vb = mean([x["variance_qlike_benchmark"] for x in forecasts]) va = mean([x["variance_qlike_alt"] for x in forecasts]) periods=[] for start,end in [(2000,2004),(2005,2009),(2010,2014),(2015,2019),(2020,2025)]: block=[x for x in forecasts if start <= int(x["month"][:4]) <= end] if not block: continue rbp=mean([x["return_sqerr_benchmark"] for x in block]); rap=mean([x["return_sqerr_alt"] for x in block]) vbp=mean([x["variance_qlike_benchmark"] for x in block]); vap=mean([x["variance_qlike_alt"] for x in block]) periods.append({"period":f"{start}-{end}","n_months":len(block),"return_oos_r2":1-rap/rbp, "variance_qlike_improvement_percent":(vbp-vap)/vbp*100}) result = {"source_url": URL, "source_description": "Kenneth R. French Data Library, daily Fama/French 3 factors; daily Mkt-RF and RF, downloaded at run time.", "return_units_in_csv": "decimal monthly excess return (0.01 = 1%)", "risk_proxy_units_in_csv": "squared decimal daily-return units, summed over each month", "data_start_month": monthly[0]["month"], "initial_estimation_months": first, "subperiod_sensitivity": periods, "sample_start": forecasts[0]["month"], "sample_end": forecasts[-1]["month"], "n_months": len(forecasts), "training_window_months_minimum": TRAIN_MONTHS, "return_benchmark": "Expanding historical mean of monthly market excess returns", "return_alternative": "Mean of the preceding 12 monthly market excess returns", "risk_proxy": "Within-month sum of squared daily market excess returns", "risk_benchmark": "Expanding historical mean of monthly realised-variance proxies", "risk_alternative": "Previous month's realised-variance proxy (persistence)", "return_loss": "Monthly squared forecast error", "return_oos_r2": 1 - ra/rb, "return_rmse_benchmark_percentage_points": math.sqrt(rb)*100, "return_rmse_alternative_percentage_points": math.sqrt(ra)*100, "return_loss_comparison": hac_tstat(ret_diffs), "variance_loss": "Normalised QLIKE: y/f - log(y/f) - 1, equivalent for rankings to Patton's QLIKE up to an outcome-only term", "variance_mean_qlike_benchmark": vb, "variance_mean_qlike_alternative": va, "variance_qlike_improvement_percent": (vb-va)/vb*100, "variance_loss_comparison": hac_tstat(vol_diffs), "caveat": "Descriptive single-market comparison; relative skill is not a commensurable measure of economic value. HAC tests are approximate; the daily squared-return proxy is noisy and is not the latent variance itself."} return forecasts, result def svg_chart(rows, result, path): W,H=1120,740; left,right=82,1070; top1,bot1=100,350; top2,bot2=430,650 ret=[x["actual_excess_return"]*100 for x in rows] rv=[math.sqrt(12*x["actual_realised_variance_proxy"])*100 for x in rows] def coords(series, top, bot): lo,hi=min(series),max(series) pad=(hi-lo)*.08 or 1 lo-=pad; hi+=pad pts=[] for i,v in enumerate(series): x=left+(right-left)*i/max(1,len(series)-1); y=bot-(v-lo)/(hi-lo)*(bot-top) pts.append(f"{x:.1f},{y:.1f}") return " ".join(pts),lo,hi p1,lo1,hi1=coords(ret,top1,bot1); p2,lo2,hi2=coords(rv,top2,bot2) # Rolling 12-month descriptive averages, not forecasts. ra=[]; va=[] for i in range(len(rows)): st=max(0,i-11); ra.append(mean(ret[st:i+1])); va.append(mean(rv[st:i+1])) pr,_a,_b=coords(ra,top1,bot1); pv,_a,_b=coords(va,top2,bot2) s=f''' US market monthly excess returns and realised volatility proxy, 2000–2025 Two time-series panels. Top: monthly market excess return in percent with a 12-month moving average. Bottom: annualised square root of the monthly realised-variance proxy in percent, with a 12-month moving average. Daily market factor data from Kenneth French's Data Library. Risk and returns are not the same forecasting target US market excess returns · monthly · 2000–2025 · {len(rows)} out-of-sample months Monthly excess return (%) Annualised realised-volatility proxy (%) monthly observation12-month mean Teal: monthly observations · white: trailing 12-month descriptive average. The lower panel is a daily-squared-return proxy, not latent true risk. 2000201020202025 2000201020202025''' path.write_text(s, encoding="utf-8") def svg_skill_chart(rows, path): """Cumulative benchmark loss minus alternative loss; positive favours the latter.""" W,H=1120,650; left,right=92,1070 panels=[("Returns: expanding-mean SSE minus trailing-12-month SSE", "return_sqerr_benchmark", "return_sqerr_alt", 115,280), ("Risk proxy: expanding-mean QLIKE minus persistence QLIKE", "variance_qlike_benchmark", "variance_qlike_alt", 390,555)] out=[f'Cumulative out-of-sample loss advantageTwo panels chart cumulative benchmark loss minus alternative forecast loss. Positive values favour the alternative. Separate vertical scales are used.Which simple forecast did better?Cumulative benchmark loss − alternative loss · positive means lower cumulative loss for the alternative'] for title,base,alt,top,bot in panels: cumulative=[]; total=0 for r in rows: total += r[base] - r[alt]; cumulative.append(total) lo=min(0.0,min(cumulative)); hi=max(0.0,max(cumulative)); pad=(hi-lo)*.08 or 1; lo-=pad; hi+=pad pts=[] for i,v in enumerate(cumulative): x=left+(right-left)*i/max(1,len(cumulative)-1); y=bot-(v-lo)/(hi-lo)*(bot-top) pts.append(f"{x:.1f},{y:.1f}") zero=bot-(0-lo)/(hi-lo)*(bot-top) out.append(f'{title}') out.append(f'{rows[0]["month"]}{rows[-1]["month"]}End: {total:.4f}') out.append('Different loss functions and units: compare each forecast only to its own benchmark, not across panels.') path.write_text("".join(out),encoding="utf-8") def main(): OUT.mkdir(parents=True, exist_ok=True); IMG.mkdir(parents=True, exist_ok=True) daily=rows_from_source(download()) monthly=monthly_aggregate(daily) forecasts,summary=experiment(monthly) with (OUT/"risk-returns-monthly-forecasts.csv").open("w",newline="",encoding="utf-8") as f: w=csv.DictWriter(f,fieldnames=list(forecasts[0]),lineterminator="\n"); w.writeheader(); w.writerows(forecasts) (OUT/"risk-returns-forecast-summary.json").write_text(json.dumps(summary,indent=2)+"\n",encoding="utf-8") svg_chart(forecasts,summary,IMG/"risk-returns-forecast-results.svg") svg_skill_chart(forecasts,IMG/"risk-returns-cumulative-loss.svg") print(json.dumps(summary,indent=2)) if __name__ == "__main__": main()