Univariate Time Series
Single-series models in two movements. First the conditional mean: Box–Jenkins SARIMA forecasting, and the change-point → Markov-switching ladder — a Bayesian AR with mean and variance shifts and a Markov-switching AR — that lets the dynamics break. Then the conditional variance, a full volatility arc: GARCH, asymmetric GJR-GARCH with Value-at-Risk backtesting, GARCH-in-mean (the risk–return tradeoff), stochastic volatility, realized volatility (HAR and Realized GARCH from intraday data), and regimes vs long memory in volatility persistence — capped by an application that prices options under GARCH and stochastic volatility. Each is built from scratch and cross-checked across Python and R.
Seasonal ARIMA (SARIMA)
Seasonal time-series modelling and forecasting — the Box–Jenkins SARIMA
framework, comparing the frequentist workflow (statsmodels / pmdarima and R's
forecast::auto.arima, with automatic Hyndman–Khandakar order selection and the
Kalman-filter likelihood) against a Bayesian structural alternative — a linear
trend + monthly seasonal effects + AR(1) errors fit by PyMC/NUTS, giving full
posterior-predictive forecasts. The two encode seasonality differently — SARIMA by seasonal
differencing (stochastic, drifting), the structural model by fixed monthly dummies
(deterministic) — both valid when the seasonal shape is stable. Five examples span the range and
build the ARIMA → changepoint → regime-switching ladder:
AirPassengers (multiplicative seasonality, ~8–9% MAPE), sunspots
(cyclicality — a variable ~11-year period from AR complex roots, a damping forecast, and a
Bayesian posterior over the cycle length), Mauna Loa CO₂ (additive trend +
seasonality, near-perfect ~0.3% MAPE), US GNP growth (the weak, irregular,
barely-forecastable business cycle — the linear benchmark for Markov-switching), and the
Nile (a structural break / intervention, with a Bayesian changepoint that
estimates the 1899 break year itself). Cross-checked across statsmodels, PyMC, and R.
Bayesian AR with Mean and Variance Shifts
A Python implementation of the McCulloch & Tsay (1993) 4-block Gibbs sampler for an AR(p) model with piecewise-constant mean and innovation variance. The key reparametrisation — raw AR residual rt ≈ κμt + et — reduces within-segment inference to conjugate Normal–Normal (mean) and Inverse-Gamma–Normal (variance) problems. Single-site sweeps on the shift indicators compare marginal likelihoods of a merge vs. a split for every active time point. Applied to: Section 1 a synthetic AR(2) with a planted mean shift (detected at P = 0.88) and variance shift; Section 2 monthly US retail gasoline price changes (1978–1991) — reproduces M&T's AR(2) coefficients ≈ (0.91, −0.30), no mean shift, and all five variance-shift episodes (Iranian Revolution, 1981 OPEC peak, 1982, 1986 oil collapse, 1989); Section 3 the full 1976–2026 BLS series, recovering variance regimes through the 2022 Russia–Ukraine invasion and the live 2026 Iran/Hormuz crisis; Section 3b a log-return refit showing the ~2000 volatility increase is real but 3–4× smaller than the nominal cent analysis suggests.
View example →Bayesian Markov-Switching AR
A Python implementation of the Hamilton (1989) K-regime switching-intercept AR(p) model via two Bayesian routes: a from-scratch Forward-Filter Backward-Sampling (FFBS) Gibbs sampler (pure NumPy) and a NumPyro NUTS fit that marginalizes the discrete states with the Hamilton forward filter. FFBS draws the entire state path S1:T jointly each sweep — exact joint draw via the Markov factorization — and recovers smoothed regime probabilities by averaging sampled paths. Regimes are identified by relabelling each sweep so intercepts are ascending. Applied to: Section 1 synthetic AR(2) with planted regimes (88.7% classification); Section 2 Hamilton's quarterly US real GNP (1951Q2–1984Q4) — intercepts −0.27 / 1.05 matching Hamilton's −0.36 / 1.16, NBER recession classification 95% (FFBS) / 91% (NumPyro); frequentist cross-validation on US gasoline log-returns with Python statsmodels and R MSwM — three independent implementations agree on AR(2) dynamics and the ~2000 volatility transition to within rounding, completing a three-way validation with the McCulloch–Tsay change-point model.
View example →Bayesian GARCH
Volatility modelling for financial returns — where SARIMA models the mean, GARCH
models the moving variance (volatility clustering). The same GARCH(1,1) is fit on daily
S&P 500 returns (1987–2009) across a full trio, and — in the from-scratch notebook —
five ways: the classical MLE plus four Bayesian methods that predate today's
HMC/NUTS samplers (Gauss–Hermite quadrature, importance sampling,
Metropolis–Hastings, and griddy Gibbs), which all recover one
posterior (α+β ≈ 0.993). The dimension-general module also spans GJR leverage and Normal/Student-t
innovations; adding fat tails gives ν ≈ 6.3 and a decisive Normal-vs-t Bayes factor
(log BF ≈ +199) read straight off the Gauss–Hermite marginal likelihood. Cross-checked against
PyMC/NUTS and R (rugarch ML + bayesGARCH MCMC) — four engines, one answer.
Bayesian Asymmetric GARCH (GJR) & Value-at-Risk
The asymmetric extension of GARCH: GJR-GARCH
adds a leverage term so bad news raises volatility more than good news.
The same five from-scratch estimators fit the extended model on S&P 500 returns; leverage is
large and unanimous (γ ≈ 0.13 — a negative shock moves next-day volatility ~15× more than a
positive one), drawn directly by the Engle–Ng news-impact curve. Because Gauss–Hermite returns
each marginal likelihood, the four models rank by Bayes factor — GARCH-N → +t →
+leverage → GJR-t (best). Cross-checked against PyMC/NUTS and R
(rugarch gjrGARCH/eGARCH). A fourth notebook applies the volatility models to
Value-at-Risk and backtests them through the 2008 crisis (Kupiec &
Christoffersen): Historical Simulation fails the independence test, RiskMetrics' Gaussian tail
under-covers, and only GARCH-t / GJR-t pass cleanly at the 1% level.
Bayesian GARCH-in-Mean — the Risk-Return Tradeoff
Does higher volatility earn a higher expected return? GARCH-in-Mean (Engle,
Lilien & Robins 1987) puts the conditional volatility into the mean equation, so the
new coefficient λ is the price of risk. Because the mean feeds
back into the variance recursion, the likelihood is a genuine forward loop (no vectorised filter).
On daily S&P 500 returns the premium has the right sign but is statistically weak — the classic
risk-return "puzzle": λ > 0 everywhere yet no interval excludes zero. Two
refinements sharpen it — fat tails absorb much of the apparent premium (it halves
under Student-t), and lower frequency helps (monthly λ is ~10× larger, a
signal-to-noise effect). Estimated four ways across three engines — from-scratch MLE + Metropolis,
PyMC/NUTS, and R rugarch (archm=TRUE) — all agreeing, extending the
GARCH / GJR family.
Bayesian Stochastic Volatility
The parameter-driven alternative to GARCH: in a stochastic-volatility model the
log-variance is its own latent AR(1) with an independent shock, so — unlike a
GARCH filter — the fitted volatility path carries a genuine posterior credible band.
The from-scratch notebook builds the Kim–Shephard–Chib (1998) sampler by hand,
exposing that SV is a state-space model: a 7-component mixture linearises the log-squared
returns, and the whole log-variance path is drawn in one block by the Kalman filter + backward
sampling (FFBS). Three engines agree on the S&P 500 posterior (φ ≈ 0.987,
ση ≈ 0.15) — from-scratch, NumPyro/NUTS, and R's stochvol (Kastner ASIS).
Extensions add SV-t (fat tails, ν ≈ 9) and leverage (ρ ≈ −0.6, the
SV analogue of GJR γ > 0), ranked by WAIC. Together with
GARCH and GJR, the
arc recovers the same three stylized facts — persistence, fat tails, leverage — written three ways.
Realized Volatility — HAR & Realized GARCH
Where GARCH and SV infer volatility from daily returns, realized volatility measures it from intraday prices (summing squared high-frequency returns). A three-phase
arc: (1) the realized measures and their two stylized facts — long memory and
log-normality; (2) Corsi's HAR-RV, a 4-parameter regression of log-RV on its
daily/weekly/monthly averages that beats a fitted GARCH out-of-sample — built three
ways (from-scratch / PyMC / R highfrequency), with HAR-CJ (jumps) and HARQ
(measurement error) extensions; and (3) Realized GARCH (Hansen–Huang–Shek 2012),
the formal bridge, plus a forecasting tournament. The verdict: realized-measure
models beat returns-only ones (GARCH,
SV) by ~25% on QLIKE — using intraday data matters more than the
model class.
Volatility Persistence: Regimes vs Long Memory
Every GARCH fit reports the same near-unit-root persistence (α+β ≈ 0.997) — but why?
Two rival explanations, which Diebold & Inoue (2001) proved are observationally confounded:
long memory (FIGARCH, fractional integration — the persistence
is real, decaying hyperbolically) vs regime switching (MS-GARCH —
the persistence is a spurious artifact of shifts in the volatility level). Both are fit from
scratch (FIGARCH by fractional differencing; MS-GARCH via the Haas–Mittnik–Paolella spec, Hamilton
filter + Kim smoother), each with an R package cross-check (rugarch,
MSGARCH). The verdict is both, but the persistence is long memory:
d = 0.52 (genuine long memory, matching the realized-volatility
evidence), and allowing regimes does not dissolve the within-regime persistence — the
regimes add a business-cycle layer on top, but the near-unit-root persistence is real.
Option Pricing under GARCH & Stochastic Volatility
The application phase of the volatility arc — the same models that fit equity returns now price options, replicating three seminal papers. Heston–Nandi (2000): a GARCH with an exponential-affine generating function, so options price in closed form by Fourier inversion (martingale identity to machine precision). Duan (1995): Monte-Carlo GARCH pricing via the locally risk-neutral valuation relationship — validates the closed form, then prices a path-dependent Asian option no closed form can. Heston (1993): the continuous-time stochastic-vol closed form via the characteristic function, reducing to Black–Scholes as vol-of-vol → 0. The unifying thread is the leverage parameter: the γ that GJR-GARCH found in returns and the ρ of stochastic volatility are the same asymmetry that tilts a symmetric option smile into the equity-index skew.
View example →