The Invariance Quest & Horizon Projection

The market-modelling foundation of the Risk and Asset Allocation arc — from scratch

"You cannot forecast a stock price or a bond yield — they wander and never repeat. You can only model the small, repeating quantity hidden inside them, and then add it up to your horizon."

Self-contained; no prior reading of Meucci's Risk and Asset Allocation (Springer 2005, ch. 3) required.

Two steps every risk model rests on

Before any of the estimation and allocation in this arc can even begin, two questions must be answered — and they are so fundamental they are usually skipped:

1. The invariance quest. What do we actually model? Not the stock price (it trends upward forever and today's price is almost exactly yesterday's — perfectly persistent), but its return. Not the bond yield, but its change. Not the VIX, but its log-change. The object we can model is the invariant: the quantity whose distribution repeats identically and independently through time (i.i.d.). Identifying it is step one, and getting it wrong dooms everything downstream.

2. Horizon projection. We estimate the invariant at a short base horizon — daily returns, say — but we invest over a long investment horizon — a month, a year. Because the invariant is i.i.d., the horizon quantity is simply the sum of $T$ independent copies. Its distribution is the $T$-fold self-convolution: its mean and variance grow linearly, and by the central limit theorem it drifts toward normal — but slowly enough that at real horizons it keeps the skew and fat tails that the industry-standard "square-root-of-time" rule throws away.

Roadmap

1. The invariance quest across three markets — equity, fixed income, volatility.
2. Horizon projection of the whole distribution — three methods that agree.
3. The central-limit drift and why "square-root-of-time" VaR is too optimistic.
4. From invariant to price and P&L at the horizon.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import horizon as hz

plt.rcParams.update({"figure.figsize": (9,4.5), "axes.grid": True, "grid.alpha": .25,
                     "axes.spines.top": False, "axes.spines.right": False, "font.size": 11})
BLUE, ORANGE, GREEN, RED, GREY, PURP = "#2b6cb0","#dd6b20","#2f855a","#c53030","#718096","#6b46c1"
np.set_printoptions(precision=3, suppress=True)

D = pd.read_csv("market_levels_daily.csv", index_col=0, parse_dates=True)
print("Daily market levels:", D.shape, "->", list(D.columns), "\n", D.tail(3).round(2))

Daily market levels: (3771, 3) -> ['SPY', 'TNX', 'VIX'] 
                SPY   TNX    VIX
date                           
2024-12-26  591.22  4.58  14.73
2024-12-27  584.99  4.62  15.95
2024-12-30  578.32  4.55  17.40

The dataset

Daily levels of three markets, 2010–2024 (3,771 trading days), so that we can hunt for each one's invariant:

Series	Market	Level (not invariant)	Invariant (i.i.d.)
SPY	US equity	price	log-return $\log(P_t/P_{t-1})$
TNX	10-yr Treasury	yield (%)	yield change $y_t-y_{t-1}$
VIX	equity implied vol	index level	log-change $\log(V_t/V_{t-1})$

The rule of thumb: prices/levels compound or trend and are highly persistent; the right invariant is the increment (return, change, or log-change) that a random walk adds each period.

1. The invariance quest

An invariant must be i.i.d.: (a) independent — no autocorrelation, so past values don't predict future ones; and (b) identically distributed — its distribution is stable across sub-samples. We test both. The clearest single diagnostic is the autocorrelation at lag 1: a persistent level has $\approx 1$, an invariant has $\approx 0$.

spy, tnx, vix = D["SPY"].values, D["TNX"].values, D["VIX"].values
series = {
    "SPY: log-price":   np.log(spy),                 "SPY: log-return":  np.diff(np.log(spy))*100,
    "TNX: yield level": tnx,                         "TNX: yield change":np.diff(tnx),
    "VIX: level":       vix,                          "VIX: log-change":  np.diff(np.log(vix))*100,
}
print("Lag-1 autocorrelation  (level ~1 = persistent, NOT invariant;  increment ~0 = invariant):\n")
for name, s in series.items():
    tag = "  <- INVARIANT" if hz.autocorr(s,1)[1] < 0.2 else "  (not invariant)"
    print("  %-20s ac(1) = %+.3f%s" % (name, hz.autocorr(s, 1)[1], tag))

