study.py

"""S&P 500: buy at the trend line (holding cash) vs bi-weekly automatic buying.Shiller monthly data, real total return, ex-ante regression, equal contributions."""import csv, json, math, sys PAYCHECK = 100.0             # real dollars per bi-weekly paycheck (June 2023 dollars)MONTHLY = PAYCHECK * 26 / 12 # 26 paychecks a year, spread over 12 monthsMIN_HISTORY = 240            # months of history needed before the trend signal is usedrows = list(csv.DictReader(open(sys.argv[1])))rows = [r for r in rows if float(r['Dividend']) > 0 and float(r['Consumer Price Index']) > 0]dates = [r['Date'][:7] for r in rows]P = [float(r['SP500']) for r in rows]D = [float(r['Dividend']) for r in rows]CPI = [float(r['Consumer Price Index']) for r in rows]Y = [float(r['Long Interest Rate']) / 100 for r in rows]n = len(rows) # nominal total return index, dividends (annual rate / 12) reinvested each monthTR = [1.0]for t in range(1, n):    TR.append(TR[-1] * (P[t] + D[t] / 12) / P[t - 1])# real total return level in end-of-sample dollarsR = [TR[t] / CPI[t] * CPI[-1] for t in range(n)]LR = [math.log(x) for x in R] # expanding window OLS of ln(real TR) on month index, using data through t onlydef trend_line():    out = [None] * n    sx = sy = sxx = sxy = 0.0    for t in range(n):        sx += t; sy += LR[t]; sxx += t * t; sxy += t * LR[t]        k = t + 1        if k >= MIN_HISTORY:            b = (k * sxy - sx * sy) / (k * sxx - sx * sx)            a = (sy - b * sx) / k            out[t] = a + b * t    return outTL = trend_line() def run(start, end, cash_mode):    """cash_mode: 'zero' = cash earns 0% nominal; 'yield' = cash earns the 10 year yield (upper bound)."""    units_d = units_t = cash_t = 0.0    contrib = 0.0    flows_d, flows_t = [], []    curve = []    invested_months = 0; months = 0; cash_share = 0.0    for t in range(start, end + 1):        # existing cash loses (or earns) value in real terms before this month's paycheck arrives        if t > start:            infl = CPI[t - 1] / CPI[t]            grow = (1 + (Y[t - 1] / 12 if cash_mode == 'yield' else 0.0)) * infl            cash_t *= grow        price = R[t]        contrib += MONTHLY        units_d += MONTHLY / price                      # bi-weekly: buy at whatever the price is        cash_t += MONTHLY                               # trend: hold the paycheck as cash...        if LR[t] <= TL[t]:                              # ...and buy only when price is at or below the line            units_t += cash_t / price; cash_t = 0.0        v_d = units_d * price        v_t = units_t * price + cash_t        flows_d.append(-MONTHLY); flows_t.append(-MONTHLY)        curve.append((dates[t], contrib, v_d, v_t, cash_t))        months += 1; invested_months += (cash_t == 0.0); cash_share += cash_t / v_t    flows_d[-1] += curve[-1][2]; flows_t[-1] += curve[-1][3]    return curve, flows_d, flows_t, invested_months / months, cash_share / months def irr(flows):    lo, hi = -0.5, 0.5   # monthly rate bounds    f = lambda r: sum(c / (1 + r) ** i for i, c in enumerate(flows))    for _ in range(200):        mid = (lo + hi) / 2        if f(mid) > 0: lo = mid        else: hi = mid    return (1 + mid) ** 12 - 1 def max_dd(values):    peak = values[0]; worst = 0.0    for v in values:        peak = max(peak, v); worst = min(worst, v / peak - 1)    return worst def summarise(start, end, cash_mode):    curve, fd, ft, inv, cs = run(start, end, cash_mode)    c, vd, vt = curve[-1][1], curve[-1][2], curve[-1][3]    return dict(contributed=c, dca=vd, trend=vt, dca_x=vd / c, trend_x=vt / c,                dca_irr=irr(fd), trend_irr=irr(ft),                dca_dd=max_dd([x[2] for x in curve]), trend_dd=max_dd([x[3] for x in curve]),                pct_invested=inv, avg_cash_share=cs, months=len(curve)), curve end = n - 1start100 = end - 100 * 12 + 1first_signal = MIN_HISTORY - 1out = {'meta': dict(first=dates[0], last=dates[-1], start=dates[start100], paycheck=PAYCHECK, monthly=MONTHLY,                    min_history_months=MIN_HISTORY, months=n)}for mode in ('zero', 'yield'):    s, curve = summarise(start100, end, mode)    out[mode] = s    if mode == 'zero':        out['curve'] = curve# rolling windows: share of start months in which trend beats bi-weekly, per horizon (final wealth only)def final_values(start, end, cash_mode):    units_d = units_t = cash_t = 0.0    for t in range(start, end + 1):        if t > start:            cash_t *= (1 + (Y[t - 1] / 12 if cash_mode == 'yield' else 0.0)) * CPI[t - 1] / CPI[t]        price = R[t]        units_d += MONTHLY / price        cash_t += MONTHLY        if LR[t] <= TL[t]:            units_t += cash_t / price; cash_t = 0.0    return units_d * R[end], units_t * R[end] + cash_troll = {}for years in (10, 20, 30, 40):    h = years * 12    for mode in ('zero', 'yield'):        wins = tot = 0; ratios = []        for st in range(first_signal, n - h):            d, tr = final_values(st, st + h - 1, mode)            tot += 1; wins += tr > d; ratios.append(tr / d)        ratios.sort()        roll[f'{years}y_{mode}'] = dict(windows=tot, win_rate=wins / tot, median_ratio=ratios[len(ratios)//2],                                       worst_ratio=ratios[0], best_ratio=ratios[-1])out['rolling'] = roll# series for the charts (zero-cash base case), real index rebased to 100 at the start of the 100 year windowbase = R[start100]out['series'] = dict(    dates=[dates[t] for t in range(start100, end + 1)],    index=[round(R[t] / base * 100, 2) for t in range(start100, end + 1)],    trend_line=[round(math.exp(TL[t]) / base * 100, 2) for t in range(start100, end + 1)],    below=[1 if LR[t] <= TL[t] else 0 for t in range(start100, end + 1)],)zc = run(start100, end, 'zero')[0]out['series']['contributed'] = [round(x[1]) for x in zc]out['series']['dca'] = [round(x[2]) for x in zc]out['series']['trend'] = [round(x[3]) for x in zc]yc = run(start100, end, 'yield')[0]out['series']['ratio_zero'] = [round(x[3] / x[2], 4) for x in zc]out['series']['ratio_yield'] = [round(x[3] / x[2], 4) for x in yc]out.pop('curve', None)sig = out['series']['below']out['signal'] = dict(share_of_months=sum(sig) / len(sig))longest = cur = 0for x in zc:    cur = cur + 1 if x[4] > 0 else 0    longest = max(longest, cur)out['signal']['longest_cash_stretch_months'] = longestout['signal']['peak_cash_in_paychecks'] = max(x[4] for x in zc) / PAYCHECK# the same test started after the 1929 to 1932 crash, to show how much of the result rests on iti33 = dates.index('1933-01')d0, t0 = final_values(i33, end, 'zero'); d1, t1 = final_values(i33, end, 'yield')out['from_1933'] = dict(start='1933-01', ratio_zero=t0 / d0, ratio_yield=t1 / d1)# before inflation: the same paychecks in the dollars of their own day, so the# money weighted return includes inflation (the gap to the real return is inflation)def nominal_irr(values_end, start, stop):    flows = [-MONTHLY * CPI[t] / CPI[stop] for t in range(start, stop + 1)]    flows[-1] += values_end    return irr(flows)out['nominal'] = dict(    dca_irr=nominal_irr(out['zero']['dca'], start100, end),    trend_irr=nominal_irr(out['zero']['trend'], start100, end),    inflation=(CPI[end] / CPI[start100]) ** (12 / (end - start100)) - 1,    contributed_nominal=sum(MONTHLY * CPI[t] / CPI[end] for t in range(start100, end + 1)),) # the thesis: what bi-weekly buying of one paycheck amount became, after inflation,# over a working life, from every start month since 1871 (no trend line needed)def dca_final(start, stop):    units = 0.0    for t in range(start, stop + 1):        units += MONTHLY / R[t]    return units * R[stop]hz = {}for years in (10, 20, 30, 40):    h = years * 12    finals = sorted(dca_final(st, st + h - 1) for st in range(0, n - h + 1))    put_in = MONTHLY * h    hz[str(years)] = dict(windows=len(finals), put_in=put_in, worst=finals[0], p10=finals[len(finals) // 10],                          median=finals[len(finals) // 2], best=finals[-1],                          below_put_in=sum(f < put_in for f in finals) / len(finals))out['horizons'] = hzjson.dump(out, open(sys.argv[2], 'w'), separators=(',', ':'))for k in ('zero', 'yield'):    print(k, {a: (round(b, 4) if isinstance(b, float) else b) for a, b in out[k].items()})for k, v in roll.items(): print(k, {a: round(b, 3) for a, b in v.items()})print(out['meta'])