Functional Volatility Surface Modelling
September 2024 – present- Extending the functional GARCH framework to a generalized autoregressive score (GAS) model to estimate and capture time-varying intraday volatility surfaces.
- Designed efficient estimation procedures using B-splines, applying
numbaJIT compilation to enable scalable modelling of volatility surfaces from granular intraday return data.
Overview
Treats continuous intraday log-return paths as functional data objects \(y_t(u)\) over the trading day \(u \in [0,1]\). The Functional GARCH(p,q) model generalizes classical volatility dynamics to an infinite-dimensional Hilbert space, capturing how shocks at any point of the day impact the entire upcoming volatility surface.
Model Definition
The functional scale model decomposes returns using a time-varying conditional variance curve \(\sigma_t^2(u)\), so we look at \[y_t(u) = \sigma_t(u)\eta_t(u)\] driven by the recursion \[\sigma_t^2 = \delta + \sum_{i=1}^{q}\alpha_i\left(y_{t-i}^2\right) + \sum_{j=1}^{p}\beta_j\left(\sigma_{t-j}^2\right)\]
where \(\delta\) is a strictly positive baseline intercept curve, and \(\alpha_i, \beta_j\) are non-negative integral kernel operators mapping past squared return curves and past volatility curves into today's volatility surface.
Functional GAS Extension
The functional scale model can be extended to a functional Generalized Autoregressive Score (GAS) model by replacing the autoregressive dependence on past squared return curves with a score-driven update. The returns are defined as \[y_t(u) = \sigma_t(u)\eta_t(u),\] where \(\sigma_t^2(u)\) denotes the time-varying conditional variance curve. The functional GAS dynamics are specified as \[\sigma_{t+1}^2 = \delta + B\sigma_t^2 + A\int \phi_K(s)y_t^2(s)\,ds,\] where \(\delta\) is a strictly positive baseline intercept curve, \(B\) is a linear persistence operator acting on the previous volatility curve, and \(A\) maps the observed squared return curve into the updated volatility curve. The basis projection term \[\int \phi_K(s)y_t^2(s)\,ds\] represents the score-driven innovation component obtained from the functional representation of the return process.
Equivalently, the recursion can be expressed using integral kernel operators as \[\sigma_{t+1}^2(u) = \delta(u) + \int B(u,v)\sigma_t^2(v)\,dv + \int A(u,s)\phi_K(s)y_t^2(s)\,ds.\] The functional GAS model therefore updates the conditional variance curve using information contained in the most recent squared return curve, while the operator \(B\) captures the persistence of past volatility curves. Compared with the functional scale model, \[\sigma_t^2 = \delta + \sum_{i=1}^{q}\alpha_i(y_{t-i}^{2}) + \sum_{j=1}^{p}\beta_j(\sigma_{t-j}^{2}),\] the GAS formulation replaces the ARCH-type feedback operators \(\alpha_i(y_{t-i}^{2})\) with a score-driven update based on the functional projection of the return innovations.
Estimation
Because standard likelihood functions cannot be directly evaluated for continuous curves, the model is estimated using Functional Quasi-Maximum Likelihood Estimation (QMLE).
Setup
git clone https://github.com/DaanZunnenberg/FunctionalScale.git
cd FunctionalScale
pip install -e . # editable install, pulls in numpy/scipy/numba/pandas/matplotlib/tqdm
pip install -e ".[dev]" # + pytest, jupyter (optional, for tests/notebooks)
Usage: functional GARCH
import numpy as np
from funcgarch import fit, garch_filter
# mY: (N, T) matrix of intraday return curves, N grid points per day, T days
mY = np.load("returns.npy")
N, T = mY.shape
M = 4 # number of Bernstein basis functions
result = fit(mY, n_grid=N, M=M) # QMLE-style estimation (scipy.optimize)
vtheta_hat = result.x
sigma2 = garch_filter(mY, n_grid=N, vtheta=vtheta_hat, M=M) # (N, T) fitted variance surface
Usage: functional GAS-GARCH
from scipy.optimize import minimize
from funcgarch import gas_garch_estimator
# vtheta = [nu, ou_scale, omega (M,), vec(B) (M*M,), vec(A) (M*M,)]
result = minimize(
gas_garch_estimator, x0=vtheta_init, args=(mY, N, M),
method="SLSQP",
)
Because the score-driven update adapts its B-spline coefficients every day rather than fitting one static operator to the whole sample, the GAS-GARCH fit tracks the true surface considerably more tightly than the plain functional GARCH fit above:
The first comparison figure shows the estimated GAS-GARCH volatility surface against the true volatility surface to assess the model’s ability to recover the underlying dynamics.
Placing the two estimators' fitted surfaces directly next to each other, rather than each against the true surface separately, makes the score-driven adaptation's smoothing effect easier to see:
The second comparison figure compares the functional GARCH and GAS-GARCH volatility surfaces, demonstrating the increased flexibility of the GAS-GARCH specification relative to the traditional functional GARCH model.
Data Flow
wrds/*.sas scripts/taq_cleaner.py funcgarch/*.py
┌─────────────────┐ ┌────────────────────┐ ┌──────────────────────┐
│ WRDS TAQ pull │ raw CSV │ DataCleaner.clean()│ mY (N,T) │ fit() / garch_filter │
│ (data_fetcher, │ ────────► │ - align to grid │ ────────► │ gas_garch_estimator │
│ taq_cleaner, │ │ - compute returns │ │ func_garch_estimator │
│ nbbo/dynamic_ │ │ - reshape to │ │ │
│ taq_minute, │ │ (N, T) matrix │ │ -> vtheta_hat, │
│ export) │ │ │ │ sigma2 surface │
└─────────────────┘ └────────────────────┘ └──────────────────────┘
Requires Python ≥ 3.9. Full theory, references, and repository layout in the README.