Learning Statistics with Python
Chapter 6 · Learning Statistics with Python
Stationarity, ACF/PACF and ARIMA as the foundations — then where time-series models still earn their keep in quant trading: volatility (ARCH/GARCH/GJR, VaR, vol targeting), structural arbitrage (cointegration and pairs, Kalman hedge ratios, regime switching), and the walk-forward discipline that underpins every ML pipeline.
Prof. Xuhu Wan
ISOM, HKUST Business School · 2026 Edition
Why this matters
Nobody on a modern desk forecasts tomorrow’s return with ARIMA. Yet every risk limit, every hedge ratio, every regime flag and every feature that goes into the ML model is a time-series model. Knowing which jobs the classical toolkit still owns is the difference between a quant and a hobbyist.
Where time-series models sit in quant trading today
Time-series models remain important in quant trading, but their role has shifted: to volatility modelling, to structural arbitrage (cointegration, Kalman filters, regime switching), and to serving as the foundation of the ML pipeline. Mean-return forecasting itself has largely been taken over by ML.
| Job | Who owns it now | Why |
|---|---|---|
| Forecast tomorrow’s mean return | ML on the cross-section (Ch. 3) | one series has almost no linear memory; the edge is in many series and many features |
| Forecast tomorrow’s variance | GARCH family (§6.5–5.6) | variance is persistent, and a 15-line MLE tracks it |
| Trade a relationship, not a direction | cointegration, Kalman, regimes (§6.3–5.4) | the object traded is a stationary spread, not a price |
| Build the ML pipeline | stationarity, lags, walk-forward (§6.1–5.2) | leakage and non-stationary features kill more models than bad algorithms do |
| Section | Concept | Tool |
|---|---|---|
| §6.1 Foundations | Stationarity, unit roots, ACF/PACF, ARIMA baseline | kpss, adfuller, plot_acf, ARIMA |
| §6.2 Why mean forecasting lost to ML | Forecast scoring, walk-forward, TS features for GBM | .shift, .rolling, HistGradientBoostingRegressor |
| §6.3 Structural arbitrage I | Cointegration, pairs trading | coint, rolling hedge ratio, z-score rule |
| §6.4 Structural arbitrage II | Kalman filter, time-varying beta, regime switching | UnobservedComponents, MarkovRegression |
| §6.5 Volatility I | ARCH, GARCH(1,1), forecasts, VaR | het_arch, scipy.optimize, stats.t |
| §6.6 Volatility II | GJR leverage, ARMA+GARCH, volatility targeting | LR test, ARIMA(...).resid, target/σ̂ |
Data: Apple minute bars, AAPL, S&P 500, GOOG, Dow, NVDA/SPY, airline passengers, FRED-MD.
You will look at one trading day of Apple twice — as prices and as returns — and predict which one a model is allowed to touch. Then two tests with opposite nulls give you a verdict you can defend, the ACF tells you how much memory a series has, and AIC picks an ARIMA — which you will hold out and score.
390 one-minute bars for 2 November 2020. Load the close price and look at it.
389 rows after dropping the first return. The price opens near 109.6, climbs above 110.5, then drifts down to 108.8 by the close. Wherever it is at 10:00 tells you a lot about where it will be at 10:01 — the level has memory.
Which series looks stationary — constant mean and variance over the day?
Mean -0.000018, std 0.00108: a flat cloud around zero. The price is the cumulative sum of these; summing a stationary series is what manufactures a wandering level.
\(X_t = a\,X_{t-1} + \varepsilon_t\). With \(a = 1\), \(\text{Var}(X_t) = t\sigma^2\) grows forever — the random walk, a unit root. With \(|a|<1\) it settles at \(\sigma^2/(1-a^2)\). Feed identical shocks into \(a = 1\) and \(a = 0.8\) (seed 1, 200 steps). Predict which ends further from zero.
Same 200 shocks: the random walk wanders to -15.8 (std 5.27); the \(a=0.8\) series never leaves [-3.1, 3.1] (std 1.15). The difference between \(a = 0.99\) and \(a = 1\) is the difference between a model and a lottery.
Trend-stationary \(y_t = g(t) + X_t\): cure by detrending. Difference-stationary \(y_t = y_{t-1}+\varepsilon_t\): cure by differencing.
\(y_t = 0.5\,t + 3\varepsilon_t\). You difference it instead of detrending. What do you get?
Detrending recovers the shock scale (2.77 ≈ 3); differencing the same series gives 4.06 ≈ \(3\sqrt2\). Differencing the random walk gives 0.93: exactly the shocks. Match the cure to the disease.
kpss(x, regression="c") tests stationarity around a constant; "ct" around a linear trend. Small p ⇒ reject stationarity. Run it on the trend series with both nulls and predict which one rejects.
Around a constant: stat 2.102, p 0.01 — rejected (the trend is not a constant mean). Around a trend: stat 0.058, p 0.1 — not rejected. KPSS reports p only inside [0.01, 0.1]; 0.01 means “≤ 0.01”.
adfuller estimates \(\Delta y_t = \alpha + \beta t + \gamma\,y_{t-1} + \dots\) and tests \(\gamma = 0\) (i.e. \(a = 1\)). Small p ⇒ reject the unit root.
The trend series fails ADF with a constant (p 0.95) but rejects the unit root once a trend is allowed (p 0.0000). The random walk: ADF p 0.81, KPSS p 0.01 — both agree it is non-stationary. AR(0.8): ADF p 0.0000, KPSS p 0.1 — both agree it is stationary.
| ADF rejects (p < 0.05) | ADF does not reject | |
|---|---|---|
| KPSS does not reject (p ≥ 0.05) | Stationary — model it | Not enough data to tell (low power) |
| KPSS rejects (p < 0.05) | Conflict: trend-stationary or structural break — try "ct" |
Unit root — difference it |
Important
Always run both. ADF alone cannot distinguish “stationary” from “too little data to reject”; KPSS alone cannot distinguish “unit root” from “deterministic trend”. The pair, with "c" and "ct", covers the cases.
For the close price: KPSS(c) p = 0.01, ADF(c) p = 0.035. What is the honest verdict?
Return: KPSS p 0.1, ADF stat -25.07, p 0.0000 — unambiguous. The price is the conflict cell; the return is the “model it” cell.
FRED-MD, monthly 1959–2022. The first row holds transformation codes — drop it. Test the CPI level and the unemployment rate.
p_cpi is done. Compute p_unrate, the ADF p-value (constant only) for UNRATE, and compare. Which one is stationary in levels?
CPI p 0.9986 — a unit root (with a trend too: "ct" gives 0.77). UNRATE p 0.011 — stationary: unemployment goes up and comes back down. A feature with a unit root goes into a model only after differencing — remember this in §6.2.
\(\rho_k = \text{corr}(y_t, y_{t-k})\). The ACF plots \(\rho_k\) against \(k\); the PACF plots the correlation at lag \(k\) after removing lags \(1,\dots,k-1\).
For monthly airline passengers (144 months, 1949–1960), which is larger: the ACF at lag 6 or at lag 12?
ACF: 0.95, 0.68, 0.76 — slow decay (trend) with a bump at 12 (season). PACF: 0.95, -0.23, 0.04 — a single spike then nothing: once you know last month, the month before adds little. AR-like level, seasonal echo.
| Pattern | ACF | PACF | Model |
|---|---|---|---|
| AR(\(p\)) | decays geometrically | cuts off after lag \(p\) | past values |
| MA(\(q\)) | cuts off after lag \(q\) | decays | past errors |
| White noise | all inside \(\pm 2/\sqrt n\) | same | nothing to model |
Dow log returns 2015–2019: all ten autocorrelations inside ±0.056 — the mean has no memory. But \(|r|\) has ACF 0.28, 0.26, 0.27, … — volatility remembers. That second panel is §6.5. It is the whole chapter in one figure.
Difference \(d\) times, then fit ARMA(\(p,q\)): \(\Delta^d y_t = c + \sum\phi_i\,\Delta^d y_{t-i} + \varepsilon_t + \sum\theta_j\,\varepsilon_{t-j}\). Log prices need \(d=1\); returns are already stationary, \(d=0\).
\(p, q \in \{0,1,2\}\), \(d=0\), chosen by AIC on the Dow returns above. Which order wins?
Winner: (0, 0, 0), AIC -8409.8. That is not a failure of ARIMA; it is a fact about the market — and the reason §6.2 exists.
Train on log passengers 1949–1959 (132 months), hold out 1960. Grid \(p, q \in \{0,1,2\}\) with \(d = 1\).
Best by AIC: (2, 1, 1), AIC −228.9. RMSE over 1960: ARIMA(2,1,1) versus the seasonal-naive forecast (1959’s value for the same month). Who wins?
ARIMA RMSE 89.4, seasonal-naive 50.7. The forecast settles between 427 and 440 with bands widening from ±84 to ±188. The bands are honest; the model is blind to the season. Order selection is not validation.
SARIMA adds a seasonal difference and a seasonal MA at lag 12: order=(0,1,1), seasonal_order=(0,1,1,12) — the classic “airline model”.
RMSE 18.6 — almost five times better than the non-seasonal fit (89.4) and nearly three times better than the seasonal-naive. Same two-parameter budget, one extra idea. AIC agrees: -441.3 versus -228.9. Where a series has structure, the classical model is unbeatable per parameter.
Next: §6.2 — score a mean forecast honestly, and see what a gradient-boosting model does with the same series.
A forecast is only as good as the baseline it beats. You will score the laziest forecast of all and discover that “predict zero” beats it; then watch the AIC-winning ARIMA lose money out of sample. The pivot: what a modern ML pipeline inherits from time series — stationary features, lags without leakage, walk-forward validation — and what a gradient-boosting model honestly delivers on one series.
With actual \(y_t\) and forecast \(\hat y_t\) over \(n\) periods:
\[\text{MAE} = \tfrac1n\sum|y_t - \hat y_t|, \quad \text{RMSE} = \sqrt{\tfrac1n\sum (y_t-\hat y_t)^2}, \quad \text{MAPE} = \tfrac{100}{n}\sum\left|\tfrac{y_t-\hat y_t}{y_t}\right|,\]
\[\text{DA} = \tfrac1n\sum \mathbf 1\{\text{sign}(y_t) = \text{sign}(\hat y_t)\}.\]
Which metric is unusable for forecasting minute-by-minute price changes?
Target Y is the next minute’s price change; the forecast Yhat is the current change.
388 forecasts: MAE 0.1157 dollars, RMSE 0.1694, directional accuracy 0.466 — worse than a coin. Minute changes are slightly negatively autocorrelated (bid–ask bounce), so copying the last change points you the wrong way.
Forecast zero for every minute. The RMSE of the zero forecast is just the std of Y around zero. Is it above or below the naive RMSE of 0.1694?
0.1057
0.1057 vs 0.1694. Copying the last change adds the variance of two changes; predicting the mean (≈ 0) adds none. Any forecast you build must beat 0.1057, not the naive number. That zero forecast is the first baseline in every ML pipeline too.
\(\hat y_{t+1} = \tfrac{1}{40}\sum_{i=0}^{39} \Delta P_{t-i}\) — .rolling(40).mean(). Predict: does averaging out the bounce beat the zero forecast (0.1057)?
On the common 349 rows: SMA RMSE 0.0993, zero 0.0986, naive 0.1428. The SMA crushes naive and ties the zero forecast; its only edge is a direction call (DA 0.521) that a zero forecast cannot make. And compare only on the same rows — the naive RMSE moved from 0.1694 to 0.1428 just by dropping 40 warm-up minutes.
\(s_t = \alpha\,x_t + (1-\alpha)\,s_{t-1}\), \(s_1 = x_1\). As \(\alpha \to 1\) it is the naive forecast; as \(\alpha \to 0\) the long-run mean.
x = [10, 20, 30], α = 0.3, adjust=False. \(s_1 = 10\), \(s_2 = 0.3·20 + 0.7·10\), \(s_3 = 0.3·30 + 0.7·s_2\). What prints?
[10.0, 13.0, 18.1]
EWMA(0.3): RMSE 0.1075, DA 0.496 — between naive and SMA, because α = 0.3 remembers only ~2 minutes and inherits the bounce. The window is the model — and you tuned it by eye. Now let a likelihood choose.
Fit ARIMA(2,0,1) — the AIC winner — on AAPL log returns 2020–2023, then feed 2024 through res.apply(new_data) for one-step-ahead forecasts without refitting. Trade the sign.
DA 0.480, RMSE 0.01451 — worse than forecasting the training mean (0.01421); Sharpe -0.20 against 1.16 for holding. The AIC winner in-sample is a loser out of sample. This is the result that moved mean forecasting off the classical desk.
ML did not replace time series; it sits on top of it. Four things every return-prediction model in Chapter 4 quietly assumed:
| Inherited discipline | From | Failure if skipped |
|---|---|---|
| Stationary features — returns, ratios, differences, never price levels | §6.1 | tree splits on a level that never recurs; the model memorises 2021 |
| Lag to avoid leakage — every feature at \(t\) uses data to \(t-1\) | .shift(1) |
in-sample accuracy 0.9, live accuracy 0.5 |
| Walk-forward validation — expanding window, refit per period | §6.2 | random K-fold lets the future train the past |
| A baseline that must be beaten — zero forecast, ARIMA(0,0,0) | §6.1–5.2 | a Sharpe with nothing to compare it to |
Tip
The ML model is the last step of a time-series pipeline. Most of the value is in the first three rows.
Target is y = r (today’s return). A colleague builds X[“std21”] = r.rolling(21).std() with no shift. What is wrong?
Eleven features, all stationary (returns, rolling moments, a ratio), all shifted. The first fold trains on 2016–2017 and tests on 2018; each later fold expands the training window by a year. That is the walk-forward split.
Same spec as Chapter 4: HistGradientBoostingRegressor(max_depth=2, learning_rate=0.05, max_iter=200, min_samples_leaf=50). Refit every January on all data to date; predict the year; trade long/flat on the sign.
Honest numbers: DA 0.494, RMSE 0.01276 — worse than the zero forecast (0.01247), correlation -0.013; the rule’s Sharpe 0.31 against 0.57 for holding the index, in the market 59 % of days. No year’s DA clears 0.54. The same walk-forward discipline that exposed ARIMA exposes the GBM. ML did not win by forecasting one series from its own past — it won by moving to the cross-section and to richer features (Ch. 3.5). The pipeline is the inheritance; the alpha is elsewhere.
In Colab
import torch, torch.nn as nn
class LSTMReg(nn.Module):
def __init__(self, n_feat, hidden=32):
super().__init__(); self.lstm = nn.LSTM(n_feat, hidden, batch_first=True); self.head = nn.Linear(hidden, 1)
def forward(self, x): out, _ = self.lstm(x); return self.head(out[:, -1])
model = LSTMReg(n_feat=11); opt = torch.optim.Adam(model.parameters(), lr=1e-3)
for epoch in range(20): # X_seq: (batch, 21 days, 11 features), y: next-day return
opt.zero_grad(); loss = ((model(X_seq) - y).pow(2)).mean(); loss.backward(); opt.step()Your LSTM beats the GBM by 4 % accuracy in-sample. What do you check first?
\(X\) Granger-causes \(Y\) if lags of \(X\) improve the forecast of \(Y\) beyond \(Y\)’s own lags — an F-test on the \(X\) coefficients in \(Y\)’s equation (a VAR is \(k\) such regressions on shared lags). Forecastability, not mechanism.
p_nvda_to_spy is done. Compute p_spy_to_nvda — the p-value that SPY’s past helps forecast NVDA — using granger_p with 2 lags.
NVDA → SPY 0.6835, SPY → NVDA 0.1257. Neither direction at 5 %: in daily returns the lead–lag is priced away within the day. On macro data (Stock–Watson, FRED-MD) Granger tests are informative; on liquid daily prices they mostly confirm efficiency.
.shift(1) before anything else, expanding-window walk-forward, and a baseline to beat.Next: §6.3 — stop forecasting a direction; trade a stationary relationship instead.
Two drunks leave the pub and wander. Each path is a random walk. But if they are tied together by a rope, the distance between them is stationary. That rope is cointegration — and a stationary spread is the one thing in this chapter that a time-series model can trade directly. You will see two pairs with no rope, one with a rope, and learn why the rope is always tighter in-sample.
Regress Apple’s log price (2015–2024) on the Dow’s log price from 1985–1995, lined up row by row. The two series have nothing to do with each other. Predict the R².
Apple 2015–24 on Dow 1985–95, 2 516 rows each. What will R² and the t-statistic look like?
R² 0.804, t 101.7 — and an ADF p of 0.2138 on the residual: the “relationship” wanders. Rule: a regression between two unit-root series is only meaningful if the residual is stationary. That is the cointegration test.
statsmodels.tsa.stattools.coint(a, b) does steps 2–3 and returns (t_stat, p_value, critical_values).
You run a plain ADF on the OLS residual and get p = 0.024, but coint() reports p = 0.080. Which is right?
Both levels keep their unit roots (NVDA p 0.23, SPY p 0.79). β = 4.38; the naive ADF on the residual says p 0.024, but Engle–Granger says p 0.080. Not cointegrated at 5 %, borderline at 10 %: the rope is loose. The residual plot shows the drunks drifting apart for months at a time.
Estimate the rope on 2023 (formation), then watch the spread in 2024 (trading) — never the other way round.
The 2023 residual already fails (ADF p 0.32), and in 2024 the spread’s mean sits at -0.227 — 1.5 formation-sigmas below zero, and it never came back. A z-score would have fired 17 “entry” days on a spread with no anchor. No rope, no pairs trade. (AAPL vs the S&P over 2015–2024 fails too: coint p = 0.25.)
Five years of daily closes, 1 259 days. Predict: will coint reject at 5 %?
Both have unit roots (p 0.85, 0.80). β = 1.271, Engle–Granger t = -3.800 against a 5 % critical value of -3.341: p = 0.0136. Cointegrated at 5 %. A deviation halves in 24.3 days — slow, but it comes back. This is the rope.
The full-sample test rejects at p = 0.014. Split into formation (first 3 years) and trading (last 2). What do you expect on the formation window alone?
Formation to Oct 2017: p 0.239; to Oct 2018: p 0.113. Neither rejects at 5 %, and the trading-window spread drifts to a mean of 0.071 (1.8 formation-sigmas) in the first case. The rope you can see in 2020 is not the rope you could have tied in 2017. So: re-estimate it as you go.
Every 63 trading days, refit \((\alpha, \beta)\) on the trailing 252 days and hold them for the next 63. The spread uses only ropes that were known at the time; the z-score uses a 40-day rolling mean and std.
16 refits, 1 007 trading days from October 2016. The hedge ratio moves from 0.569 to 1.601 — the rope is not a constant. 110 days have |z| > 2: those are the candidates. Now the rule.
Short the spread (short GOOG, long β·SPX) when \(z > 2\); long it when \(z < -2\); close when \(|z|\) falls below 0.5. Yesterday’s position earns today’s spread change.
In a trade 35 % of days, 52 entries and exits. Pairs: cumulative 30.5 %, Sharpe 0.75, max drawdown -11.6 %. Holding GOOG: 68.6 %, Sharpe 0.63, drawdown -36.8 %. Less money, better money — and market-neutral: the book made it through March 2020 with a third of the drawdown. Before costs.
Change entry to 1.5 (keep exit at 0.5) and read off sh_rule. More trades — but is it a better book, or just a busier one?
Entry 1.5: in a trade 47 % of days, 80 round trips, cumulative 53.0 %, Sharpe 1.09, drawdown -11.6 %. Better on paper — with 80 crossings the transaction costs you ignored are now 50 % larger. One tuned threshold on one pair is a hypothesis, not a book.
coint, with its harsher critical values — makes a level regression meaningful.Next: §6.4 — the hedge ratio as a state that moves every day (Kalman), and the market’s regime as a hidden state (Markov switching).
A price is a hidden “true level” plus noise; a hedge ratio is a hidden slope that drifts; the market is in a hidden calm or stormy regime. You never observe any of them — but each new data point lets you update your best guess. Six lines of Kalman, eight lines of a time-varying beta, and one MarkovRegression call that cuts a drawdown by three quarters.
Local level model: \[\mu_t = \mu_{t-1} + \eta_t,\ \eta_t\sim N(0,Q) \qquad y_t = \mu_t + \varepsilon_t,\ \varepsilon_t \sim N(0,R).\]
Keep two numbers: the estimate \(\hat\mu_{t|t}\) and its variance \(P_{t|t}\).
Observations become much noisier (\(R\) ↑). The Kalman gain \(K_t\)…
Start \(\hat\mu = 0\), \(P = 1\), with \(Q = 0.5\), \(R = 1\). Observe \(y = 2\). Predict: \(P \to 1.5\); \(K = 1.5/2.5\); \(\hat\mu \to 0 + K\cdot 2\). What are K and mu (3 dp)?
0.6 1.2
[1.2]: 60 % of the surprise (\(2 - 0\)) is absorbed. New estimate = old prediction + gain × prediction error. Six lines; that is the whole filter.
sm.tsa.UnobservedComponents(y, level="local level") estimates \(Q\) (sigma2.level) and \(R\) (sigma2.irregular) by maximum likelihood, then filters.
\(\sqrt R\) = $0.048 — about a nickel of bid–ask bounce per minute; \(\sqrt Q\) = $0.096 — the true level moves twice as much. Steady-state gain ≈ 0.83: the filter follows the data closely because the level genuinely wanders. The six-line filter matches statsmodels to 0.005.
In §6.3 the rope moved from 0.57 to 1.60 and you re-tied it every quarter. A Kalman filter re-ties it every day: let \(\beta_t\) itself be a random walk.
\[\beta_t = \beta_{t-1} + \eta_t,\ \eta_t \sim N(0,Q) \qquad y_t = \beta_t\,x_t + \varepsilon_t,\ \varepsilon_t \sim N(0,R).\]
The only change: the observation is \(\beta_t x_t\), so the gain is \(K_t = P\,x_t/(x_t^2 P + R)\) and the update is \(\beta + K_t(y_t - \beta x_t)\).
Set \(Q = 0\). What does the filtered \(\beta_t\) become?
Daily returns in percent, 2023–2024. Static OLS gives β = 2.00 in 2023 and 2.64 in 2024 — one number per year. Set \(R\) to the OLS residual variance and the signal-to-noise ratio \(Q/R = 0.001\).
Quarter-end betas 1.72, 2.28, 2.05, 1.73, 2.66, 2.67, 3.18, 1.73 around a static 2.307: the filter sees the 2024 rise and the late-2024 fall that an annual OLS would report a year late. Honest number: the hedge error is 2.475 with yesterday’s Kalman beta versus 2.463 with the look-ahead static one and 3.068 unhedged. On 500 days of a single beta, the gain is timeliness, not variance — no window to choose, no refit schedule.
Calm and stormy are not two values of \(\sigma\) on a continuum; they are two states the market jumps between. Hamilton (1989): a Markov chain nobody observes drives the parameters.
\[r_t = \mu_{S_t} + \sigma_{S_t}\,\varepsilon_t, \qquad P(S_t = j \mid S_{t-1} = i) = p_{ij}, \qquad S_t \in \{0, 1\}.\]
Six parameters: two means, two variances, two staying probabilities. The filter returns \(P(S_t = \text{calm} \mid \text{data to } t)\) every day; the smoother uses the whole sample.
If the staying probabilities are \(p_{00} = 0.98\) (calm) and \(p_{11} = 0.96\) (stormy), the expected durations are…
Calm: sd 0.602 % a day, mean +0.107, expected duration 54.8 days. Stormy: sd 1.778 %, mean -0.096, duration 26.6 days — three times the volatility, a negative drift, and the desk spends 68 % of its days in the calm state. One fit, 0.3 seconds.
34 of 120 months are mostly stormy: January 2016, February and the fourth quarter of 2018, March–June 2020, all of 2022 — and August 2024. The model has no calendar; it found the bear markets from the variance alone.
Hold the index only when yesterday’s filtered probability of calm exceeds 0.5 (no look-ahead — the smoother is for the picture, the filter is for trading).
In the market 67.9 % of days. Regime filter: cumulative 76.4 %, Sharpe 0.81, max drawdown -11.3 %. Buy-and-hold: 105.0 %, 0.59, -41.4 %. It gave up a quarter of the return to cut the worst drawdown by three quarters. Like every volatility tool in this chapter it forecasts turbulence, not direction.
Next: §6.5 — the regime model said variance is what changes. Model it directly: ARCH and GARCH.
Returns look like white noise — until you square them. You will predict which ACF is flat and which is loud, turn the picture into a \(p\)-value, write a GARCH(1,1) estimator in 15 lines, forecast tomorrow’s variance and back-test a Value-at-Risk through March 2020.
For daily Dow Jones returns, which is more predictable from yesterday?
The daily standard deviation was 0.417 % in 2017 and 6.381 % in March 2020 — a 15× swing in scale, a 234× swing in variance. The worst day, −13.84 % on 2020-03-16, sits inside a cluster of other huge days. Quiet begets quiet; storms beget storms.
The white-noise band is \(\pm 2/\sqrt{1479} = \pm 0.052\).
Which pattern will acf(r) and acf(r**2) show at lags 1–5?
Squared returns: ACF 0.442, 0.551, 0.325, 0.324, 0.297 — every lag far outside ±0.052, and \(Q(10) = 1929\). Raw returns: −0.195, 0.157, −0.015, −0.080, 0.079 — small, though \(Q(10)=321\) is not zero either (the March 2020 reversals).
\[\varepsilon_t = \sigma_t z_t,\ z_t \sim N(0,1),\qquad \text{ARCH}(q):\ \sigma_t^2 = \omega + \sum_{i=1}^{q}\alpha_i \varepsilon_{t-i}^2,\qquad \text{GARCH}(1,1):\ \sigma_t^2 = \omega + \alpha\,\varepsilon_{t-1}^2 + \beta\,\sigma_{t-1}^2 .\]
\(\varepsilon_t\) stays white noise (\(\mathbb E[\varepsilon_t\mid\text{past}] = 0\)) while \(\varepsilon_t^2\) follows an AR(\(q\)). Unconditional variance \(\omega/(1-\alpha-\beta)\); persistence \(\alpha+\beta\).
Compared with ARCH(\(q\)), the \(\beta\) term makes the ACF of \(\varepsilon_t^2\)…
Engle’s ARCH-LM test for order \(q\) regresses:
\(R^2 = 0.362\) on a regression that “should” have \(R^2 = 0\): \(LM = 533.2\) against a \(\chi^2_5\) whose 1 % critical value is 15.1. het_arch reports the same 533.2 — you just re-implemented it. (On 1 479 i.i.d. normals the same code gives \(LM \approx 5\), \(p = 0.42\).)
Load appl.csv, build percent log returns for 2018–2020 into ra, and run arch_lm(ra.values, q=5). The starter leaves LM_a as a placeholder; replace the marked line so that LM_a is Apple’s ARCH-LM statistic.
Expected once fixed: \(LM = 132.6\), \(p < 10^{-4}\) on 756 days, with squared-return ACF 0.319, 0.302, 0.159. Apple clusters like the index.
In Colab
\(\hat\alpha = 0.1000\) and \(\hat\beta = 0.8800\) to four decimals. What is the most likely explanation?
arch rounds to two decimals\(r_t = \mu + \varepsilon_t\). Gaussian log-likelihood of one observation: \(-\tfrac12\left[\ln(2\pi\sigma_t^2) + \varepsilon_t^2/\sigma_t^2\right]\). Sum, negate, minimise — in percent, so \(\omega\) is of order 0.01, not \(10^{-6}\).
Equity indices typically show \(\alpha \approx 0.05\)–\(0.2\) and \(\alpha+\beta \approx 0.95\)–\(0.99\). Will alpha + beta land above or below 0.95?
0.975 — above
\(\hat\mu=0.082\), \(\hat\omega=0.037\), \(\hat\alpha=0.214\), \(\hat\beta=0.761\). Persistence 0.975 — a shock to variance halves in 27 days. Long-run sd 1.217 % vs sample 1.240 %. Well under a second.
From 0.43 % on 2017-11-09 to 8.88 % on 2020-03-17 — 141 % annualised. Volatility is not a number; it is a time series, and the model tracks it one day at a time.
Tomorrow is exact: \(\hat\sigma^2_{T+1} = \omega + \alpha\varepsilon_T^2 + \beta\sigma_T^2\). Beyond that, \(\mathbb{E}[\varepsilon^2_{T+h}] = \sigma^2_{T+h}\), so \(\hat\sigma^2_{T+h} = \omega + (\alpha+\beta)\,\hat\sigma^2_{T+h-1}\).
On 2020-11-13, \(\varepsilon_T = 1.283\) and \(\sigma_T^2 = 1.913\). What is \(\hat\sigma^2_{T+1}\) to 3 decimals?
1.845
1.845 → 1.810 → 1.704 → 1.560, sliding toward the long-run 1.48. A GARCH forecast is a decay path, not a level: the further out, the less today matters.
\(\text{VaR}_p = -\left(\hat\mu + \hat\sigma_{T+1}\, q_p\right)\), a positive loss. With \(q_{0.01} = -2.326\), \(\hat\mu = 0.082\), \(\hat\sigma_{T+1} = \sqrt{1.845} = 1.358\):
Standardised residuals still have excess kurtosis 2.43; the fitted \(\nu = 12.7\). Tomorrow’s 1 % VaR: 3.08 % (normal) vs 3.23 % (\(t\)). In-sample the normal 1 % VaR is breached 2.5 % of days, the \(t\) 2.2 % — both too often. At 5 % both are close (5.8 %, 6.1 %). Fat tails survive GARCH; they are just thinner.
Fit on 2015–2018 (1 006 days), freeze the parameters, run the recursion forward through 2019–2020 (473 days), count breaches of the 1 % VaR. Compare with an unconditional VaR from the training mean and sd.
Over 2019–2020 the 1 % VaR should be breached about 5 times (473 × 0.01). Which is closer to 5?
1 % VaR: GARCH 18 breaches (3.8 %, Kupiec \(p<0.001\)), unconditional 31 (6.6 %). 5 % VaR: GARCH 32 (6.8 %, \(p = 0.09\) — not rejected), unconditional 47 (9.9 %). GARCH-normal is still too thin at 1 %; the fixed-sd model is simply wrong. The unconditional line is crossed ten times in March 2020 alone; GARCH’s is crossed once — it was already at \(\sigma \approx 6\) %.
minimize with bounds. Dow 2015–2020: \(\alpha=0.214\), \(\beta=0.761\), half-life 27 days. Work in percent or the optimiser stalls at its starting values.Next: §6.6 — does bad news raise volatility more than good news? And what happens when you size positions by \(\hat\sigma_t\)?
One extra parameter lets the model react differently to losses and gains. You will test whether it is needed, stack GARCH on an ARMA mean, and then do what every systematic fund does with \(\hat\sigma_t\): not forecast direction, but scale the position — and see what that does to a drawdown.
Symmetric GARCH cannot tell them apart: \(\varepsilon_{t-1}^2\) is the same. Glosten, Jagannathan and Runkle (1993) add an indicator:
\[\sigma_t^2 = \omega + \left(\alpha + \gamma\,\mathbf 1[\varepsilon_{t-1}<0]\right)\varepsilon_{t-1}^2 + \beta\sigma_{t-1}^2 .\]
For a stock index, the sign of \(\gamma\) is expected to be:
Same likelihood, one more parameter. Setting the bound on \(\gamma\) to \((0,0)\) gives plain GARCH(1,1) — so one function fits both and the likelihood-ratio test is free.
\(LR = 2(\ell_{GJR} - \ell_{GARCH})\), \(\chi^2_1\) under \(H_0: \gamma=0\) (5 % critical value 3.84). Log-likelihoods −1815.05 and −1836.45. What is \(LR\)?
42.8
\(\hat\gamma = 0.254\), \(LR = 42.8\), \(p = 6\times10^{-11}\). A −2 % day feeds 0.335 of its square into tomorrow’s variance; a +2 % day only 0.081 — a 4.2× asymmetry. Persistence is unchanged at 0.974. The S&P 500 gives the same story (\(\hat\gamma = 0.272\), \(LR = 38.9\)).
Plot tomorrow’s variance as a function of today’s shock, holding \(\sigma_{t-1}^2\) at its long-run level.
After −3 %: variance 4.11; after +3 %: 1.82. The symmetric model splits the difference and is wrong on both sides — too calm after crashes, too nervous after rallies.
Everything so far assumed \(r_t = \mu + \varepsilon_t\). If the mean has ARMA structure, fit it first and hand the residuals to GARCH.
Two-step ARMA+GARCH means:
\(\hat\phi = -0.195\) (se 0.008) — but look at the standardised residuals: ACF of \(z\) at lag 1 is +0.151. The AR(1) was pulled by March 2020’s giant reversals; once each day is weighted by \(1/\sigma_t\) the mean dynamics look different. ACF of \(z^2\) (0.03, −0.01, 0.02) is clean — the variance is done. That mismatch is why joint estimation (armagarch, weighted by \(\sigma_t\)) is preferred when the mean matters.
Apple 2020–2024, 1 258 days. Fit GARCH(1,1), sort days into quartiles of the forecast \(\hat\sigma_t\) (known at the close of \(t-1\)), and compare next-day returns.
Compared with the calmest quartile, the most volatile quartile will show:
Calm quartile: mean 0.104 %/day, sd 1.35; stormy: mean 0.049, sd 2.91. Return per unit of daily risk falls from 0.077 to 0.017. Apple’s \(\alpha=0.086\), \(\beta=0.885\) — more persistent than the Dow.
Position \(=1\) when \(\hat\sigma_t <\) thr, else cash. The starter has thr = 100 — always invested, i.e. buy-and-hold. Set a real threshold (try 2.0; the median \(\hat\sigma\) is 1.70).
Change thr so the rule is in the market on fewer than 99 % of days, then read off cumulative return, annualised Sharpe and maximum drawdown versus buy-and-hold.
At thr = 2.0: in the market 73 % of days; cumulative 95 % vs 123 % for buy-and-hold, Sharpe 0.89 vs 0.78, max drawdown −25.5 % vs −37.7 %. Less money, better risk-adjusted money. An on/off switch is crude, though — the industry version is continuous.
Instead of in-or-out, hold \(w_t = \min\!\left(\dfrac{\sigma^\ast}{\hat\sigma_t},\, w_{\max}\right)\) of the index: more when it is calm, less when it is stormy, so that realised volatility stays near \(\sigma^\ast\). This is how CTAs, risk-parity funds and most systematic books scale every position.
Target \(\sigma^\ast = 1\) % a day, cap \(w_{\max} = 1.5\). Buy-and-hold over these five years realised 1.35 % a day with a −41 % drawdown.
Compared with buy-and-hold, the vol-targeted book should show…
Realised sd 0.994 % against the 1 % target (buy-and-hold 1.350). Vol-target: cumulative 60.7 %, Sharpe 0.77, max drawdown -22.6 %; buy-and-hold 59.9 %, 0.56, -41.4 %. Same money, half the drawdown. The weight ran from 0.13 in March 2020 to the 1.5 cap in quiet months. Nothing here forecasts direction — the whole edge is the variance forecast from §6.5.
Warning
From September 2011 the Swiss National Bank held EUR/CHF above 1.20, and for three years the exchange rate behaved like a textbook stationary series pinned to its floor: ADF rejected, KPSS did not, the spread to 1.20 mean-reverted within days. Traders and retail brokers treated the floor as a cointegrating rope between the franc and the euro, ran carry and short-volatility positions against it, and sized them on a variance estimated from those three quiet years.
On 15 January 2015 the SNB abandoned the floor without warning. EUR/CHF fell from 1.20 to 0.85 within minutes — a move of roughly 30 %, about 60 standard deviations of the pre-announcement daily change. FXCM, then the largest US retail FX broker, faced US$225 million of client losses in excess of collateral and needed a US$300 million rescue loan the next day; Alpari UK entered insolvency; Everest Capital’s US$830 million Global fund was wiped out.
Lesson for this chapter: a stationarity test tells you about the sample you fed it. A regime held in place by a policy is stationary until the policy changes, and no ADF statistic, GARCH \(\hat\sigma_t\) or smoothed regime probability can see the meeting at which it does. Test the spread, then ask who is holding the rope.
To: Head of Systematic Trading From: <Your name>, quant research Subject: Where to deploy the time-series budget: volatility-targeting overlay vs. a pairs book Date: 2026-03-02
Recommendation: Deploy the volatility-targeting overlay on the index book now; run the pairs book as a capped pilot (≤ 5 % of risk) until it survives two more formation windows.
Evidence: - Overlay, S&P 500 2020–24, target 1 %/day, cap 1.5×: Sharpe 0.77 vs 0.56, max drawdown −22.6 % vs −41.4 %, cumulative return unchanged (60.7 % vs 59.9 %). Regime filter alternative: Sharpe 0.81, drawdown −11.3 %, but gives up a quarter of the return. - Pairs, GOOG/S&P 2016–20, walk-forward hedge ratio, entry 2 / exit 0.5: Sharpe 0.75, drawdown −11.6 %, in a trade 35 % of days, 52 crossings before costs. Full-sample
cointp = 0.014 — but p = 0.24 and 0.11 on the formation windows a trader would actually have had.Caveats: the overlay depends only on variance persistence (\(\alpha+\beta \approx 0.97\)), which has held in every equity market since 1987; the pairs book depends on a rope whose in-sample evidence is stronger than its out-of-sample evidence, and its Sharpe doubles or halves with one threshold. Neither is net of costs or slippage.
Next step: overlay — add GJR (bad-news \(\gamma\) = 0.25) and refit weekly. Pairs — extend to 20 candidate pairs, require formation-window p < 0.05 and a half-life under 30 days, and stop trading any pair whose rolling ADF p exceeds 0.10.
arch_model(r) on decimal returns — the stalled fit in the notebook (alpha=0.1000, beta=0.8800). “Scale returns to percent first, and tell me the optimiser’s convergence flag.”coint p = 0.01 and then back-tests a z-score rule on the same window has used the future twice. Ask it to state the formation window, the trading window, and the hedge ratio it used on each day — and to refuse a Sharpe on any spread whose formation-window test did not reject.Pitfall: when a copilot reports a regime model or a Kalman beta, ask whether it used the smoothed (two-sided) or filtered (one-sided) estimate in the back-test. The smoothed one knows the future.
| Role | Concept | Tool |
|---|---|---|
| Foundation | Stationarity, unit roots, decision table | kpss (null: stationary), adfuller (null: unit root), "c" / "ct" |
| Foundation | Memory diagnostics, ARIMA baseline | plot_acf, plot_pacf, ±\(2/\sqrt n\); ARIMA(y, order).fit().aic, .get_forecast, .apply |
| ML base | Forecast scoring | MAE / RMSE / DA on the same rows; the zero forecast |
| ML base | Features without leakage, walk-forward | .shift(1), .rolling(w), expanding yearly refit of HistGradientBoostingRegressor |
| Structural arb | Cointegration and pairs | coint(a, b); formation vs trading window; rolling hedge ratio; z-score entry/exit |
| Structural arb | Hidden states | UnobservedComponents(level="local level"); 8-line time-varying beta; MarkovRegression(k_regimes=2, switching_variance=True) |
| Volatility | ARCH/GARCH, forecasts, VaR | het_arch; numpy recursion + optimize.minimize; \(\hat\sigma^2_{T+h}\) decay path; stats.norm.ppf / stats.t.ppf; Kupiec |
| Volatility | Leverage, vol targeting | GJR indicator + LR test; ARIMA(...).resid → GARCH; \(w_t = \min(\sigma^\ast/\hat\sigma_t, w_{\max})\) |
The message, once more: the classical models lost the mean-forecasting job to ML — and kept the variance, the spread, the regime, and the pipeline.
The walk-forward GBM on the S&P 500 scored DA 0.494 and a Sharpe below buy-and-hold, yet Chapter 4.5’s cross-sectional model earned a positive spread. What, precisely, is different about the cross-section that makes the same algorithm work? Name two things a single series cannot provide.
GOOG/S&P 500 passed the Engle–Granger test on the full sample (p = 0.014) but not on either formation window (p = 0.24, 0.11). A colleague proposes trading it anyway because “the rope was there all along”. Write the two-sentence reply, with numbers, and say what evidence would change your mind.
The regime filter and the vol-targeting overlay both cut the S&P 500 drawdown by more than half, using only a variance model. Under what return-generating process would either of them raise the total return, not just the Sharpe ratio — and does the quartile table on Apple support that process?
The Kalman beta tracked NVDA’s 2024 swing (2.7 → 3.2 → 1.7) a year before annual OLS could, yet the hedge error was the same to two decimals. When does timeliness in a hedge ratio pay, and when is it just noise? Propose one experiment on the NVDA/SPY data that would tell the two apart.
Prof. Xuhu Wan · HKUST ISOM · Learning Statistics with Python