The Invariance Quest & Horizon Projection — the R engine¶

The independent R implementation, in base R only (acf, Box.test, fft). It reproduces the two foundational steps: find the i.i.d. invariant, then project it to the horizon by convolution — and confirms that the naive "square-root-of-time" VaR is too optimistic. Self-contained; no prior reading required.

The rule: model the i.i.d. increment (return, yield change), never the persistent level; then the horizon quantity is the sum of $T$ independent copies, whose distribution is the $T$-fold self-convolution.

In [1]:
options(repr.plot.width=11, repr.plot.height=4.6)
BLUE<-"#2b6cb0"; ORANGE<-"#dd6b20"; GREEN<-"#2f855a"; RED<-"#c53030"; GREY<-"#718096"

D <- read.csv("market_levels_daily.csv", row.names=1)
spy <- D$SPY; tnx <- D$TNX; vix <- D$VIX
ret <- diff(log(spy))*100                        # SPY daily log-return (the invariant), %
cat(sprintf("Daily market levels: %d days, columns %s\n", nrow(D), paste(colnames(D), collapse=", ")))

Daily market levels: 3771 days, columns SPY, TNX, VIX

The dataset¶

Daily levels of three markets, 2010–2024 (3,771 days): SPY (equity price → invariant = log-return), TNX (10-yr yield → invariant = yield change), VIX (implied vol → invariant = log-change). Levels are persistent; the increments are i.i.d.

1. The invariance quest¶

An invariant is i.i.d.: near-zero autocorrelation and a distribution stable over time. The lag-1 autocorrelation is the quick test — $\approx 1$ for a persistent level, $\approx 0$ for the invariant.

In [2]:
ac1 <- function(x) acf(x, lag.max=1, plot=FALSE)$acf[2]
series <- list("SPY: log-price"=log(spy), "SPY: log-return"=ret,
               "TNX: yield level"=tnx, "TNX: yield change"=diff(tnx),
               "VIX: level"=vix, "VIX: log-change"=diff(log(vix))*100)
cat("Lag-1 autocorrelation (level ~1 persistent; increment ~0 invariant):\n\n")
for (nm in names(series)) {
  a <- ac1(series[[nm]]); tag <- if (abs(a) < 0.2) "  <- INVARIANT" else "  (not invariant)"
  cat(sprintf("  %-20s ac(1) = %+.3f%s\n", nm, a, tag))
}
# visual: price acf (persistent) vs return acf (flat)
par(mfrow=c(1,2))
plot(0:20, acf(log(spy), lag.max=20, plot=FALSE)$acf, type="h", lwd=3, col=RED, ylim=c(-.3,1.05),
     xlab="lag", ylab="autocorrelation", main="SPY log-PRICE (persistent)"); abline(h=0)
plot(0:20, acf(ret, lag.max=20, plot=FALSE)$acf, type="h", lwd=3, col=BLUE, ylim=c(-.3,1.05),
     xlab="lag", ylab="autocorrelation", main="SPY log-RETURN (invariant)"); abline(h=0)
par(mfrow=c(1,1))
lb <- Box.test(ret, lag=20, type="Ljung-Box")
# Box.test computes 1 - pchisq(Q, df), which cancels to exactly 0 for tails this small;
# take the upper tail directly so the p-value survives in double precision.
lb_p <- pchisq(lb$statistic, 20, lower.tail=FALSE)
cat(sprintf("\nLjung-Box(20) on returns: Q = %.0f on 20 df, p = %.2g -- Q is NOT tiny: returns show mild\n",
            lb$statistic, lb_p))
cat("volatility clustering, the honest caveat that motivates GARCH; returns are only APPROXIMATELY i.i.d.\n")

