Research code for functional GARCH and functional GAS-GARCH models applied to intraday volatility surfaces. The models treat within-day return curves as functional observations, extending classical GARCH/GAS dynamics to function space (
Let
The conditional variance curve on day
Following Cerovecki et al. (2018), the model specifies:
where multiplication is pointwise, the intercept curve
$$\alpha(x) = \int_0^1 K_\alpha(u, v), x(v), dv, \qquad K_\alpha(u,v) \geq 0$$
A unique strictly stationary, non-anticipative solution exists if the top Lyapunov exponent of the companion operator sequence satisfies
where
To estimate the infinite-dimensional operators, we expand them using a basis of non-negative Bernstein polynomials
where
The parameters are estimated via Quasi-Maximum Likelihood:
Note: The current implementation (
loss_funcingarch.py) uses a Bernstein-projected MSE approximation rather than the formal QMLE criterion. The inner products are evaluated via sample means on a uniform grid of$N$ points, optimized usingscipy.optimize(SLSQP).
The Generalized Autoregressive Score (GAS) extension maps the infinite-dimensional update to a score-driven recursion on a low-dimensional cubic B-spline coefficient vector
Let
The coefficients evolve via the GAS recursion:
Evaluated on a discrete grid
Let
where
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)Requires Python >= 3.9.
importnumpyasnpfromfuncgarchimportfit, garch_filter# mY: (N, T) matrix of intraday return curves — N grid points per day, T daysmY=np.load("returns.npy")
N, T=mY.shapeM=4# number of Bernstein basis functionsresult=fit(mY, n_grid=N, M=M) # QMLE-style estimation (scipy.optimize)vtheta_hat=result.xsigma2=garch_filter(mY, n_grid=N, vtheta=vtheta_hat, M=M) # (N, T) fitted variance surfacefromscipy.optimizeimportminimizefromfuncgarchimportgas_estimator# vtheta = [nu, ou_scale, omega (M,), vec(B) (M*M,), vec(A) (M*M,)]result=minimize(
gas_estimator, x0=vtheta_init, args=(mY, N, M),
method="SLSQP",
)# 1. Pull raw TAQ data via SAS (see wrds/data_fetcher.sas)# 2. Clean into an intraday return matrix
python -m scripts.taq_cleaner --parent-path /path/to/taq/exportsfromdata.taq_cleanerimportDataCleanercleaner=DataCleaner(ticker="SPY")
mY=cleaner.clean(parent_path="/path/to/taq/exports") # -> (N, T) return matrixwrds/*.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_estimator │
│ taq_cleaner, │ │ - compute returns │ │ bspline_garch_ │
│ nbbo/dynamic_ │ │ - reshape to │ │ estimator │
│ taq_minute, │ │ (N, T) matrix │ │ -> vtheta_hat, │
│ export) │ │ │ │ sigma2 surface │
└────────────────┘ └───────────────────┘ └─────────────────────┘
Dependency direction: funcgarch (core package) has no dependency on data/ or wrds/ — it only
consumes an (N, T) return matrix. The data/taq_cleaner.py pipeline depends on raw output from the
wrds/*.sas scripts. examples/*.ipynb depend on both funcgarch and pre-cleaned data.