Lag-1 autocorrelation  (level ~1 = persistent, NOT invariant;  increment ~0 = invariant):

  SPY: log-price       ac(1) = +0.999  (not invariant)
  SPY: log-return      ac(1) = -0.099  <- INVARIANT
  TNX: yield level     ac(1) = +0.997  (not invariant)
  TNX: yield change    ac(1) = -0.015  <- INVARIANT
  VIX: level           ac(1) = +0.965  (not invariant)
  VIX: log-change      ac(1) = -0.071  <- INVARIANT

# Visual proof for equity: the price is a smooth trend; the return is a structureless cloud
lp = np.log(spy); ret = np.diff(lp)*100
fig, ax = plt.subplots(2, 3, figsize=(14, 7))
ax[0,0].plot(lp, color=RED, lw=.7); ax[0,0].set_title("SPY log-PRICE (persistent, trending)")
ax[0,1].bar(range(21), hz.autocorr(lp,20), color=RED); ax[0,1].set_title("autocorrelation of log-price"); ax[0,1].set_ylim(-.3,1.05)
ax[0,2].scatter(lp[:-1], lp[1:], s=3, alpha=.3, color=RED); ax[0,2].set_title("$x_t$ vs $x_{t-1}$: a line (dependent)")
ax[1,0].plot(ret, color=BLUE, lw=.5); ax[1,0].set_title("SPY log-RETURN (the invariant)")
ax[1,1].bar(range(21), hz.autocorr(ret,20), color=BLUE); ax[1,1].set_title("autocorrelation of returns"); ax[1,1].set_ylim(-.3,1.05)
ax[1,2].scatter(ret[:-1], ret[1:], s=3, alpha=.3, color=BLUE); ax[1,2].set_title("$x_t$ vs $x_{t-1}$: a cloud (independent)")
plt.tight_layout(); plt.show()
print("The log-price forms a line under its own lag (each value ~ the last); the return forms a")
print("featureless cloud. Only the cloud can be modelled as i.i.d. draws from a fixed distribution.")

The log-price forms a line under its own lag (each value ~ the last); the return forms a
featureless cloud. Only the cloud can be modelled as i.i.d. draws from a fixed distribution.

The same test, for all three markets

The panel above is equity only. The line-versus-cloud contrast is the argument, so it is worth seeing it hold where the asset class changes — a yield and a volatility index behave nothing like a share price, and the invariant is a different transform in each case.

# The lag scatter for every market: level on top (a line), increment below (a cloud)
mkts = [("SPY", np.log(spy),  np.diff(np.log(spy))*100, "log-price", "log-return %",  RED),
        ("TNX", tnx,          np.diff(tnx),             "yield level", "yield change", ORANGE),
        ("VIX", vix,          np.diff(np.log(vix))*100, "level",       "log-change %",  PURP)]

fig, ax = plt.subplots(2, 3, figsize=(14, 7.2))
for j, (nm, lev, inc, lname, iname, col) in enumerate(mkts):
    a0 = hz.autocorr(lev, 1)[1]; a1 = hz.autocorr(inc, 1)[1]
    ax[0, j].scatter(lev[:-1], lev[1:], s=3, alpha=.28, color=col, lw=0)
    ax[0, j].set_title(f"{nm} {lname}   ac(1) = {a0:+.3f}", fontsize=10)
    ax[1, j].scatter(inc[:-1], inc[1:], s=3, alpha=.28, color=col, lw=0)
    ax[1, j].set_title(f"{nm} {iname}   ac(1) = {a1:+.3f}", fontsize=10)
    for r in (0, 1):
        ax[r, j].set_xlabel("$x_{t-1}$"); ax[r, j].set_ylabel("$x_t$")
ax[0, 0].annotate("LEVELS: a line", (.04, .90), xycoords="axes fraction", fontsize=11,
                  weight="bold", color=GREY)
ax[1, 0].annotate("INCREMENTS: a cloud", (.04, .90), xycoords="axes fraction", fontsize=11,
                  weight="bold", color=GREY)
plt.tight_layout(); plt.show()

