The Invariance Quest & Horizon Projection

Python · PyMC · R  ·  Download projection module

The Two Steps

The market-modelling foundation of the Risk and Asset Allocation arc, and two questions so basic they are usually skipped. What do we actually model? Not the stock price — it trends forever and today's value is almost exactly yesterday's — but its return. Not the bond yield, but its change. Not the VIX, but its log-change. The modellable object is the invariant: the quantity whose distribution repeats identically and independently through time. And over what period? We estimate the invariant daily but invest over months or years, so it has to be projected to the investment horizon. Get either step wrong and everything downstream inherits the error.

SeriesMarketLevel — persistentInvariant — i.i.d.
SPYUS equity price, ρ1=+0.999\rho_1=+0.999 log-return, ρ1=0.099\rho_1=-0.099
TNX10-year Treasury yield, ρ1=+0.997\rho_1=+0.997 yield change, ρ1=0.015\rho_1=-0.015
VIXimplied volatility index level, ρ1=+0.965\rho_1=+0.965 log-change, ρ1=0.071\rho_1=-0.071

Finding the invariant

The hunt is settled empirically rather than by assertion, across three different markets — equity, fixed income and volatility — using lag-1 autocorrelation as the discriminator. Levels are overwhelmingly persistent (SPY log-price +0.999+0.999, the 10-year yield +0.997+0.997, VIX +0.965+0.965); their increments are not (0.099-0.099, 0.015-0.015, 0.071-0.071). Plotted against its own lag, a log-price forms a straight line while the return forms a featureless cloud — and only a cloud can be treated as i.i.d. draws. The notebooks draw that lag scatter for all three markets rather than for equity alone, which is what makes the point an empirical finding rather than an equity habit: the transform differs by asset class (a log-difference for the index and for VIX, a plain difference for the yield, since a yield is already a rate) but the test that decides whether it worked is identical, and all three pass it.

Testing the assumption instead of assuming it

The notebooks then do something the textbook version usually skips: they test the assumption they depend on, and report that it fails. A Ljung-Box test on 20 lags of SPY returns gives Q=217Q=217 against 20 degrees of freedom, p4×1035p\approx4\times10^{-35} — decisive evidence of serial dependence. Returns are therefore only approximately i.i.d.: the dependence is in the volatility, not the mean, which is precisely the opening for the GARCH and stochastic-volatility models elsewhere in the collection. Stating that plainly is more useful than a clean assumption quietly assumed.

Projecting to the Horizon

Because the invariant is i.i.d., the horizon quantity is a sum of TT independent copies, and its distribution is the TT-fold self-convolution. The notebooks compute it three independent ways — exactly by characteristic function via FFT, approximately by projecting moments and inverting with Cornish–Fisher, and by direct simulation — and all three agree. What the picture shows is the central-limit drift: skewness decays like 1/T1/\sqrt{T} and excess kurtosis like 1/T1/T, so the sharply peaked, fat-tailed daily distribution becomes nearly normal by a year. The interesting risk lives in between.

XT=t=1TXtφXT(ω)=[φX(ω)]T,skewT1/2,  kurtT1X_{T}=\sum_{t=1}^{T}X_t\quad\Longrightarrow\quad \varphi_{X_T}(\omega)=\big[\varphi_X(\omega)\big]^{T},\qquad \text{skew}\sim T^{-1/2},\ \ \text{kurt}\sim T^{-1}

Why square-root-of-time is too optimistic