Lag-1 autocorrelation (level ~1 persistent; 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

Ljung-Box(20) on returns: Q = 217 on 20 df, p = 4.3e-35 -- Q is NOT tiny: returns show mild

volatility clustering, the honest caveat that motivates GARCH; returns are only APPROXIMATELY i.i.d.
No description has been provided for this image

2. Horizon projection by convolution¶

Since the invariant is i.i.d., the $T$-day return is the sum of $T$ copies, whose density is the $T$-fold self-convolution — computed exactly with the FFT: $\mathcal F^{-1}[\mathcal F[p]^T]$. We check it against simulation, and watch the distribution drift toward normal (the central limit theorem).

In [3]:
moms <- function(x){ m<-mean(x); s<-sd(x); z<-(x-m)/s
  list(mean=m, sd=s, skew=mean(z^3), exkurt=mean(z^4)-3) }
project_fft <- function(x, T, bins=1500){
  lo<-min(x); hi<-max(x); dx<-(hi-lo)/(bins-1)
  h<-hist(x, breaks=seq(lo-dx/2, hi+dx/2, length.out=bins+1), plot=FALSE)$counts
  p<-h/sum(h); L<-T*(bins-1)+1
  ph<-fft(c(p, rep(0, L-length(p))))
  dens<-pmax(Re(fft(ph^T, inverse=TRUE))/L, 0)
  grid<-T*lo + (0:(L-1))*dx; dens<-dens/(sum(dens)*dx)
  list(grid=grid, dens=dens)
}
sim_h <- function(x, T, n) rowSums(matrix(sample(x, n*T, replace=TRUE), n, T))

m <- moms(ret)
cat(sprintf("Daily invariant: mean %.3f  sd %.3f  skew %.2f  excess-kurtosis %.1f\n\n", m$mean, m$sd, m$skew, m$exkurt))
set.seed(1); par(mfrow=c(1,3))
for (T in c(1, 21, 252)) {
  f <- project_fft(ret, T); s <- sim_h(ret, T, 100000)
  muT <- T*m$mean; sdT <- sqrt(T)*m$sd
  hist(s, breaks=100, freq=FALSE, col=rgb(.87,.42,.13,.35), border=NA, xlim=c(muT-4*sdT, muT+4*sdT),
       main=sprintf("horizon = %d day(s)", T), xlab="cumulative return (%)")
  lines(f$grid, f$dens, col=BLUE, lwd=2)
  curve(dnorm(x, muT, sdT), add=TRUE, col=GREY, lwd=2, lty=2)
}
par(mfrow=c(1,1))
cat("FFT convolution (blue) matches simulation (bars); the normal (dashed) fits only at long horizons.\n")

Daily invariant: mean 0.051  sd 1.078  skew -0.72  excess-kurtosis 11.5


FFT convolution (blue) matches simulation (bars); the normal (dashed) fits only at long horizons.
No description has been provided for this image

3. The central-limit drift and square-root-of-time VaR¶

Skewness decays like $1/\sqrt T$ and excess kurtosis like $1/T$: the distribution becomes normal, slowly. The square-root-of-time rule assumes normality immediately, so it understates the fat left tail — and hence the VaR — at short-to-medium horizons.

In [4]:
cf_z <- function(a, sk, ek){ z<-qnorm(a); z + (z^2-1)*sk/6 + (z^3-3*z)*ek/24 - (2*z^3-5*z)*sk^2/36 }
var_sqrt <- function(T, a=0.01) T*m$mean + qnorm(a)*sqrt(T)*m$sd
var_fft  <- function(T, a=0.01){ f<-project_fft(ret,T); cdf<-cumsum(f$dens)*(f$grid[2]-f$grid[1])
  suppressWarnings(approx(cdf, f$grid, a)$y) }
Ts <- c(1,2,5,10,21,42,63,126,252)
vs <- sapply(Ts, var_sqrt); vf <- sapply(Ts, var_fft)
par(mfrow=c(1,2))
plot(Ts, m$skew/sqrt(Ts), type="b", pch=19, col=BLUE, log="x", xlab="horizon T (days)", ylab="moment",
     main="Non-normality decays (CLT drift)", ylim=range(c(m$skew/sqrt(Ts), m$exkurt/Ts)))
lines(Ts, m$exkurt/Ts, type="b", pch=15, col=RED); abline(h=0, col=GREY)
legend("topright", c("skewness ~1/sqrt(T)","excess kurtosis ~1/T"), col=c(BLUE,RED), pch=c(19,15), bty="n")
plot(Ts, vs, type="b", pch=15, col=GREY, log="x", xlab="horizon T (days)", ylab="1% VaR (cumulative %)",
     main="Square-root rule too optimistic", ylim=range(c(vs,vf)))
lines(Ts, vf, type="b", pch=19, col=BLUE)
legend("bottomleft", c("square-root-of-time","exact projection"), col=c(GREY,BLUE), pch=c(15,19), bty="n")
par(mfrow=c(1,1))
for (T in c(5,21,63)) cat(sprintf("  T=%3dd: sqrt-rule VaR %.2f%%  exact %.2f%%  -> understated by %.2f%%\n",
                                   T, var_sqrt(T), var_fft(T), var_fft(T)-var_sqrt(T)))

  T=  5d: sqrt-rule VaR -5.35%  exact -6.23%  -> understated by -0.88%
  T= 21d: sqrt-rule VaR -10.42%  exact -11.46%  -> understated by -1.04%
  T= 63d: sqrt-rule VaR -16.69%  exact -17.57%  -> understated by -0.88%
No description has been provided for this image

4. Summary¶

  • The invariance quest, in base R: levels are persistent (autocorrelation $\approx1$), increments are invariant (autocorrelation $\approx0$) — across equity, fixed income, and volatility.
  • Returns are only approximately i.i.d. (mild volatility clustering, non-trivial Ljung-Box) — the honest caveat behind GARCH/SV.
  • Horizon projection by FFT convolution matches simulation; non-normality decays (skew $\sim1/\sqrt T$, kurtosis $\sim1/T$), so the square-root-of-time VaR understates risk at short-to-medium horizons.

Identical conclusions to the Python and PyMC engines — the foundational two-step recipe (find the invariant, project to the horizon) that underpins the entire Risk and Asset Allocation arc.