print("Three asset classes, three different transforms -- a log-difference for the equity index,")
print("a plain difference for the yield, a log-difference for the volatility index -- and the same")
print("verdict in all three. The level is its own best predictor and cannot be treated as repeated")
print("draws; the increment carries no such memory and can. The transform is chosen to suit the")
print("series (a yield is already a rate, so it is differenced rather than log-differenced), but")
print("the test that decides whether it worked is identical.")

Three asset classes, three different transforms -- a log-difference for the equity index,
a plain difference for the yield, a log-difference for the volatility index -- and the same
verdict in all three. The level is its own best predictor and cannot be treated as repeated
draws; the increment carries no such memory and can. The transform is chosen to suit the
series (a yield is already a rate, so it is differenced rather than log-differenced), but
the test that decides whether it worked is identical.

# Identically distributed? Compare the return distribution across the two halves of the sample
h1, h2 = ret[:len(ret)//2], ret[len(ret)//2:]
fig, ax = plt.subplots(1, 2, figsize=(12, 4.2))
ax[0].hist(h1, bins=80, density=True, alpha=.5, color=BLUE, label="first half")
ax[0].hist(h2, bins=80, density=True, alpha=.5, color=ORANGE, label="second half")
ax[0].set_xlim(-6,6); ax[0].set_title("SPY return distribution is stable across sub-samples"); ax[0].legend()
sm = hz.split_sample_moments(ret, 4)
ax[1].plot(sm[:,0], "o-", color=BLUE, label="mean"); ax[1].plot(sm[:,1], "s-", color=ORANGE, label="std")
ax[1].set_xticks(range(4)); ax[1].set_xticklabels(["Q1","Q2","Q3","Q4"]); ax[1].set_title("moments across 4 blocks"); ax[1].legend()
plt.tight_layout(); plt.show()
from scipy.stats import chi2 as _chi2
_q = hz.ljung_box(ret, 20); _p = float(_chi2.sf(_q, 20))
print(f"Ljung-Box(20) on returns: Q = {_q:.0f} on 20 df (p = {_p:.2g}) -- Q is NOT tiny: returns show")
print("mild dependence (volatility clustering), which is exactly why GARCH / stochastic-vol models")
print("exist. Returns are APPROXIMATELY i.i.d. -- good enough to be the invariant, and the honest")
print("caveat behind treating them as such.")

Ljung-Box(20) on returns: Q = 217 on 20 df (p = 4.3e-35) -- Q is NOT tiny: returns show
mild dependence (volatility clustering), which is exactly why GARCH / stochastic-vol models
exist. Returns are APPROXIMATELY i.i.d. -- good enough to be the invariant, and the honest
caveat behind treating them as such.

2. Horizon projection: adding the invariant up

Having found the invariant (SPY daily log-return), we now ask what it looks like over an investment horizon of $T$ days. Since it is i.i.d., the horizon log-return is the sum of $T$ independent copies, and there are three ways to get its distribution — which should all agree:

• Exact (characteristic function). The sum's density is the $T$-fold self-convolution of the invariant's density, computed in one shot as $\mathcal F^{-1}\!\big[\mathcal F[p]^{T}\big]$ (an FFT). No distributional assumption.
• Moments + Cornish–Fisher. The horizon moments follow simple rules — $\text{mean}\!\to\!T\mu$, $\text{var}\!\to\!T\sigma^2$, $\text{skew}\!\to\!\text{skew}/\sqrt T$, $\text{ex-kurt}\!\to\!\text{ex-kurt}/T$ — and a Cornish–Fisher expansion turns them into quantiles.
• Simulation. Bootstrap: sum $T$ resampled daily returns, many times.

m = hz.sample_moments(ret)
print("Daily invariant: mean %.3f%%  std %.3f%%  skew %.2f  excess-kurtosis %.1f (fat-tailed, left-skewed)\n"
      % (m["mean"], m["std"], m["skew"], m["exkurt"]))
rng = np.random.default_rng(0)
fig, ax = plt.subplots(1, 3, figsize=(15, 4.2))
for a, T, lab in zip(ax, [1, 21, 252], ["1 day", "1 month (21d)", "1 year (252d)"]):
    g, d = hz.project_density_fft(ret, T, bins=1500)
    sim = hz.simulate_horizon(ret, T, 150000, rng)
    a.plot(g, d, color=BLUE, lw=2, label="exact (FFT convolution)")
    a.hist(sim, bins=120, density=True, alpha=.35, color=ORANGE, label="simulation")
    mu, sd = T*m["mean"], np.sqrt(T)*m["std"]
    xx = np.linspace(g.min(), g.max(), 300)
    a.plot(xx, np.exp(-0.5*((xx-mu)/sd)**2)/(sd*np.sqrt(2*np.pi)), "--", color=GREY, label="normal")
    a.set_xlim(mu-4*sd, mu+4*sd); a.set_title("horizon = %s" % lab); a.set_xlabel("cumulative return (%)")
ax[0].legend(fontsize=8)
plt.suptitle("Projecting the invariant: exact convolution vs simulation vs normal", y=1.03)
plt.tight_layout(); plt.show()
print("All three methods agree. At 1 day the distribution is sharply peaked and fat-tailed; by 1 year")
print("it is nearly normal -- the central-limit drift. The interesting risk lives in between.")

Daily invariant: mean 0.051%  std 1.078%  skew -0.72  excess-kurtosis 11.5 (fat-tailed, left-skewed)

All three methods agree. At 1 day the distribution is sharply peaked and fat-tailed; by 1 year
it is nearly normal -- the central-limit drift. The interesting risk lives in between.

3. The central-limit drift, and why "square-root-of-time" is too optimistic

The horizon moments show the CLT precisely: skewness decays like $1/\sqrt T$ and excess kurtosis like $1/T$. The distribution becomes normal — eventually. The industry shortcut assumes it is normal immediately and just scales volatility by $\sqrt T$ (the "square-root-of-time" rule). That ignores the fat left tail at short-to-medium horizons and therefore understates risk exactly when it matters. We compare the 1%-VaR from the square-root rule against the exact projection across horizons.

Ts = np.array([1,2,5,10,21,42,63,126,252])
skew_T = [m["skew"]/np.sqrt(T) for T in Ts]; kurt_T = [m["exkurt"]/T for T in Ts]
var_sqrt = [hz.var_sqrt_rule(m, T, 0.01) for T in Ts]
var_cf   = [hz.var_cornish_fisher(m, T, 0.01) for T in Ts]
var_fft  = [hz.var_from_density(*hz.project_density_fft(ret, T, 1500), 0.01) for T in Ts]

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].plot(Ts, skew_T, "o-", color=BLUE, label="skewness $\\sim 1/\\sqrt{T}$")
ax[0].plot(Ts, kurt_T, "s-", color=RED,  label="excess kurtosis $\\sim 1/T$")
ax[0].axhline(0, color=GREY, lw=.8); ax[0].set_xscale("log"); ax[0].set_xlabel("horizon T (days)")
ax[0].set_title("Non-normality decays with horizon (the CLT drift)"); ax[0].legend()
ax[1].plot(Ts, var_sqrt, "s-", color=GREY, label="square-root-of-time (normal)")
ax[1].plot(Ts, var_fft,  "o-", color=BLUE, label="exact projection")
ax[1].plot(Ts, var_cf,   "^--", color=GREEN, label="Cornish-Fisher")
ax[1].set_xscale("log"); ax[1].set_xlabel("horizon T (days)"); ax[1].set_ylabel("1% VaR (cumulative %)")
ax[1].set_title("Square-root-of-time VaR is too optimistic at short horizons"); ax[1].legend()
plt.tight_layout(); plt.show()
for T in [5, 21, 63]:
    vs, vf = hz.var_sqrt_rule(m,T,0.01), hz.var_from_density(*hz.project_density_fft(ret,T,1500),0.01)
    print("  T=%3dd:  sqrt-rule VaR %.2f%%   exact VaR %.2f%%   -> naive rule understates loss by %.2f%%"
          % (T, vs, vf, vf - vs))

  T=  5d:  sqrt-rule VaR -5.35%   exact VaR -6.23%   -> naive rule understates loss by -0.88%
  T= 21d:  sqrt-rule VaR -10.42%   exact VaR -11.46%   -> naive rule understates loss by -1.04%
  T= 63d:  sqrt-rule VaR -16.69%   exact VaR -17.57%   -> naive rule understates loss by -0.88%