That decay is exactly why the industry-standard square-root-of-time rule misleads. Scaling a daily VaR by T\sqrt{T} implicitly assumes a normality that has not yet arrived, so it misprices whatever non-normality is still present at short and medium horizons. On the equity invariant that means it understates the loss by around a percentage point across the range: 5.35%-5.35\% against 6.23%-6.23\% at five days, 10.42%-10.42\% against 11.46%-11.46\% at a month, 16.69%-16.69\% against 17.57%-17.57\% at a quarter. Python and R agree on every one of those figures to the basis point. But the direction of the error is not a constant, and running the comparison on the other two invariants shows it flipping. The rule is biased whichever way the invariant is skewed. Equity returns are left-skewed (0.72-0.72) — the bad tail is the long one — so a normal approximation truncates it and the rule understates the loss, by 16%16\% at a one-week horizon. Yield changes are almost symmetric (+0.02+0.02) and nearly normal already, so the rule is close to right, off by 3%3\%. VIX log-changes are right-skewed (+1.18+1.18), because volatility spikes upward and drifts back down, so their left tail is thinner than normal and the same rule overstates the loss by 7%7\%. The familiar summary — square-root-of-time is too optimistic — turns out to be a statement about equities rather than about the rule. The kurtosis side tells the same story from the other direction: all three decay like 1/T1/T as advertised, but from very different heights, so the horizon at which approximately normal becomes safe is about 23 days for equity, 14 for VIX and under a week for yields. It is a property of the series, not a constant — which is why the projection step cannot be skipped.

Risk with error bars

A third notebook re-runs the whole pipeline in PyMC and adds what a plug-in estimate cannot express: uncertainty about the parameters themselves. Fitting a Student-tt to the daily invariant puts the tail parameter at ν2.8\nu\approx2.8 with a 94% interval of roughly 2.52.53.13.1 — the data simply cannot pin down tail thickness, so any risk number depending on the tail must inherit that. Re-fitting the AR(1) coefficient as a Bayesian check returns ϕ0.098\phi\approx-0.098, matching the sample autocorrelation and confirming the invariant is serially clean in the mean. The payoff comes at the horizon: the one-month VaR carries a 94% interval more than a percentage point wide, and the Bayesian predictive VaR (11.50%-11.50\%) is more extreme than the plug-in (11.43%-11.43\%), because integrating parameter uncertainty fattens the tail. Reporting a single plug-in number quietly overstates how much is actually known.

Notebooks

Downloads

Projection Module — Source Code

"""
horizon.py -- The invariance quest and horizon projection.

Backs the notebooks in  "The Invariance Quest & Horizon Projection".
Follows Meucci, "Risk and Asset Allocation" (Springer 2005), chapter 3-A/B.

THE TWO STEPS OF MARKET MODELLING
---------------------------------
1. INVARIANCE QUEST.  Before you can estimate anything you must find the quantity
   whose distribution repeats identically through time -- the "invariant".  You
   cannot model a stock PRICE (it trends and is perfectly persistent) but you can
   model its log-RETURN; you cannot model a bond YIELD but you can model its daily
   CHANGE.  An invariant is identically and independently distributed (i.i.d.):
   its autocorrelation is ~0 and its distribution is stable across sub-samples.

2. HORIZON PROJECTION.  You estimate the invariant at a short base horizon (say
   daily), but you invest over a long one (say a year).  Because the invariant is
   i.i.d., the horizon quantity is the SUM of T independent copies, whose
   distribution is the T-fold self-convolution.  Its moments scale in a precise
   way and, by the CLT, it drifts toward normal -- but at short-to-medium horizons
   it keeps the skew and fat tails that a naive "square-root-of-time" rule ignores.

CONVENTIONS
-----------
x : (n,) a candidate invariant (e.g. daily log-returns).  T : projection horizon
(number of base periods).  All in the invariant's own units.
"""

import numpy as np
from scipy.stats import skew, kurtosis, norm


# --------------------------------------------------------------------------- #
#  Step 1 -- the invariance quest                                               #
# --------------------------------------------------------------------------- #

def autocorr(x, lags=20):
    """Sample autocorrelation function up to `lags`. An invariant has acf ~ 0."""
    x = np.asarray(x, float) - np.mean(x)
    n = len(x); denom = np.dot(x, x)
    return np.array([1.0 if k == 0 else np.dot(x[:-k], x[k:]) / denom for k in range(lags + 1)])


def ljung_box(x, lags=20):
    """Ljung-Box Q statistic testing joint zero autocorrelation (independence)."""
    n = len(x); ac = autocorr(x, lags)[1:]
    q = n * (n + 2) * np.sum(ac ** 2 / (n - np.arange(1, lags + 1)))
    return q          # compare to chi-square(lags); large Q -> reject i.i.d.