Do the three invariants project the same way?

They are all invariants, so they all obey the same convolution arithmetic. But they do not start from the same shape, and the central-limit drift is a statement about the starting shape decaying — so it is worth asking whether the square-root-of-time verdict reached on equity carries across.

# The CLT drift, and the sqrt-rule error, for each invariant
invs = [("SPY log-return", np.diff(np.log(spy))*100, RED),
        ("TNX yield change", np.diff(tnx),            ORANGE),
        ("VIX log-change",  np.diff(np.log(vix))*100, PURP)]
Ts3 = np.array([1,2,5,10,21,42,63,126,252])

fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))
print(f"{'invariant':18s} {'skew':>8} {'exkurt':>9}   {'sqrt-rule vs exact 1% VaR':>34}")
print("-" * 78)
for nm, v, col in invs:
    mm = hz.sample_moments(v)
    ax[0].plot(Ts3, [abs(mm["skew"])/np.sqrt(T) for T in Ts3], "o-", color=col, label=nm)
    ax[1].plot(Ts3, [mm["exkurt"]/T for T in Ts3], "s-", color=col, label=nm)
    gaps = []
    for T in Ts3:
        vs = hz.var_sqrt_rule(mm, T, 0.01)
        vf = hz.var_from_density(*hz.project_density_fft(v, T, 1500), 0.01)
        gaps.append(100*(vf - vs)/abs(vs))
    ax[2].plot(Ts3, gaps, "^-", color=col, label=nm)
    g5 = gaps[Ts3.tolist().index(5)]
    verdict = "understates the loss" if g5 < 0 else "OVERstates the loss"
    print(f"{nm:18s} {mm['skew']:8.3f} {mm['exkurt']:9.3f}   at T=5d: {g5:+6.1f}%  ({verdict})")