def split_sample_moments(x, k=2):
    """Mean/std of x over k equal consecutive sub-samples. A stable (invariant)
    series has near-identical moments across blocks."""
    x = np.asarray(x, float); parts = np.array_split(x, k)
    return np.array([[p.mean(), p.std()] for p in parts])


# --------------------------------------------------------------------------- #
#  Moments and their horizon scaling                                            #
# --------------------------------------------------------------------------- #

def sample_moments(x):
    """(mean, std, skewness, excess-kurtosis) of the invariant."""
    x = np.asarray(x, float)
    return dict(mean=x.mean(), std=x.std(ddof=1),
                skew=float(skew(x)), exkurt=float(kurtosis(x)))   # kurtosis() is excess by default


def project_moments(m, T):
    """Moments of the sum of T i.i.d. copies of the invariant.
       mean_T = T*mean;  var_T = T*var;  skew_T = skew/sqrt(T);  exkurt_T = exkurt/T.
    The last two shrink toward 0 -- the central-limit drift to normality."""
    return dict(mean=T * m["mean"], std=np.sqrt(T) * m["std"],
                skew=m["skew"] / np.sqrt(T), exkurt=m["exkurt"] / T)


# --------------------------------------------------------------------------- #
#  Step 2 -- horizon projection of the whole distribution                       #
# --------------------------------------------------------------------------- #

def project_density_fft(x, T, bins=2000):
    """EXACT i.i.d. horizon distribution via the characteristic function.
    Bin the invariant into a density p, then the T-fold self-convolution is
    IFFT( FFT(p)^T ) -- the horizon quantity's density on a grid. Returns
    (grid, density)."""
    x = np.asarray(x, float)
    lo, hi = x.min(), x.max()
    dx = (hi - lo) / (bins - 1)
    p, _ = np.histogram(x, bins=bins, range=(lo - dx / 2, hi + dx / 2))
    p = p.astype(float); p /= p.sum()
    L = T * (bins - 1) + 1
    ph = np.fft.rfft(p, n=L)
    dens = np.fft.irfft(ph ** T, n=L)
    dens = np.maximum(dens, 0.0)
    grid = T * lo + np.arange(L) * dx
    dens /= (dens.sum() * dx)
    return grid, dens


def simulate_horizon(x, T, n, rng):
    """Bootstrap: n samples of the horizon sum = sum of T i.i.d. draws of x."""
    idx = rng.integers(0, len(x), size=(n, T))
    return np.asarray(x)[idx].sum(axis=1)


# --------------------------------------------------------------------------- #
#  Quantiles / VaR at the horizon                                               #
# --------------------------------------------------------------------------- #

def cornish_fisher_z(alpha, skew, exkurt):
    """Cornish-Fisher adjusted standard normal quantile for skew & excess
    kurtosis -- a moment-based correction to the Gaussian quantile."""
    z = norm.ppf(alpha)
    return (z + (z ** 2 - 1) * skew / 6
            + (z ** 3 - 3 * z) * exkurt / 24
            - (2 * z ** 3 - 5 * z) * skew ** 2 / 36)


def var_cornish_fisher(m, T, alpha=0.01):
    """Horizon VaR from projected moments + Cornish-Fisher (captures skew/kurtosis)."""
    mT = project_moments(m, T)
    zcf = cornish_fisher_z(alpha, mT["skew"], mT["exkurt"])
    return mT["mean"] + zcf * mT["std"]


def var_sqrt_rule(m, T, alpha=0.01):
    """The naive Gaussian 'square-root-of-time' VaR: normal, moments scaled."""
    return T * m["mean"] + norm.ppf(alpha) * np.sqrt(T) * m["std"]


def var_from_density(grid, dens, alpha=0.01):
    """Quantile of a density given on a grid (for the exact FFT projection)."""
    cdf = np.cumsum(dens) * (grid[1] - grid[0])
    return float(np.interp(alpha, cdf, grid))

References