for a, t in zip(ax, ["|skewness| decays like $1/\\sqrt{T}$", "excess kurtosis decays like $1/T$",
                     "square-root rule error vs exact"]):
    a.set_xscale("log"); a.set_xlabel("horizon T (days)"); a.set_title(t, fontsize=10); a.legend(fontsize=8)
ax[2].axhline(0, color="k", lw=1)
ax[2].set_ylabel("% error in the 1% VaR")
plt.tight_layout(); plt.show()

print()
print("The three invariants are not variations on one shape -- they are different shapes, and the")
print("skewness column is where it matters. Equity returns are LEFT-skewed (-0.72): the bad tail is")
print("the long one, so a normal approximation cuts it off and the square-root rule understates the")
print("loss. VIX log-changes are RIGHT-skewed (+1.18) -- volatility spikes upward and drifts down --")
print("so their LEFT tail is thinner than normal and the same rule OVERstates the loss. Yield changes")
print("are almost symmetric (+0.02) and nearly normal already, so the rule is close to right.")
print()
print("That is worth stating plainly, because 'square-root-of-time is too optimistic' is the standard")
print("summary and it is only true for the standard case. The rule is not biased in a fixed direction;")
print("it is biased in whichever direction the invariant is skewed, and by roughly the amount the")
print("invariant departs from normality. On this data it errs by -16% for equity, -3% for yields and")
print("+7% for the volatility index at a one-week horizon.")
print()
print("The kurtosis panel shows the other half. All three decay like 1/T as promised, but from very")
print("different heights: equity starts at 11.5 excess kurtosis and needs ~23 days to fall below 0.5,")
print("VIX starts at 6.9 and needs ~14, while yield changes start at 1.8 and are there within a week.")
print("The horizon at which 'approximately normal' becomes safe is a property of the series, not a")
print("constant -- which is exactly why the projection step cannot be skipped.")

invariant              skew    exkurt            sqrt-rule vs exact 1% VaR
------------------------------------------------------------------------------

SPY log-return       -0.721    11.517   at T=5d:  -16.4%  (understates the loss)

TNX yield change      0.023     1.838   at T=5d:   -3.2%  (understates the loss)
VIX log-change        1.184     6.916   at T=5d:   +6.8%  (OVERstates the loss)

The three invariants are not variations on one shape -- they are different shapes, and the
skewness column is where it matters. Equity returns are LEFT-skewed (-0.72): the bad tail is
the long one, so a normal approximation cuts it off and the square-root rule understates the
loss. VIX log-changes are RIGHT-skewed (+1.18) -- volatility spikes upward and drifts down --
so their LEFT tail is thinner than normal and the same rule OVERstates the loss. Yield changes
are almost symmetric (+0.02) and nearly normal already, so the rule is close to right.

That is worth stating plainly, because 'square-root-of-time is too optimistic' is the standard
summary and it is only true for the standard case. The rule is not biased in a fixed direction;
it is biased in whichever direction the invariant is skewed, and by roughly the amount the
invariant departs from normality. On this data it errs by -16% for equity, -3% for yields and
+7% for the volatility index at a one-week horizon.

The kurtosis panel shows the other half. All three decay like 1/T as promised, but from very
different heights: equity starts at 11.5 excess kurtosis and needs ~23 days to fall below 0.5,
VIX starts at 6.9 and needs ~14, while yield changes start at 1.8 and are there within a week.
The horizon at which 'approximately normal' becomes safe is a property of the series, not a
constant -- which is exactly why the projection step cannot be skipped.

4. From invariant to price and P&L

The invariant and its projection live in return space; an investor cares about money. The final step ("pricing") maps the projected invariant back to the quantity of interest. For a stock held over $T$ days, the horizon simple return is $e^{S_T}-1$ where $S_T$ is the projected cumulative log-return, and the horizon price is $P_0\,e^{S_T}$. We push the exact projected distribution through this transform to get the honest distribution of the horizon P&L and its Value-at-Risk in dollars.

T = 252; P0 = spy[-1]
g, d = hz.project_density_fft(ret/100, T, bins=1500)          # back to raw log-return units
sim_logret = hz.simulate_horizon(ret/100, T, 200000, rng)
price_T = P0*np.exp(sim_logret); simple_ret = np.exp(sim_logret) - 1
var_ret = np.percentile(simple_ret, 1); var_price = P0*np.exp(np.percentile(sim_logret, 1))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ax[0].hist(simple_ret*100, bins=120, density=True, color=BLUE, alpha=.7)
ax[0].axvline(var_ret*100, color=RED, lw=2, label="1%% VaR = %.1f%%" % (var_ret*100))
ax[0].set_xlabel("1-year simple return (%)"); ax[0].set_title("Horizon return distribution"); ax[0].legend()
ax[1].hist(price_T, bins=120, density=True, color=GREEN, alpha=.7)
ax[1].axvline(P0, color="k", ls="--", label="today $%.0f" % P0)
ax[1].axvline(var_price, color=RED, lw=2, label="1%% VaR price $%.0f" % var_price)
ax[1].set_xlabel("SPY price in 1 year ($)"); ax[1].set_title("Horizon price distribution (lognormal-ish)"); ax[1].legend()
plt.tight_layout(); plt.show()
print("1-year outlook: expected simple return %.1f%%, but a 1-in-100 year loses %.1f%% (price to $%.0f)."
      % (100*simple_ret.mean(), -100*var_ret, var_price))
print("This is the full pipeline: identify the invariant -> project it to the horizon -> price the P&L.")

1-year outlook: expected simple return 15.4%, but a 1-in-100 year loses 24.1% (price to $439).
This is the full pipeline: identify the invariant -> project it to the horizon -> price the P&L.

5. Summary

• The invariance quest is step one of all market modelling: model the i.i.d. increment (return, yield change, log-vol change), never the persistent level. Autocorrelation and sub-sample stability tell them apart, cleanly, across equity, fixed income and volatility.
• An honest caveat surfaced by the tests: returns are only approximately i.i.d. — they cluster in volatility — which is precisely the opening for the GARCH and stochastic-volatility models elsewhere in the collection.
• Horizon projection adds the invariant up: the horizon distribution is the $T$-fold convolution, obtainable exactly by characteristic function, approximately by moments + Cornish–Fisher, or by simulation — all three agree.
• Non-normality decays (skew $\sim1/\sqrt T$, kurtosis $\sim1/T$), so the naive square-root-of-time VaR is too optimistic at short-to-medium horizons; the exact projection captures the fat left tail it misses.
• Pricing closes the loop: push the projected invariant through the payoff to get the horizon P&L and its risk in dollars.

This is the scaffolding beneath the whole arc: every estimator we built assumes an invariant, and every risk number implicitly lives at some horizon. Getting these two steps right is what makes the rest meaningful.