3 GARCH volatility forecasting in R
3.1 Volatility as a forecasting problem
An investor observing a quiet market today may expect relatively small price movements tomorrow. After several turbulent sessions, the same investor may revise that assessment even without changing the expected direction of the market. This chapter develops a model for that changing assessment of risk. Where the previous chapters forecast the level of a monthly sales series, we now ask how much a daily financial return is likely to vary, given what has already happened in the market. The implementation builds on the GARCH example in Lozano (2025) and the risk-management applications discussed by Hull (2022) and McNeil et al. (2015).
The starting observation is that large financial returns tend to occur close together in time. A turbulent session is often followed by further large movements, while quiet periods can persist for many days. A standard deviation calculated from the entire sample averages over both kinds of period. To describe the risk faced on a particular date, we instead allow variance to change as new returns are observed. A GARCH model provides a rule for making these updates and carrying them forward into a forecast.
Let \(u_t\) be the daily return and let \(\mathcal{F}_{t-1}\) denote the information available at the end of day \(t-1\). The object of interest is
\[ v_t = \operatorname{Var}(u_t \mid \mathcal{F}_{t-1}). \]
The conditioning information in this expression fixes the timing of the question: \(v_t\) is the variance assigned to day \(t\) using information available at the end of day \(t-1\). Its square root, \(\sigma_t=\sqrt{v_t}\), is the conditional volatility, expressed in the same units as the return. By constructing the model step by step, we can follow how the data lead to these daily risk estimates and eventually to a forecast:
\[ \text{returns} \longrightarrow \text{GARCH recursion} \longrightarrow \text{likelihood} \longrightarrow \text{estimated volatility} \longrightarrow \text{volatility forecast}. \]
The first arrow begins with a measurement choice. We convert index levels to returns so that a price movement is measured relative to the amount invested. Those returns become inputs to the GARCH recursion, which calculates today’s variance from yesterday’s squared return and variance. At first, however, the parameters of that rule are unknown. Different choices produce different variance paths, and we need a way to decide which path best accounts for the observations. The likelihood supplies that comparison by measuring how plausible the observed returns are under each proposed path. Estimation therefore moves repeatedly between the recursion and the likelihood until the numerical search finds the best-fitting parameters.
Once those parameters have been selected, the recursion gives a fitted variance for each date; taking square roots produces the estimated volatility series in the diagram. Extending this series into the future requires one further change. Historical updates can use returns that have already occurred, whereas future updates depend on shocks that have yet to arrive. Replacing their squared values by the expectations implied by the model gives the forecast. The distinction between observed shocks and expected future shocks will explain both the shape of that forecast and how it changes when new data become available. Before developing the updating rule, we need to construct the returns on which it operates.
3.2 Returns and data preparation
Our application uses the S&P 500 index, with prices requested from 2015-07-10 to a download cutoff of 2020-07-10. Because the index level changes over time, a movement of a given number of points can represent very different changes in the value of an investment. We account for this by dividing the change in the closing level, \(S_t-S_{t-1}\), by the preceding close, \(S_{t-1}\). The resulting simple daily return is
\[ u_t = \frac{S_t-S_{t-1}}{S_{t-1}} = \frac{S_t}{S_{t-1}}-1. \]
For example, a rise of 10 points from 100 to 110 gives a 10% return, while the same rise from 1,000 to 1,010 gives a 1% return. The absolute changes are identical, but the first movement is ten times as large relative to the starting value. Measuring returns therefore lets us compare risk across assets with different prices and across periods in which the same index trades at different levels. Throughout the calculations, we express these returns as decimals, so that 0.01 consistently represents 1%.
Working with returns also brings the statistical measurement closer to the financial question. Price levels often trend and retain the accumulated effects of past shocks; their dispersion over a long sample can consequently reflect years of growth as well as short-run uncertainty. Daily returns isolate the relative movements whose variability we want to forecast. This transformation still leaves an important feature to explain: the size of those movements can change markedly from one period to another. GARCH models that changing conditional variance. Whether the returns are stationary and whether this particular model captures their behavior remain questions to assess from the data.
Code
sp500 <- tq_get("^GSPC", from = "2015-07-10", to = "2020-07-10") |>
select(date, close) |>
rename(price = close) |>
mutate(
return = price / lag(price) - 1
)
sp500 |>
select(date, price, return) |>
slice(c(1:6, (n() - 1):n())) |>
kable(
caption = "S&P 500 prices and simple daily returns.",
col.names = c("Date", "S&P 500 level", "\\(u_t\\)"),
escape = FALSE,
digits = 8,
format.args = list(scientific = FALSE),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Date | S&P 500 level | \(u_t\) |
|---|---|---|
| 2015-07-10 | 2076.62 | NA |
| 2015-07-13 | 2099.60 | 0.01106605 |
| 2015-07-14 | 2108.95 | 0.00445316 |
| 2015-07-15 | 2107.40 | -0.00073499 |
| 2015-07-16 | 2124.29 | 0.00801468 |
| 2015-07-17 | 2126.64 | 0.00110618 |
| 2020-07-08 | 3169.94 | 0.00782746 |
| 2020-07-09 | 3152.05 | -0.00564361 |
The table contains prices from 2015-07-10 through 2020-07-09. Its first row has no return because calculating a change requires a preceding close. We remove that initial missing value when constructing the return vector, while retaining the original prices for the calculation tables below. The resulting return sample will serve throughout as our training data. Later, we will evaluate forecasts using subsequent observations that have played no part in estimating the parameters.
Plotting these returns reveals why a single measure of dispersion would be an incomplete description of risk. Around the 2020 market stress episode, large positive and negative movements appear close together, in contrast to the much smaller movements observed during quieter periods. An updating rule for variance should be able to respond to this clustering.
Code
3.3 The GARCH(1,1) model
To describe the clustering visible in the return plot, we let a large movement raise the variance assigned to the following day and allow some of that increase to carry forward. The GARCH(1,1) model of Bollerslev (1986) implements both ideas in one equation. It extends the ARCH model of Engle (1982) by including the preceding conditional variance alongside the preceding squared shock:
\[ v_t = \omega + \alpha u_{t-1}^2 + \beta v_{t-1}. \]
The intercept and the weights on these two sources of information are the three parameters we will estimate:
| Parameter | Interpretation |
|---|---|
| \(\omega\) | Baseline variance component. |
| \(\alpha\) | Reaction to the most recent squared return. |
| \(\beta\) | Persistence from the previous conditional variance. |
The use of a squared return follows from the question we are asking. A large gain and a large loss both indicate substantial variation, even though their directions differ. Under the zero-conditional-mean assumption used here, \(u_t\) is also the unexpected return, or shock, and \(u_t^2\) measures its magnitude without retaining its sign. Returns of 2% and -2%, for instance, both give \(u_t^2=0.0004\), four times the squared value of a 1% return. Through \(\alpha u_{t-1}^2\), a large movement in either direction can therefore raise the next variance estimate. A single squared return is a noisy indication of risk, however, so the equation also draws on the previous variance instead of relying entirely on that one observation. The equal response to positive and negative shocks of the same size is a feature of this symmetric model.
The preceding variance enters through \(\beta v_{t-1}\), allowing the model to retain some of its earlier assessment of risk. This provides one source of persistence, but the squared-return term adds another: when variance is high, larger squared returns are also expected, and these feed back into subsequent updates. Taking expectations over future shocks combines the two effects into \(\alpha+\beta\). If this sum is close to one, an increase in variance is expected to fade slowly. Along an observed path, further shocks can interrupt that decline, so persistence describes expected propagation rather than a smooth recovery that must occur after every turbulent day.
For this updating rule to produce meaningful variances and a finite long-run level, we impose the restrictions
\[ \omega>0, \qquad \alpha \ge 0, \qquad \beta \ge 0, \qquad \alpha+\beta<1. \]
The first three conditions protect the sign of the calculated variance. With a positive intercept and nonnegative weights, the terms in the recursion cannot combine to produce a negative value. Allowing \(\omega<0\) would remove that protection on a quiet day, when the other contributions are small; allowing \(\alpha<0\) would make a sufficiently large squared shock push the next variance below zero. A negative \(\beta\) would create a similar problem when the preceding variance is large. We permit either weight to equal zero, which simply removes its channel, but require the intercept to remain strictly positive. With \(\alpha+\beta<1\), setting the intercept to zero would imply a degenerate zero long-run variance.
Positivity alone does not determine whether the process settles around a finite average level. To examine that question, take unconditional expectations in the recursion. Under the zero-mean assumption, the average squared return equals the average conditional variance, denoted by \(\bar v\). The equation then becomes \(\bar v=\omega+(\alpha+\beta)\bar v\), which gives
\[ \bar{v} = \frac{\omega}{1-\alpha-\beta}. \]
This expression explains the restriction on the sum of the weights. For \(\omega>0\), a positive, finite long-run variance requires \(\alpha+\beta<1\). At a sum of one, expected variance fails to return to a finite level, and above one its expected dynamics are explosive. The finite-second-moment condition is the one needed for our variance-targeting estimator and mean-reverting forecasts below. Strict stationarity is a distinct, more general issue; treatments of that distinction and of financial time-series modeling more broadly can be found in Francq and Zakoian (2019), Tsay (2010), and Ruppert (2011).
3.4 Conditional variance recursion
The equation gives us a rule for one day, but applying it to the sample requires a sequence of calculations. Once we have obtained \(v_t\) from \(u_{t-1}^2\) and \(v_{t-1}\), that new variance becomes an input to the calculation of \(v_{t+1}\). We combine it with \(u_t^2\) and continue forward in the same way. This dependence of each calculation on the preceding result is what makes the model recursive.
To begin the sequence, we need a variance value that the data do not directly provide. We use the first available squared return to initialize the variance for the following return date. Counting from the first row of prices makes the timing clear: day 1 supplies the initial close, day 2 supplies the first return, and day 3 receives the initialized variance \(v_3=u_2^2\). The full GARCH update can then begin on day 4, using the return and initialized variance from day 3. Because a return and a positive variance are available together from day 3 onward, those are also the observations that enter the likelihood.
The R function implements this sequence using the following objects:
| Mathematical object | R object | Description |
|---|---|---|
| \(u_t\) | u[t] |
Daily return, also the shock under the zero-mean assumption. |
| \(v_t\) | v[t] |
Conditional variance assigned to that return date. |
| \(\omega\) | omega |
Positive intercept in the variance recursion. |
| \(\alpha\) | alpha |
Weight on the preceding squared return. |
| \(\beta\) | beta |
Weight on the preceding conditional variance. |
In this mapping, t indexes the return vector after the first missing return has been removed. Its first element corresponds to price-row day 2. The calculation tables below retain the original price-row day numbers so that the initial missing values remain visible.
The function first fills a variance vector with missing values and sets v[2] equal to u[1]^2. The loop can then begin at t = 3, using u[t - 1]^2 and v[t - 1] in exactly the order required by the equation. These indices refer to the shortened return vector, so this first loop iteration corresponds to day 4 in the original price data. We now have a function that produces a complete variance path for any admissible parameter vector. The remaining question is how to choose that vector from the data.
3.5 Maximum likelihood estimation
Suppose we try a parameter vector that assigns very low variance to a period of large returns. Those observations would be difficult to explain under the proposed model. Raising the variance would make them more plausible, but a rule that assigns high variance on every date would also fit quiet periods poorly. Maximum likelihood formalizes this comparison by evaluating how much probability density the model assigns near each observed return, conditional on the information preceding it. We use density because returns are continuous: the probability of any one exact value is zero, while the probability of a small interval around it depends on the density there.
To make this comparison numerical, we need a distribution whose spread is determined by the GARCH variance. We assume
\[ u_t \mid \mathcal{F}_{t-1} \sim \mathcal{N}(0,v_t), \]
so that a small \(v_t\) concentrates the normal density around zero and a larger \(v_t\) spreads it over a wider range of possible returns. The zero mean lets us focus on this changing spread; in an application with a separate mean model, such as a constant or an ARMA equation, its residuals would take the place of \(u_t\) in the variance calculation.
Multiplying the conditional densities over the sample gives the likelihood of the observed sequence, conditional on our initialization. We work with its logarithm because it turns this product into a sum without changing which parameter vector is preferred. This also avoids the numerical difficulty of multiplying many small densities. For estimation, it is convenient to begin with the following proportional daily contribution:
\[ \ell_t^\ast = -\log(v_t) - \frac{u_t^2}{v_t}. \]
The trade-off described above is visible in the two terms. If a large return occurs when \(v_t\) is small, the ratio \(u_t^2/v_t\) becomes large and reduces the contribution through its negative sign. Increasing \(v_t\) eases that penalty, but also lowers \(-\log(v_t)\), so making every variance arbitrarily large cannot solve the fitting problem. The optimizer must balance these effects across the sample using one common set of GARCH parameters.
Adding the daily contributions produces \(\ell^*=\sum_t\ell_t^\ast\). The asterisk reminds us that this score has been simplified for estimation. To report the complete Gaussian log-likelihood, we restore the normal-density constant and the factor of one half:
\[ \ell = -\frac{1}{2} \sum_t \left[ \log(2\pi) + \log(v_t) + \frac{u_t^2}{v_t} \right]. \]
Both expressions use the same variance path and the same dates. Their dependence on \((\omega,\alpha,\beta)\) comes through the variances generated by the recursion; the additional constant contains no unknown parameter. Consequently, if \(n\) is the number of included observations, we can obtain the full Gaussian total directly from the proportional total:
\[ \ell = \frac{\ell^\ast}{2} - \frac{n}{2}\log(2\pi). \]
Dividing every candidate score by two and subtracting the same constant preserves its ordering. The parameters that maximize \(\ell^*\) therefore also maximize \(\ell\), even though the reported totals differ. Either total can be positive because a density can exceed one; the relevant comparison is which parameter vector gives the higher value on the same data and scale.
That equivalence is enough for parameter estimation, but information criteria combine the likelihood with a penalty whose scale is fixed. The formulas
\[ AIC=-2\ell+2k, \qquad BIC=-2\ell+k\log(n) \]
use \(k\) to count the estimated quantities entering the penalty. Substituting \(\ell^*\) for \(\ell\) would double the relative weight assigned to fit while leaving that penalty unchanged. We therefore use the complete Gaussian value for both criteria. Although its constant cancels in comparisons on identical observations, retaining it also makes the reported likelihood conform to the usual definition.
The code reflects this distinction. The estimation functions nll_garch() and nll_garch_vt() return \(-\ell^*\) because R’s optim() searches for a minimum. After estimation, a returned value called nll can be put on the Gaussian scale using
We also calculate the Gaussian likelihood directly from the returns and fitted variances in gaussian_loglik(). Agreement between that calculation and the conversion provides a check on the implementation. Both use only dates with valid variances, since the shifted initialization leaves the first return without an associated variance value. Keeping this convention identical for the two estimators ensures that differences in fit come from their parameters rather than their observation sets.
During the search, a vector that violates the parameter restrictions receives a very large objective value. This keeps such a vector from being accepted as a better fit simply because the variance calculation has become invalid. The resulting conditional-likelihood procedure is developed more formally in Francq and Zakoian (2019); here we can trace each of its steps in the two functions below.
Code
gaussian_loglik <- function(u, v) {
valid <- is.finite(u) & is.finite(v) & v > 0
if (!any(valid)) {
stop("No valid observations are available for the Gaussian likelihood.")
}
list(
value = -0.5 * sum(
log(2 * pi) + log(v[valid]) + (u[valid]^2) / v[valid]
),
n = sum(valid),
valid = valid
)
}
nll_garch <- function(par, u) {
omega <- par[1]
alpha <- par[2]
beta <- par[3]
if (omega <= 0 || alpha < 0 || beta < 0 || (alpha + beta) >= 1) {
return(1e12)
}
v <- garch_variance(u, omega, alpha, beta)
valid <- !is.na(v) & v > 0
if (!any(valid)) {
return(1e12)
}
-sum(-log(v[valid]) - (u[valid]^2) / v[valid])
}3.6 Full MLE results
We first allow the likelihood to choose all three parameters jointly. Returning to the model equation makes clear what is being estimated:
\[ v_t=\omega+\alpha u_{t-1}^2+\beta v_{t-1}. \]
For each trial choice of \((\omega,\alpha,\beta)\), the recursion produces a variance path and the likelihood assigns it a score. We use L-BFGS-B to search for the vector with the lowest negative score. Its name refers to the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with bound constraints. As a quasi-Newton method, it learns about the objective’s curvature from changes in its slope and uses that information to choose the next search direction. It keeps a limited record of recent steps, avoiding the need to construct an exact second-derivative matrix at every iteration. In our call to optim(), the slopes themselves are approximated numerically because we do not supply a derivative function.
Its bound constraints are useful here because the parameters cannot take arbitrary values. We specify a positive lower bound for \(\omega\) and nonnegative lower bounds for \(\alpha\) and \(\beta\). However, a separate bound on each coefficient still permits some combinations with \(\alpha+\beta\ge1\), so the penalty in the objective handles that joint restriction. The reported convergence code tells us whether the search met its stopping rule. We will assess the solution further by checking its admissibility and comparing its likelihood with the variance-targeting fit.
Code
start <- c(omega = 4e-6, alpha = 0.2, beta = 0.7)
fit_full <- optim(
par = start,
fn = nll_garch,
u = u,
method = "L-BFGS-B",
lower = c(1e-12, 0, 0),
upper = c(Inf, 1 - 1e-6, 1 - 1e-6),
control = list(
parscale = c(1e-6, 0.1, 0.1),
factr = 1e7,
pgtol = 1e-10
)
)
if (fit_full$convergence != 0) {
stop("The scaled Full MLE optimization did not converge.")
}
theta_full <- fit_full$par
omega_full <- theta_full["omega"]
alpha_full <- theta_full["alpha"]
beta_full <- theta_full["beta"]
v_full <- garch_variance(u, omega_full, alpha_full, beta_full)
likelihood_full <- gaussian_loglik(u, v_full)
valid_full <- likelihood_full$valid
n_full <- likelihood_full$n
loglik_full <- likelihood_full$value
nll_full <- nll_garch(theta_full, u)
loglik_full_from_conversion <-
-nll_full / 2 - n_full / 2 * log(2 * pi)
stopifnot(isTRUE(all.equal(
loglik_full,
loglik_full_from_conversion,
tolerance = 1e-10
)))
k_full <- 3
aic_full <- -2 * loglik_full + 2 * k_full
bic_full <- -2 * loglik_full + log(n_full) * k_full
ab_full <- alpha_full + beta_full
half_life_full <- ifelse(ab_full < 1, log(0.5) / log(ab_full), NA_real_)
full_summary <- tibble(
metric = c(
"omega",
"alpha",
"beta",
"alpha + beta",
"half-life in trading days",
"logLik",
"n_obs",
"AIC",
"BIC"
),
value = c(
omega_full,
alpha_full,
beta_full,
ab_full,
half_life_full,
loglik_full,
n_full,
aic_full,
bic_full
)
)
full_summary |>
kable(
caption = "Full maximum likelihood GARCH(1,1) estimates with Gaussian log-likelihood.",
digits = 6,
col.names = c("Metric", "Value"),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Metric | Value |
|---|---|
| omega | 0.000004 |
| alpha | 0.223788 |
| beta | 0.747579 |
| alpha + beta | 0.971368 |
| half-life in trading days | 23.860302 |
| logLik | 4263.605601 |
| n_obs | 1257.000000 |
| AIC | -8521.211202 |
| BIC | -8505.801752 |
The estimates in the table also explain the scaling used in the optimization call. The intercept is of order \(10^{-6}\), while the two weights are of order tenths, so a numerically sensible step in one coordinate can be far too large or too small in another. The parscale argument allows the search to account for these different magnitudes. It ends with convergence code 0, and gaussian_loglik() then evaluates the fitted path on the scale used for logLik, AIC, and BIC. Inserting the coefficients back into the model gives the daily updating rule
\[ \hat v_t= 0.0000039818 + 0.22378836\,u_{t-1}^2 + 0.74757931\,\hat v_{t-1}. \]
We can now follow this rule through the first few observations. In the table below, day \(i\) counts from the first price in sp500, so the return on day 2 is u[1] and the initialized variance on day 3 is v_full[2]. The dashes show where an earlier observation is still missing. Once both the current return and its variance are available, they give the daily contribution in the last column, \(\ell_i^*=-\log(\hat v_i)-u_i^2/\hat v_i\). Adding these contributions over all defined rows, including those omitted from the display, reproduces the total score optimized above. Returns are shown in decimals, with variances in the corresponding squared units.
In the numerical substitutions below, equality signs describe the calculation carried out at full precision. Displayed inputs and results are rounded, so recalculating from the printed digits may differ in the last decimal places.
Code
garch_calculation <- function(prices, variances) {
daily <- prices |>
mutate(day = row_number(), variance = c(NA_real_, variances)) |>
select(date, day, price, return, variance)
daily$contribution <- NA_real_
valid <- is.finite(daily$return) & is.finite(daily$variance) & daily$variance > 0
daily$contribution[valid] <- -log(daily$variance[valid]) -
daily$return[valid]^2 / daily$variance[valid]
daily
}
compact_garch_calculation <- function(daily) {
display <- daily |>
mutate(
date = as.character(date), day = as.character(day),
price = sprintf("%.2f", price),
return = ifelse(is.na(return), "--", sprintf("%.10f", return)),
variance = ifelse(is.na(variance), "--", sprintf("%.10f", variance)),
contribution = ifelse(is.na(contribution), "--", sprintf("%.6f", contribution))
)
omitted <- as_tibble(setNames(rep(list("..."), ncol(display)), names(display)))
bind_rows(head(display, 6), omitted, tail(display, 2))
}
full_daily <- garch_calculation(sp500, v_full)
stopifnot(isTRUE(all.equal(
sum(full_daily$contribution, na.rm = TRUE), -nll_full, tolerance = 1e-10
)))
compact_garch_calculation(full_daily) |>
kable(
caption = "Daily GARCH(1,1) variance and likelihood calculations: full maximum likelihood.",
col.names = c("Date", "Day", "S&P 500", "\\(u_i\\)", "\\(\\hat v_i\\)", "\\(\\ell_i^*\\)"),
escape = FALSE, row.names = FALSE, align = c("l", rep("r", 5))
) |>
kable_styling(full_width = FALSE, latex_options = "HOLD_position")| Date | Day | S&P 500 | \(u_i\) | \(\hat v_i\) | \(\ell_i^*\) |
|---|---|---|---|---|---|
| 2015-07-10 | 1 | 2076.62 | -- | -- | -- |
| 2015-07-13 | 2 | 2099.60 | 0.0110660492 | -- | -- |
| 2015-07-14 | 3 | 2108.95 | 0.0044531592 | 0.0001224574 | 8.845808 |
| 2015-07-15 | 4 | 2107.40 | -0.0007349861 | 0.0000999663 | 9.205273 |
| 2015-07-16 | 5 | 2124.29 | 0.0080146804 | 0.0000788355 | 8.633348 |
| 2015-07-17 | 6 | 2126.64 | 0.0011061830 | 0.0000772927 | 9.452080 |
| ... | ... | ... | ... | ... | ... |
| 2020-07-08 | 1258 | 3169.94 | 0.0078274619 | 0.0001672895 | 8.329538 |
| 2020-07-09 | 1259 | 3152.05 | -0.0056436062 | 0.0001427553 | 8.631268 |
Day 3 is initialized as \(\hat v_3=u_2^2=(0.0110660492)^{2} = 0.0001224574\). The parameters first enter the variance calculation on day 4. Reading the preceding return and variance from the table gives
\[ \begin{aligned} \hat v_4 &=\hat\omega+\hat\alpha u_3^2+\hat\beta\hat v_3\\ &=0.0000039818 +0.22378836(0.0044531592)^2\\ &\quad+0.74757931(0.0001224574)\\ &= 0.0000999663. \end{aligned} \]
We then carry this newly calculated variance forward one day:
\[ \begin{aligned} \hat v_5 &=\hat\omega+\hat\alpha u_4^2+\hat\beta\hat v_4\\ &=0.0000039818 +0.22378836(-0.0007349861)^2\\ &\quad+0.74757931(0.0000999663)\\ &= 0.0000788355. \end{aligned} \]
The resulting values reproduce the fifth column for days 4 and 5, subject to the rounding used in the displayed substitutions; R retains full precision throughout. Notice that each update uses only the preceding row. The current return enters afterwards, when we evaluate the likelihood contribution for the variance already assigned to that date. This is how the same calculation can serve both estimation and forecasting without using a return before it has been observed.
The estimated weights tell us how long the effect of a large movement is expected to remain. While \(\alpha\) determines the immediate response and \(\beta\) carries forward the preceding variance, their sum governs the expected decay toward the long-run level. To express that persistence in trading days, we calculate the half-life:
\[ \text{half-life} = \frac{\log(0.5)}{\log(\alpha+\beta)}. \]
With an estimated persistence of 0.9714, it takes approximately 23.9 trading days for an initial excess variance to fall by half in expectation. For instance, a variance initially 0.00010 above its long-run level would be expected to lie about 0.00005 above it after one half-life. The quantity being halved is the deviation from the long-run variance; total volatility is obtained by taking the square root of the full variance at each date. New shocks can alter the actual path, but the half-life still gives a useful financial reading of the fitted dynamics: a longer half-life means that an episode of high risk continues to influence forecasts for longer. We next examine how these dynamics change when estimation is anchored to a specified long-run variance.
3.7 Variance targeting
The first estimator allowed all three parameters to adjust to the likelihood. Variance targeting retains the same GARCH(1,1) equation,
\[ v_t=\omega+\alpha u_{t-1}^2+\beta v_{t-1}, \]
but ties its long-run variance to a value \(V_L\) chosen before optimizing the weights. The stationary-variance identity then determines the intercept:
\[ \omega = V_L(1-\alpha-\beta). \]
We use the training sample variance of returns for \(V_L\). For each proposed pair \((\alpha,\beta)\), the equation above supplies \(\omega\), after which the same recursion and likelihood can be evaluated. The search therefore has two free coordinates. Its reduction in dimension comes from anchoring the model’s long-run level to the sample estimate, while allowing the likelihood to choose how quickly the variance responds and how persistent those responses are.
Code
nll_garch_vt <- function(par, u, long_run_variance) {
alpha <- par[1]
beta <- par[2]
if (alpha < 0 || beta < 0 || (alpha + beta) >= 1) {
return(1e12)
}
omega <- long_run_variance * (1 - alpha - beta)
v <- garch_variance(u, omega, alpha, beta)
valid <- !is.na(v) & v > 0
if (!any(valid)) {
return(1e12)
}
-sum(-log(v[valid]) - (u[valid]^2) / v[valid])
}
fit_vt <- optim(
par = c(alpha = 0.2, beta = 0.7),
fn = nll_garch_vt,
u = u,
long_run_variance = long_run_variance,
method = "L-BFGS-B",
lower = c(0, 0),
upper = c(1 - 1e-6, 1 - 1e-6)
)
if (fit_vt$convergence != 0) {
stop("The variance targeting optimization did not converge.")
}
alpha_vt <- fit_vt$par["alpha"]
beta_vt <- fit_vt$par["beta"]
omega_vt <- long_run_variance * (1 - alpha_vt - beta_vt)
v_vt <- garch_variance(u, omega_vt, alpha_vt, beta_vt)
likelihood_vt <- gaussian_loglik(u, v_vt)
valid_vt <- likelihood_vt$valid
n_vt <- likelihood_vt$n
loglik_vt <- likelihood_vt$value
nll_vt <- nll_garch_vt(fit_vt$par, u, long_run_variance)
loglik_vt_from_conversion <-
-nll_vt / 2 - n_vt / 2 * log(2 * pi)
stopifnot(isTRUE(all.equal(
loglik_vt,
loglik_vt_from_conversion,
tolerance = 1e-10
)))
theta_vt_in_full <- c(
omega = unname(omega_vt),
alpha = unname(alpha_vt),
beta = unname(beta_vt)
)
nll_vt_from_full <- nll_garch(theta_vt_in_full, u)
stopifnot(isTRUE(all.equal(
nll_vt,
nll_vt_from_full,
tolerance = 1e-10
)))
if (!identical(which(valid_full), which(valid_vt))) {
stop("Full MLE and variance targeting use different likelihood observations.")
}
if (loglik_full + 1e-8 < loglik_vt) {
stop("Full MLE has a lower Gaussian log-likelihood than variance targeting.")
}
k_vt <- 2
aic_vt <- -2 * loglik_vt + 2 * k_vt
bic_vt <- -2 * loglik_vt + log(n_vt) * k_vt
ab_vt <- alpha_vt + beta_vt
half_life_vt <- ifelse(ab_vt < 1, log(0.5) / log(ab_vt), NA_real_)
vt_summary <- tibble(
metric = c(
"omega",
"alpha",
"beta",
"alpha + beta",
"half-life in trading days",
"logLik",
"n_obs",
"AIC",
"BIC"
),
value = c(
omega_vt,
alpha_vt,
beta_vt,
ab_vt,
half_life_vt,
loglik_vt,
n_vt,
aic_vt,
bic_vt
)
)
vt_summary |>
kable(
caption = "Variance targeting GARCH(1,1) estimates with Gaussian log-likelihood.",
digits = 6,
col.names = c("Metric", "Value"),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Metric | Value |
|---|---|
| omega | 0.000004 |
| alpha | 0.226349 |
| beta | 0.747038 |
| alpha + beta | 0.973387 |
| half-life in trading days | 25.696953 |
| logLik | 4263.596586 |
| n_obs | 1257.000000 |
| AIC | -8523.193171 |
| BIC | -8512.920205 |
With \(V_L=0.0001490365\), the optimized weights imply \(\hat\omega=V_L(1-\hat\alpha-\hat\beta) =0.0000039664\). Although the table reports three coefficients, only \(\alpha\) and \(\beta\) have been chosen directly by the optimizer. Their values and the long-run input together give the fitted rule
\[ \hat v_t= 0.0000039664 + 0.22634899\,u_{t-1}^2 + 0.74703765\,\hat v_{t-1}. \]
Applying this rule to the same price and return sequence lets us see where the two procedures begin to differ. Their initialized variance on day 3 is identical, since it depends only on the first available return. Differences first appear on day 4, when the fitted coefficients enter the recursive update. Keeping the dates and initialization fixed makes the following table directly comparable with its full maximum likelihood counterpart.
Code
vt_daily <- garch_calculation(sp500, v_vt)
stopifnot(isTRUE(all.equal(
sum(vt_daily$contribution, na.rm = TRUE), -nll_vt, tolerance = 1e-10
)))
compact_garch_calculation(vt_daily) |>
kable(
caption = "Daily GARCH(1,1) variance and likelihood calculations: variance targeting.",
col.names = c("Date", "Day", "S&P 500", "\\(u_i\\)", "\\(\\hat v_i\\)", "\\(\\ell_i^*\\)"),
escape = FALSE, row.names = FALSE, align = c("l", rep("r", 5))
) |>
kable_styling(full_width = FALSE, latex_options = "HOLD_position")| Date | Day | S&P 500 | \(u_i\) | \(\hat v_i\) | \(\ell_i^*\) |
|---|---|---|---|---|---|
| 2015-07-10 | 1 | 2076.62 | -- | -- | -- |
| 2015-07-13 | 2 | 2099.60 | 0.0110660492 | -- | -- |
| 2015-07-14 | 3 | 2108.95 | 0.0044531592 | 0.0001224574 | 8.845808 |
| 2015-07-15 | 4 | 2107.40 | -0.0007349861 | 0.0000999353 | 9.205582 |
| 2015-07-16 | 5 | 2124.29 | 0.0080146804 | 0.0000787441 | 8.633562 |
| 2015-07-17 | 6 | 2126.64 | 0.0011061830 | 0.0000773307 | 9.451596 |
| ... | ... | ... | ... | ... | ... |
| 2020-07-08 | 1258 | 3169.94 | 0.0078274619 | 0.0001684868 | 8.325010 |
| 2020-07-09 | 1259 | 3152.05 | -0.0056436062 | 0.0001437005 | 8.626136 |
Using the variance-targeting coefficients, the first recursive update is
\[ \begin{aligned} \hat v_4 &=\hat\omega+\hat\alpha u_3^2+\hat\beta\hat v_3\\ &=0.0000039664 +0.22634899(0.0044531592)^2\\ &\quad+0.74703765(0.0001224574)\\ &= 0.0000999353. \end{aligned} \]
For day 5, we use the return on day 4 and the variance just obtained:
\[ \begin{aligned} \hat v_5 &=\hat\omega+\hat\alpha u_4^2+\hat\beta\hat v_4\\ &=0.0000039664 +0.22634899(-0.0007349861)^2\\ &\quad+0.74703765(0.0000999353)\\ &= 0.0000787441. \end{aligned} \]
These substitutions reproduce the fifth column, up to display rounding, using exactly the same timing as the full MLE calculation. What has changed is the restriction on the long-run level. Fixing that level makes the estimation problem smaller and can improve numerical stability, provided that the sample variance is a reasonable anchor for the period under study. The comparison below examines how much fit is lost by imposing this restriction.
It also requires us to state how the estimated quantities are counted in AIC and BIC. Our main calculation treats \(V_L\) as fixed during optimization and uses \(k=2\) for the two fitted weights. Yet \(V_L\) itself was obtained from the training sample, so we also report a count of three that includes this sample-based input. Showing both conventions allows us to separate the empirical difference in fit from the choice of complexity penalty.
Before comparing the objectives, it is useful to place the two daily variance paths alongside the same observed returns. The table below keeps the dates, day numbers, and price levels used in the individual calculation tables, and places the two estimates of \(\hat v_i\) in adjacent columns. Their shared initialization makes day 3 identical. From day 4 onward, the small differences reflect the fitted coefficients, since both procedures receive exactly the same return history. These columns report conditional variances; taking their square roots would give conditional volatilities.
Code
daily_comparison <- compact_garch_calculation(full_daily) |>
select(date, day, price, return, variance_full = variance) |>
mutate(variance_targeting = compact_garch_calculation(vt_daily)$variance)
stopifnot(identical(full_daily$date, vt_daily$date))
daily_comparison |>
kable(
caption = "Daily GARCH(1,1) conditional variances: full maximum likelihood and variance targeting.",
col.names = c("Date", "Day", "S&P 500", "\\(u_i\\)",
"\\(\\hat v_i\\) (Full MLE)", "\\(\\hat v_i\\) (Variance targeting)"),
escape = FALSE, row.names = FALSE, align = c("l", rep("r", 5))
) |>
kable_styling(full_width = FALSE, latex_options = "HOLD_position")| Date | Day | S&P 500 | \(u_i\) | \(\hat v_i\) (Full MLE) | \(\hat v_i\) (Variance targeting) |
|---|---|---|---|---|---|
| 2015-07-10 | 1 | 2076.62 | -- | -- | -- |
| 2015-07-13 | 2 | 2099.60 | 0.0110660492 | -- | -- |
| 2015-07-14 | 3 | 2108.95 | 0.0044531592 | 0.0001224574 | 0.0001224574 |
| 2015-07-15 | 4 | 2107.40 | -0.0007349861 | 0.0000999663 | 0.0000999353 |
| 2015-07-16 | 5 | 2124.29 | 0.0080146804 | 0.0000788355 | 0.0000787441 |
| 2015-07-17 | 6 | 2126.64 | 0.0011061830 | 0.0000772927 | 0.0000773307 |
| ... | ... | ... | ... | ... | ... |
| 2020-07-08 | 1258 | 3169.94 | 0.0078274619 | 0.0001672895 | 0.0001684868 |
| 2020-07-09 | 1259 | 3152.05 | -0.0056436062 | 0.0001427553 | 0.0001437005 |
3.8 Objective surface and estimator comparison
Having reduced the optimization to two free parameters, we can visualize the search itself. Each point in the \((\alpha,\beta)\) plane determines an intercept through variance targeting and hence a complete variance path and likelihood value. The contours join pairs giving the same fit, while the marked solution identifies the pair returned by optim(). We plot the Gaussian log-likelihood so that the contour labels use the same scale as the reported results. Converting from the proportional score preserves both the contour shapes and the optimum’s location; the selected axis limits let us examine their behavior near that solution.
Code
alpha_grid <- seq(0.22, 0.30, length.out = 100)
beta_grid <- seq(0.66, 0.76, length.out = 100)
objective_surface <- outer(alpha_grid, beta_grid, Vectorize(function(alpha, beta) {
if (alpha < 0 || beta < 0 || (alpha + beta) >= 0.999) {
return(NA_real_)
}
nll_garch_vt(c(alpha, beta), u = u, long_run_variance = long_run_variance)
}))
proportional_levels <- c(10700, 10750, 10800, 10820, 10830, 10835)
loglik_surface <-
-objective_surface / 2 - n_vt / 2 * log(2 * pi)
contour_levels <-
proportional_levels / 2 - n_vt / 2 * log(2 * pi)
contour_lines <- grDevices::contourLines(
x = beta_grid,
y = alpha_grid,
z = t(loglik_surface),
levels = contour_levels
)
contour_tbl <- bind_rows(lapply(seq_along(contour_lines), function(id) {
tibble(
contour_id = id,
beta = contour_lines[[id]]$x,
alpha = contour_lines[[id]]$y,
log_likelihood = contour_lines[[id]]$level
)
}))
label_candidates <- contour_tbl |>
filter(
beta >= 0.662,
beta <= 0.756,
alpha >= 0.225,
alpha <= 0.285
)
level_tbl <- tibble(
log_likelihood = contour_levels,
label_beta_target = c(0.715, 0.730, 0.743, 0.748, 0.752, 0.755)
)
contour_label_tbl <- label_candidates |>
inner_join(level_tbl, by = "log_likelihood") |>
group_by(log_likelihood) |>
slice_min(abs(beta - label_beta_target), n = 1, with_ties = FALSE) |>
ungroup() |>
mutate(label = scales::number(log_likelihood, accuracy = 0.1, big.mark = ""))
solution_tbl <- tibble(
beta = beta_vt,
alpha = alpha_vt
)
contour_tbl |>
ggplot(aes(beta, alpha)) +
geom_path(aes(group = contour_id), color = "#1B4E8A", linewidth = 0.5) +
geom_label(
data = contour_label_tbl,
aes(label = label),
color = "#2F2F2F",
fill = "white",
label.size = 0,
size = 3,
alpha = 0.9
) +
geom_hline(yintercept = alpha_vt, color = "#B22222", linetype = 2) +
geom_vline(xintercept = beta_vt, color = "#B22222", linetype = 2) +
geom_point(
data = solution_tbl,
color = "#B22222",
size = 2.8
) +
coord_cartesian(xlim = c(0.66, 0.76), ylim = c(0.22, 0.30), expand = FALSE) +
labs(
x = expression(beta),
y = expression(alpha),
caption = "Contour labels show Gaussian log-likelihood values; higher values indicate better fit."
) +
theme_minimal(base_size = 12) +
theme(legend.position = "none")
As the contour labels increase toward the red point, the fit improves and the negative score minimized by optim() decreases. The shape also helps us understand the search: a long, flat contour would indicate that some changes in \(\alpha\) can be offset by changes in \(\beta\) with little loss of fit. To assess the cost of anchoring the long-run variance, we now compare this solution with the one obtained when \(\omega\) is estimated freely.
Code
comparison_tbl <- bind_rows(
full_summary |> mutate(model = "Full MLE"),
vt_summary |> mutate(model = "Variance targeting")
) |>
pivot_wider(names_from = model, values_from = value) |>
select(metric, `Full MLE`, `Variance targeting`)
comparison_tbl |>
kable(
caption = "Full MLE and variance targeting comparison.",
digits = 6,
col.names = c("Metric", "Full MLE", "Variance targeting"),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Metric | Full MLE | Variance targeting |
|---|---|---|
| omega | 0.000004 | 0.000004 |
| alpha | 0.223788 | 0.226349 |
| beta | 0.747579 | 0.747038 |
| alpha + beta | 0.971368 | 0.973387 |
| half-life in trading days | 23.860302 | 25.696953 |
| logLik | 4263.605601 | 4263.596586 |
| n_obs | 1257.000000 | 1257.000000 |
| AIC | -8521.211202 | -8523.193171 |
| BIC | -8505.801752 | -8512.920205 |
Code
optimization_audit_tbl <- tibble(
model = c("Full MLE", "Variance targeting"),
convergence_code = c(fit_full$convergence, fit_vt$convergence),
proportional_criterion = c(-fit_full$value, -fit_vt$value),
gaussian_logLik = c(loglik_full, loglik_vt)
)
optimization_audit_tbl |>
kable(
caption = "Optimization and likelihood audit.",
digits = 6,
col.names = c(
"Model",
"Convergence code",
"Proportional criterion",
"Gaussian logLik"
),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Model | Convergence code | Proportional criterion | Gaussian logLik |
|---|---|---|---|
| Full MLE | 0 | 10837.42 | 4263.606 |
| Variance targeting | 0 | 10837.40 | 4263.597 |
Code
aic_vt_k3 <- -2 * loglik_vt + 2 * 3
bic_vt_k3 <- -2 * loglik_vt + log(n_vt) * 3
parameter_count_tbl <- tibble(
specification = c(
"Full MLE",
"Variance targeting: conditional count",
"Variance targeting: count sample-based V_L"
),
quantities_counted = c(
"omega, alpha, beta",
"alpha, beta",
"V_L, alpha, beta"
),
k = c(3, 2, 3),
AIC = c(aic_full, aic_vt, aic_vt_k3),
BIC = c(bic_full, bic_vt, bic_vt_k3)
)
parameter_count_tbl |>
kable(
caption = "AIC and BIC under two parameter counts for variance targeting.",
digits = 6,
col.names = c(
"Specification",
"Quantities counted",
"k",
"AIC",
"BIC"
),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Specification | Quantities counted | k | AIC | BIC |
|---|---|---|---|---|
| Full MLE | omega, alpha, beta | 3 | -8521.211 | -8505.802 |
| Variance targeting: conditional count | alpha, beta | 2 | -8523.193 | -8512.920 |
| Variance targeting: count sample-based V_L | V_L, alpha, beta | 3 | -8521.193 | -8505.784 |
The two likelihoods are evaluated on the same 1257 observations, so their difference measures fit on a common basis. Full MLE reaches 4 263.606, compared with 4 263.597 for variance targeting. The gain from estimating the intercept freely is therefore 0.009. The ordering also provides a numerical check: the variance-targeting parameter vector is available to the unrestricted model, so the unrestricted maximum cannot be lower. Both searches report convergence code zero, and the direct Gaussian evaluation agrees with the conversion from the proportional score. Together, these checks support reading the small difference as a feature of the fitted models rather than a mismatch in calculations.
The similarity in fit carries over to the economic interpretation. Full MLE has persistence 0.971 and a half-life of 23.9 trading days, compared with 0.973 and 25.7 days under variance targeting. Both therefore describe a persistent response to volatility shocks, with the freely estimated model allowing that response to decay slightly faster. Anchoring the long-run variance changes the fit modestly without changing this broad reading of the dynamics.
The AIC and BIC comparison adds a penalty for the estimated quantities, which is where the treatment of \(V_L\) enters. Full MLE estimates three coefficients directly and therefore has \(k=3\). Variance targeting optimizes two weights and obtains \(\omega=V_L(1-\alpha-\beta)\), giving \(k=2\) when the long-run input is treated as fixed during that step. Including the fact that \(V_L\) was itself estimated from these returns gives a count of three instead. In either case, \(k\) counts quantities entering the penalty; it is unrelated to the lag orders denoted by \((1,1)\).
Under the conditional count of two, the smaller penalty gives variance targeting the lower AIC and BIC. Counting \(V_L\) puts both fits at \(k=3\), so Full MLE’s slightly higher likelihood then gives it the lower criteria, by 0.018 points each. Changing the count leaves the fitted coefficients, likelihood, variance paths, and forecasts untouched: the switch in preference comes entirely from the penalty. Reporting both conventions makes that dependence visible and puts the very small equal-count difference in perspective.
The information criteria thus give little reason to expect substantially different descriptions of historical risk. To see what the estimates imply in financial terms, we turn from the parameter comparison to the volatility paths and forecasts they generate.
3.9 Conditional volatility through time and forecasts
Taking the square root of each fitted variance translates the updating rule into a daily measure of risk that can be read alongside the returns. We use the variance-targeting fit for this illustration because its long-run level is anchored explicitly to the sample variance, giving a natural reference for the plot. The full MLE path can be reconstructed in the same way. Such time-varying risk measures underpin the risk-management, scenario-analysis, and derivatives applications discussed by Hull (2022) and McNeil et al. (2015).
Code
volatility_tbl <- returns |>
mutate(
variance_vt = v_vt,
volatility_vt = sqrt(variance_vt),
sample_volatility = sqrt(long_run_variance)
)
volatility_tbl |>
select(date, return, variance_vt, volatility_vt) |>
slice(c(1:6, (n() - 1):n())) |>
kable(
caption = "Conditional variance and volatility from variance targeting.",
col.names = c("Date", "\\(u_t\\)", "\\(\\hat v_t\\)", "\\(\\sqrt{\\hat v_t}\\)"),
escape = FALSE,
digits = 8,
format.args = list(scientific = FALSE),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Date | \(u_t\) | \(\hat v_t\) | \(\sqrt{\hat v_t}\) |
|---|---|---|---|
| 2015-07-13 | 0.01106605 | NA | NA |
| 2015-07-14 | 0.00445316 | 0.00012246 | 0.01106605 |
| 2015-07-15 | -0.00073499 | 0.00009994 | 0.00999677 |
| 2015-07-16 | 0.00801468 | 0.00007874 | 0.00887379 |
| 2015-07-17 | 0.00110618 | 0.00007733 | 0.00879379 |
| 2015-07-20 | 0.00077123 | 0.00006201 | 0.00787479 |
| 2020-07-08 | 0.00782746 | 0.00016849 | 0.01298025 |
| 2020-07-09 | -0.00564361 | 0.00014370 | 0.01198752 |
Code
volatility_tbl |>
ggplot(aes(date, 100 * volatility_vt)) +
geom_line(color = "#1B4E8A", linewidth = 0.6) +
geom_hline(
aes(yintercept = 100 * sample_volatility),
color = "#B22222",
linetype = 2
) +
labs(
x = "Date",
y = "Daily volatility",
caption = "Dashed line: sample daily volatility."
) +
scale_y_continuous(labels = scales::label_percent(accuracy = 0.1, scale = 1)) +
theme_minimal(base_size = 12)
The horizontal line averages over the full sample, whereas the fitted series responds to the changing size of the observed returns. Its sharp increase during the 2020 market stress episode shows how the recursion translates successive large movements into a higher assessment of risk. During quieter periods, smaller squared returns gradually bring that assessment down.
To move from filtering the historical data to forecasting, first recall the GARCH updating rule at an ordinary date \(t\):
\[ v_t=\omega+\alpha u_{t-1}^2+\beta v_{t-1}. \]
The variance assigned to day \(t\) uses the return and variance from day \(t-1\). At the final observation, \(T\), advancing this same rule by one period means replacing \(t\) by \(T+1\): the preceding return becomes \(u_T\) and the preceding variance becomes \(v_T\). With fitted parameters, we already have the return and the fitted variance needed for this update. Writing its result as \(\hat v_{T+1\mid T}\) makes the information date explicit: this is tomorrow’s variance forecast using information through today. Thus the first forecast requires no assumption about a future realized shock:
\[ \hat{v}_{T+1\mid T} = \omega + \alpha u_T^2 + \beta v_T. \]
The next step introduces uncertainty because \(u_{T+1}\) has not yet been observed. Under the zero-mean model, its expected square equals its conditional variance. Using that expectation in the recursion combines the two weights into \(\alpha+\beta\), and repeating the argument gives the forecast at any horizon \(h\):
\[ \hat{v}_{T+h\mid T} = \bar{v} + (\alpha+\beta)^{h-1} \left( \hat{v}_{T+1\mid T}-\bar{v} \right), \qquad h \ge 1. \]
This expression begins at the one-step forecast and reduces its deviation from the long-run variance, \(\bar v\), by a factor of \(\alpha+\beta\) for each additional trading day. In the calculation below, the parameters and final variance are replaced by their fitted values. Taking square roots then puts the forecasts on the same volatility scale as the historical plot.
Code
last_return <- tail(u, 1)
last_variance <- tail(v_vt, 1)
v_next <- omega_vt + alpha_vt * last_return^2 + beta_vt * last_variance
long_run_vt <- omega_vt / (1 - alpha_vt - beta_vt)
horizon <- 20
forecast_variance <- long_run_vt +
(alpha_vt + beta_vt)^(seq_len(horizon) - 1) * (v_next - long_run_vt)
forecast_tbl <- tibble(
horizon = seq_len(horizon),
variance_forecast = forecast_variance,
volatility_forecast = sqrt(forecast_variance)
)
forecast_tbl |>
kable(
caption = "Twenty-day conditional variance and volatility forecasts from variance targeting.",
col.names = c("Horizon", "\\(\\hat v_{T+h\\mid T}\\)", "\\(\\sqrt{\\hat v_{T+h\\mid T}}\\)"),
escape = FALSE,
digits = 8,
format.args = list(scientific = FALSE),
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Horizon | \(\hat v_{T+h\mid T}\) | \(\sqrt{\hat v_{T+h\mid T}}\) |
|---|---|---|
| 1 | 0.00011853 | 0.01088694 |
| 2 | 0.00011934 | 0.01092416 |
| 3 | 0.00012013 | 0.01096028 |
| 4 | 0.00012090 | 0.01099532 |
| 5 | 0.00012165 | 0.01102932 |
| 6 | 0.00012237 | 0.01106232 |
| 7 | 0.00012308 | 0.01109435 |
| 8 | 0.00012378 | 0.01112543 |
| 9 | 0.00012445 | 0.01115560 |
| 10 | 0.00012510 | 0.01118489 |
| 11 | 0.00012574 | 0.01121333 |
| 12 | 0.00012636 | 0.01124095 |
| 13 | 0.00012696 | 0.01126776 |
| 14 | 0.00012755 | 0.01129380 |
| 15 | 0.00012812 | 0.01131909 |
| 16 | 0.00012868 | 0.01134365 |
| 17 | 0.00012922 | 0.01136750 |
| 18 | 0.00012975 | 0.01139068 |
| 19 | 0.00013026 | 0.01141319 |
| 20 | 0.00013076 | 0.01143506 |
Code
set.seed(20260923)
simulation_count <- 50000L
simulated_variance <- matrix(NA_real_, nrow = horizon, ncol = simulation_count)
simulated_variance[1, ] <- v_next
for (h in 2:horizon) {
simulated_return <- sqrt(simulated_variance[h - 1, ]) * rnorm(simulation_count)
simulated_variance[h, ] <- omega_vt + alpha_vt * simulated_return^2 +
beta_vt * simulated_variance[h - 1, ]
}
forecast_interval <- forecast_tbl |>
mutate(
lower = apply(sqrt(simulated_variance), 1, quantile, probs = 0.025),
upper = apply(sqrt(simulated_variance), 1, quantile, probs = 0.975),
trading_day = nrow(returns) + horizon
)
# Trading-day positions avoid inventing future prices or calendar observations.
history_plot <- volatility_tbl |>
mutate(trading_day = row_number())
forecast_line <- bind_rows(
tibble(trading_day = nrow(returns), volatility_forecast = sqrt(last_variance)),
forecast_interval |> select(trading_day, volatility_forecast)
)
year_ticks <- history_plot |>
group_by(year = format(date, "%Y")) |>
slice_head(n = 1) |>
ungroup()
zoom_history <- history_plot |> filter(date >= as.Date("2020-01-01"))
month_ticks <- zoom_history |>
group_by(month = format(date, "%Y-%m")) |>
slice_head(n = 1) |>
ungroup()
panel_levels <- c("Full history", "From January 2020")
history_panels <- bind_rows(
history_plot |> mutate(panel = panel_levels[1]),
zoom_history |> mutate(panel = panel_levels[2])
) |>
mutate(panel = factor(panel, levels = panel_levels))
interval_panels <- bind_rows(lapply(panel_levels, function(view) {
forecast_interval |> mutate(panel = factor(view, levels = panel_levels))
}))
line_panels <- bind_rows(lapply(panel_levels, function(view) {
forecast_line |> mutate(panel = factor(view, levels = panel_levels))
}))
forecast_end <- nrow(returns) + horizon
ggplot() +
geom_line(
data = history_panels,
aes(trading_day, 100 * volatility_vt, color = "Historical"), linewidth = 0.6
) +
geom_ribbon(
data = interval_panels,
aes(trading_day, ymin = 100 * lower, ymax = 100 * upper),
fill = "#CC7A00", alpha = 0.20
) +
geom_line(
data = line_panels,
aes(trading_day, 100 * volatility_forecast, color = "Forecast"), linewidth = 0.9
) +
geom_vline(xintercept = nrow(returns), linetype = 3, color = "gray50") +
geom_hline(
yintercept = 100 * sqrt(long_run_vt),
color = "#B22222", linetype = 2
) +
scale_color_manual(values = c("Historical" = "#1B4E8A", "Forecast" = "#CC7A00")) +
scale_x_continuous(
breaks = function(limits) {
ticks <- if (diff(limits) > 500) year_ticks$trading_day else month_ticks$trading_day
c(ticks, forecast_end)
},
labels = function(ticks) {
date_labels <- if (diff(range(ticks)) > 500) "%Y" else "%b"
labels <- format(history_plot$date[match(ticks, history_plot$trading_day)], date_labels)
labels[ticks == forecast_end] <- "T + 20"
labels
}
) +
facet_wrap(~ panel, ncol = 1, scales = "free_x") +
labs(
x = "Trading-day sequence (T is the final training date)",
y = "Daily volatility", color = NULL,
caption = "Dashed horizontal line: long-run volatility. Dotted vertical line: forecast origin."
) +
scale_y_continuous(labels = scales::label_percent(accuracy = 0.1, scale = 1)) +
theme_minimal(base_size = 12) +
theme(legend.position = "right")
The upper panel places the forecast in the full historical record, while the lower panel enlarges the period from January 2020 onward so that the orange segment and its shaded interval can be examined. Both panels use the same vertical scale and the same values; only the horizontal window changes. Observations are spaced by trading day, with years or months marking their position and \(T+20\) indicating the final forecast horizon.
In the enlarged view, the smooth forecast is easier to distinguish from the irregular historical path. Each historical update incorporates a newly observed shock, so the fitted volatility rises and falls as market movements change. At \(T\), all later shocks are still unknown. Averaging over them leaves the regular decay toward the long-run variance described by the forecast equation. Because the last fitted variance is already relatively close to that level, the projected volatility changes only gradually over these 20 days. Its smoothness reflects the information available when it is made.
The plotted line takes the square root of expected future variance, giving the standard deviation of the future daily return conditional on information at \(T\) under the zero-mean model. That choice keeps the forecast on the same percentage scale as historical volatility. It also identifies precisely what is being plotted: taking an expectation and taking a square root in the opposite order would generally give a different quantity.
The smooth line leaves open how much conditional volatility could vary as future shocks arrive. To explore that question, we simulate 50,000 possible paths using the fitted parameters and independent standard-normal shocks, following the approach in Sheppard (n.d.). On each path, a shock is multiplied by its conditional standard deviation to generate a return, which then feeds into the next variance update. Collecting the resulting volatilities at a given horizon lets us locate their 2.5th and 97.5th percentiles. Repeating this across horizons produces the shaded pointwise 95% interval. The seed fixes the simulation draws for reproducibility and plays no role in fitting or selecting the model.
Every path begins with the same one-step variance because the final training return and variance are already known. The interval consequently has zero width at that step. Once simulated shocks enter subsequent updates, the paths separate and their distribution can become asymmetric, which is why we construct the band from percentiles. Since the coefficients are held fixed, it represents uncertainty from future shocks under the fitted model and excludes uncertainty about the coefficient estimates. Its 95% coverage applies separately at each horizon; it neither describes an interval for returns nor guarantees coverage of a whole volatility path. These distinctions will also be important when we compare the forecast with what happens after the training period.
3.10 Standardized Residual Diagnostics
Before evaluating the forecasts, we can use the training sample to check whether the fitted model has accounted for the dependence that motivated it. The return plot showed periods of unusually large movements; if the variance updates capture that behavior, returns measured relative to their fitted risk should exhibit less clustering. This suggests dividing each observed return by the conditional standard deviation assigned to it:
\[ z_t = \frac{u_t}{\sqrt{\hat{v}_t}}. \]
Here \(u_t\) is the residual from our zero-conditional-mean equation, and division by \(\sqrt{\hat v_t}\) removes its estimated time-varying scale. A value of \(z_t=2\) means that the return was two fitted standard deviations above zero; \(z_t=-2\) means two below. A movement during a quiet period can therefore be judged on the same basis as one during a turbulent period. This construction uses the return because the data contain no observed “true conditional volatility” from which we could subtract the fitted volatility. Under the correctly specified model with known parameters, the standardized shocks would have conditional mean zero and variance one, and the additional normality assumption would make them standard normal.
Standardization gives us two related ways to examine what the model has left unexplained. If \(z_t\) remains autocorrelated, earlier returns still contain linear information about later ones, which can point to a weakness in the zero-mean specification. If its square remains autocorrelated, large standardized movements still tend to follow other large movements, suggesting that the variance equation has left some clustering behind. Examining both series separates these mean and variance concerns, which could easily be confused by looking only at the original returns.
For each series, we use the Ljung-Box statistic to assess the first 20 autocorrelations jointly. A small p-value gives evidence against their being zero, complementing the individual lag patterns in the plots. We apply the same calculation on the same valid dates to both estimators, using Box.test() with fitdf = 0. The reported p-values should be read as approximate diagnostics: the usual chi-squared reference distribution does not fully account for estimating a GARCH model, and formal residual inference may require a model-specific adjustment or bootstrap. Even favorable results would address only the dependence examined here, leaving normality, independence, and forecasting performance to separate assessments.
Code
diag_tbl <- bind_rows(
tibble(date = returns$date, return = u, variance = v_full, model = "Full MLE"),
tibble(date = returns$date, return = u, variance = v_vt, model = "Variance targeting")
) |>
filter(is.finite(variance), variance > 0) |>
mutate(
z = return / sqrt(variance), z_squared = z^2,
model = factor(model, levels = c("Full MLE", "Variance targeting"))
) |>
pivot_longer(c(z, z_squared), names_to = "series", values_to = "value")
diag_tests <- diag_tbl |>
group_by(model, series) |>
group_modify(function(.x, .y) {
test <- stats::Box.test(.x$value, lag = 20, type = "Ljung-Box", fitdf = 0)
tibble(statistic = unname(test$statistic), p_value = test$p.value)
}) |>
ungroup()
diag_tests |>
mutate(series = ifelse(series == "z", "\\(z_t\\)", "\\(z_t^2\\)")) |>
kable(
caption = "Ljung-Box diagnostics at 20 lags for both GARCH estimators (approximate p-values).",
col.names = c("Estimator", "Series", "Statistic", "p-value"), escape = FALSE,
digits = 4,
row.names = FALSE
) |>
kable_styling(latex_options = "HOLD_position")| Estimator | Series | Statistic | p-value |
|---|---|---|---|
| Full MLE | \(z_t\) | 13.0077 | 0.8771 |
| Full MLE | \(z_t^2\) | 14.0276 | 0.8291 |
| Variance targeting | \(z_t\) | 12.9968 | 0.8775 |
| Variance targeting | \(z_t^2\) | 14.0188 | 0.8295 |
Code
acf_tbl <- diag_tbl |>
group_by(model, series) |>
group_modify(function(.x, .y) {
ac <- stats::acf(.x$value, plot = FALSE, lag.max = 20)
tibble(lag = as.numeric(ac$lag[-1]), acf = as.numeric(ac$acf[-1]),
limit = 1.96 / sqrt(nrow(.x)))
}) |>
ungroup() |>
mutate(series = factor(series, levels = c("z", "z_squared"),
labels = c("Standardized returns", "Squared standardized returns")))
acf_tbl |>
ggplot(aes(lag, acf)) +
geom_col(width = 0.08, fill = "#1B4E8A") +
geom_hline(data = distinct(acf_tbl, model, series, limit),
aes(yintercept = limit), linetype = 2, color = "#B22222") +
geom_hline(data = distinct(acf_tbl, model, series, limit),
aes(yintercept = -limit), linetype = 2, color = "#B22222") +
facet_grid(series ~ model) +
scale_x_continuous(breaks = c(1, 5, 10, 15, 20)) +
labs(x = "Lag", y = "ACF") +
theme_minimal(base_size = 12) +
theme(
panel.spacing.x = grid::unit(1.2, "cm"),
panel.spacing.y = grid::unit(0.6, "cm"),
panel.border = element_rect(color = "grey65", fill = NA, linewidth = 0.5),
strip.background = element_rect(fill = "grey95", color = NA)
)
The columns place the two estimators side by side, while the rows distinguish standardized returns from their squares. Their similar patterns agree with the similar parameter estimates found earlier. Some individual spikes cross the dashed pointwise reference limits, which can happen by chance when many lags are examined. The joint Ljung-Box checks help put those crossings in context: all four approximate p-values exceed 0.05, with the smallest equal to 0.829. At this 20-lag horizon, neither estimator leaves clear evidence of the linear dependence being tested in standardized returns or their squares.
These results support the fitted dynamics on the training sample, but we still need to examine predictions on later observations. The dashed lines here assess sampling variation in autocorrelations; the band in Section 3.9 describes possible future conditional volatilities. Neither calculation, by itself, tells us how accurate a variance forecast will be when confronted with new returns.
3.11 Out-of-sample forecast evaluation
So far, the observed returns have served both to estimate the parameters and to assess the residuals produced by that fit. A forecast faces a further test: it must describe observations that were unavailable when it was made. We create that separation by stopping estimation at the existing final training date, 2020-07-09, and evaluating forecasts against the next 20 S&P 500 returns. The additional data enter only this evaluation; all earlier estimates and calculations continue to use the original sample.
The practical question is what each estimator would have predicted before each test return became available, and how those predictions compare with the movements subsequently observed. We examine three forecasting policies for each estimator. The first commits to all twenty forecasts at the end of training. The second updates tomorrow’s forecast each day using the latest return while retaining the original coefficients. The third also re-estimates the coefficients before issuing each daily forecast. Keeping these policies distinct lets us separate the contribution of fresh return information from the additional effect of learning new parameter values.
We obtain the next 20 closing observations from the same provider and calculate the first test return relative to the final training close. This retains the price change across the boundary while keeping its return entirely in the test period. The first two policies hold both parameter vectors and the variance-targeting input \(V_L\) at their original training values. The third uses an expanding estimation sample: a test return can enter estimation only after its date has passed. Thus the return being predicted remains excluded from the information used for that prediction, even when earlier test returns have become available for re-estimation.
Code
training_end <- max(sp500$date)
test_prices <- tq_get(
"^GSPC", from = as.character(training_end + 1),
to = as.character(training_end + 60)
) |>
select(date, close) |>
rename(price = close) |>
arrange(date) |>
filter(date > training_end) |>
slice_head(n = 20)
stopifnot(
nrow(test_prices) == 20L, !anyDuplicated(test_prices$date),
all(is.finite(test_prices$price)), all(test_prices$price > 0),
all(diff(test_prices$date) > 0)
)
test_returns <- test_prices |>
mutate(
return = price / lag(price, default = tail(sp500$price, 1)) - 1,
horizon = row_number(), squared_return = return^2
)
tibble(
sample = c("Training returns", "Test returns"),
first_date = c(min(returns$date), min(test_returns$date)),
last_date = c(max(returns$date), max(test_returns$date)),
observations = c(nrow(returns), nrow(test_returns))
) |>
kable(caption = "Chronological separation of training and test returns.",
col.names = c("Sample", "First date", "Last date", "Returns"), row.names = FALSE) |>
kable_styling(full_width = FALSE, latex_options = "HOLD_position")| Sample | First date | Last date | Returns |
|---|---|---|---|
| Training returns | 2015-07-13 | 2020-07-09 | 1258 |
| Test returns | 2020-07-10 | 2020-08-06 | 20 |
The dates in the table show the resulting split. The test runs from 2020-07-10 through 2020-08-06, covering 20 trading observations rather than 20 calendar days. Its boundary follows the last price actually included in training, so the download cutoff cannot silently shift a return from one sample to the other.
We begin by asking how the arrival of these returns changes the model’s assessment of risk. Keep the 20 forecasts originally made at \(T\) unchanged, and alongside them calculate a second path that incorporates each test return once it has been observed. Denoting this updated path by \(\tilde v_{T+h}\) distinguishes it from the fixed-origin forecast \(\hat v_{T+h\mid T}\). For either fitted parameter vector, the update is
\[ \begin{aligned} \tilde v_{T+1}&=\hat v_{T+1\mid T},\\ \tilde v_{T+h}&=\hat\omega+\hat\alpha u_{T+h-1}^{2} +\hat\beta\tilde v_{T+h-1},\qquad h=2,\ldots,20. \end{aligned} \]
For the first test day, both paths use the final training return and variance, so their values coincide. For the second day, the updated path can use the first test return, whereas the original forecast must still use its expectation at \(T\). The same distinction continues through the window. Although we reconstruct the updated path after observing the data, every update uses returns only through the preceding day. It is therefore also the sequence of one-step forecasts that could have been issued as the test unfolded, without changing the coefficients. Comparing it with the fixed-origin path shows the effect of new information within the model; both paths remain model-based estimates of an unobserved variance.
This timing also explains why the two lines can separate substantially. For the second test day, subtracting the original forecast from the updated variance gives
\[ \tilde v_{T+2}-\hat v_{T+2\mid T} =\hat\alpha\left(u_{T+1}^2-\hat v_{T+1\mid T}\right). \]
The original forecast used the expected squared return for day \(T+1\). Once that return is observed, the update replaces the expectation by its realized square. A smaller squared return therefore lowers the next day’s updated variance relative to the original forecast; a larger one raises it. On subsequent days, the model also carries forward the variance changes caused by earlier surprises. A sequence of relatively quiet returns can therefore pull the updated path well below the original forecast even though the coefficients remain unchanged. Conversely, the fixed-origin path keeps moving toward its fitted long-run level because it receives no new information. This distinction applies separately to both estimators.
Code
test_garch_paths <- function(omega, alpha, beta, last_variance, model) {
h <- test_returns$horizon
first_variance <- omega + alpha * tail(u, 1)^2 + beta * last_variance
long_run <- omega / (1 - alpha - beta)
fixed_origin <- long_run + (alpha + beta)^(h - 1) * (first_variance - long_run)
filtered <- numeric(length(h))
filtered[1] <- first_variance
for (i in 2:length(h)) {
filtered[i] <- omega + alpha * test_returns$return[i - 1]^2 + beta * filtered[i - 1]
}
test_returns |>
mutate(model = model, fixed_origin = as.numeric(fixed_origin), filtered = filtered)
}
test_paths <- bind_rows(
test_garch_paths(omega_full, alpha_full, beta_full, tail(v_full, 1), "Full MLE"),
test_garch_paths(omega_vt, alpha_vt, beta_vt, tail(v_vt, 1), "Variance targeting")
) |>
mutate(model = factor(model, levels = c("Full MLE", "Variance targeting")))
stopifnot(isTRUE(all.equal(
test_paths$fixed_origin[test_paths$model == "Variance targeting"],
as.numeric(forecast_variance), tolerance = 1e-12
)))Code
test_paths |>
pivot_longer(c(fixed_origin, filtered), names_to = "path", values_to = "variance") |>
mutate(path = factor(path, levels = c("fixed_origin", "filtered"),
labels = c("Forecast at T", "Updated with observed returns"))) |>
ggplot(aes(date, 100 * sqrt(variance), color = path)) +
geom_line(linewidth = 0.8) +
facet_wrap(~ model) +
scale_color_manual(values = c("#CC7A00", "#1B4E8A")) +
scale_x_date(date_breaks = "1 week", date_labels = "%b %d") +
scale_y_continuous(labels = scales::label_percent(accuracy = 0.01, scale = 1)) +
labs(x = "Test date", y = "Daily volatility", color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom")
The orange line in each panel preserves what could be forecast at \(T\), while the blue line shows how that estimator responds to the returns as they arrive. For continuity with Section 3.9, both paths are shown as daily volatility percentages: we take the square root of each variance and multiply by 100. A value of 1% therefore has the same meaning in these figures. The shorter test window permits a closer view of the movement, so its axis limits differ from those of the full historical plot while the units remain the same. This display transformation leaves the underlying variance forecasts and all subsequent loss calculations unchanged. Within each panel, the divergence between the lines comes entirely from the information update; their similarity across panels reflects how close the two estimated models are.
By the last test day, the updated daily volatility is approximately 0.755% for Full MLE and 0.757% for variance targeting. The original forecasts for that day remain at 1.126% and 1.144%. The returns observed over the window thus lead both models to a lower end-of-period assessment of risk. This illustrates their reaction to new information, but evaluating forecast accuracy requires a separate observable outcome. Using the filtered path itself as the outcome would merely compare two outputs of the same fitted model.
The third policy adds parameter re-estimation to this daily information update. Before predicting day \(t\), we extend the estimation sample through \(t-1\), fit the model, and reconstruct its variance path with those newly estimated coefficients. Denote the coefficients by \(\hat\omega_{t-1}\), \(\hat\alpha_{t-1}\), and \(\hat\beta_{t-1}\), where the subscript identifies the last observation available for estimation. The resulting one-day forecast is
\[ \hat v_{t\mid t-1}^{\mathrm{refit}} =\hat\omega_{t-1}+\hat\alpha_{t-1}u_{t-1}^2 +\hat\beta_{t-1}\hat v_{t-1}^{\mathrm{refit}}. \]
Here \(\hat v_{t-1}^{\mathrm{refit}}\) is the final filtered variance obtained by running the recursion over the expanded sample with the new coefficients. Reconstructing this path keeps the variance state consistent with the model just estimated. For variance targeting, we also recalculate \(V_L\) from the expanded sample and derive the new intercept from that value. Full MLE continues to estimate all three coefficients jointly. Neither procedure uses the return on day \(t\) until after forecasting that day.
These policies answer different operational questions. The original twenty-day forecast describes what can be planned at the initial date, before any test returns arrive. Daily updating describes a risk assessment that can be revised as the market evolves while keeping a previously calibrated model. Daily re-estimation also allows the calibration to change, at the cost of repeating the optimization and introducing estimation variation. Comparing the first two policies measures the benefit of fresh information in this window; comparing the last two isolates the additional effect of re-estimating the model on an expanding sample.
Code
daily_refits <- list()
for (h in seq_len(nrow(test_returns))) {
past <- c(u, head(test_returns$return, h - 1L))
if (h == 1L) {
pf <- unname(theta_full)
pv <- c(unname(omega_vt), unname(alpha_vt), unname(beta_vt))
} else {
full_fit <- optim(c(4e-6, 0.2, 0.7), nll_garch, u = past,
method = "L-BFGS-B", lower = c(1e-12, 0, 0),
upper = c(Inf, 1 - 1e-6, 1 - 1e-6),
control = list(parscale = c(1e-6, 0.1, 0.1), factr = 1e7, pgtol = 1e-10))
target <- var(past)
vt_fit <- optim(c(0.2, 0.7), nll_garch_vt, u = past,
long_run_variance = target, method = "L-BFGS-B",
lower = c(0, 0), upper = c(1 - 1e-6, 1 - 1e-6))
# Retry an unsuccessful search using the previous fit and finer numerical steps.
if (vt_fit$convergence != 0) {
vt_fit <- optim(unname(pv[2:3]), nll_garch_vt, u = past,
long_run_variance = target, method = "L-BFGS-B",
lower = c(0, 0), upper = c(1 - 1e-6, 1 - 1e-6),
control = list(parscale = c(0.1, 0.1), ndeps = c(1e-5, 1e-5), maxit = 2000))
}
stopifnot(full_fit$convergence == 0, vt_fit$convergence == 0)
pf <- unname(full_fit$par)
pv <- c(target * (1 - sum(vt_fit$par)), unname(vt_fit$par))
stopifnot(nll_garch(pf, past) <= nll_garch(pv, past) + 1e-6)
}
for (model in c("Full MLE", "Variance targeting")) {
p <- if (model == "Full MLE") pf else pv
stopifnot(all(p >= 0), p[1] > 0, sum(p[2:3]) < 1)
last_v <- tail(garch_variance(past, p[1], p[2], p[3]), 1)
prediction <- p[1] + p[2] * tail(past, 1)^2 + p[3] * last_v
daily_refits[[length(daily_refits) + 1L]] <- tibble(
date = test_returns$date[h], model = model, forecast = prediction,
squared_return = test_returns$squared_return[h],
policy = "Daily re-estimation", training_n = length(past),
omega = p[1], alpha = p[2], beta = p[3])
}
}
daily_refits <- bind_rows(daily_refits)
forecast_comparison <- test_paths |>
select(date, model, squared_return, fixed_origin, filtered) |>
pivot_longer(c(fixed_origin, filtered), names_to = "policy", values_to = "forecast") |>
mutate(policy = ifelse(policy == "fixed_origin", "Original 20-day forecast", "Daily update, fixed parameters")) |>
bind_rows(daily_refits |> select(date, model, squared_return, policy, forecast)) |>
mutate(policy = factor(policy, levels = c("Original 20-day forecast",
"Daily update, fixed parameters", "Daily re-estimation")))
stopifnot(all(is.finite(forecast_comparison$forecast)), all(forecast_comparison$forecast > 0))
first_day_check <- forecast_comparison |>
filter(date == min(date)) |>
group_by(model) |>
summarise(gap = max(forecast) - min(forecast), .groups = "drop")
stopifnot(all(first_day_check$gap < 1e-12))The expression head(test_returns$return, h - 1L) is the chronological safeguard in this calculation: before test day \(h\), only the preceding \(h-1\) test returns can be appended to training. On the first test day, all three policies consequently give the same forecast for each estimator. Later re-estimations retain the original initialization rule and check convergence, parameter restrictions, and likelihood ordering. If the initial variance-targeting search fails to converge, the code retries from the previous day’s estimate with finer numerical steps; an unsuccessful retry stops the calculation. These numerical checks concern estimation, without using test losses to select parameter values.
To evaluate all three policies, we turn to the squared return \(u_{T+h}^2\). Under the zero-conditional-mean assumption, its conditional expectation equals the day’s latent variance, making it a useful observable proxy. Any one squared return can nevertheless lie far above or below that variance. A quiet realization is possible on a risky day, just as a large movement can occur on a day assigned relatively low risk. We must consequently choose a loss function suited to a noisy proxy and interpret averages over only 20 observations cautiously. The role of that choice is examined in Patton (2011).
For this second comparison, the next figure uses variance units because the observable outcome is a squared return. Both the forecast variance and the squared return are multiplied by 10,000 to express them in squared percentage points. For example, a daily volatility of 1% corresponds to a variance of 1 squared percentage point; a volatility of 2% corresponds to 4 squared percentage points. Each panel shows the variance forecast issued under each policy and plots the squared return actually observed on each test date. All six panels use the same observed returns, so differences in forecast performance arise from the forecast lines. The variation of the points around those lines also shows why a single squared return provides such an imprecise assessment of that day’s underlying risk.
Code
forecast_comparison |>
ggplot(aes(date)) +
geom_line(aes(y = 10000 * forecast, color = "Variance forecast"), linewidth = 0.8) +
geom_line(aes(y = 10000 * squared_return, color = "Observed squared return"), linewidth = 0.35) +
geom_point(aes(y = 10000 * squared_return, color = "Observed squared return"), size = 1.5) +
facet_grid(model ~ policy) +
scale_color_manual(values = c("Variance forecast" = "#CC7A00", "Observed squared return" = "#1B4E8A")) +
scale_x_date(date_breaks = "2 weeks", date_labels = "%b %d") +
labs(x = "Test date", y = "Daily variance (squared percentage points)", color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom", panel.spacing = grid::unit(0.7, "cm"),
strip.text.x = element_text(size = 10),
panel.border = element_rect(color = "grey65", fill = NA, linewidth = 0.5))
The squared returns vary much more abruptly than the forecast paths because each point records one realized daily movement, whereas each forecast describes an expected squared movement using the information available when that forecast was issued. The updated variances underlying the preceding volatility figure smooth these movements because the recursion combines the previous squared return with the previous variance. They also respond with a one-day lag: a return observed today enters tomorrow’s conditional variance. These differences in timing and meaning account for the different shapes of the figures. The loss calculations below aggregate the discrepancies in this second comparison on the same test dates. Each forecast is aligned with the return it was issued to predict; the observed series is connected by a thin line to make its chronological sequence easier to follow.
For each estimator and policy, let \(f_h\) be its variance forecast for test day \(T+h\), issued using the information allowed by that policy, and let \(q_h=u_{T+h}^2\) be the subsequently observed proxy. We score their discrepancy using two average losses, QLIKE and mean squared error:
\[ \begin{aligned} \mathrm{QLIKE}&=\frac{1}{20}\sum_{h=1}^{20} \left[\log(f_h)+\frac{q_h}{f_h}\right],\\ \mathrm{MSE}&=\frac{1}{20}\sum_{h=1}^{20}(q_h-f_h)^2. \end{aligned} \]
The form of QLIKE should look familiar from the likelihood calculation. Assigning a large variance raises \(\log(f_h)\), while assigning too little variance to a large observed squared return raises \(q_h/f_h\). Their sum balances the two costs, now judging forecasts already issued rather than choosing parameters to fit the sample. Terms depending only on the proxy can be omitted because they would add the same amount to each model’s loss on the same dates. This version also remains defined when a return is zero. Its value can be negative, so performance is read by comparing which loss is lower, regardless of the sign.
MSE supplies a more direct distance measure by squaring the gap between the forecast variance and the squared return. It gives substantial weight to large discrepancies and has units equal to the fourth power of decimal returns. To make the small values readable, the table multiplies MSE by \(10^8\), equivalent to calculating it from variances in squared percentage points. Under a conditionally unbiased proxy and the required moment conditions, both criteria have useful expected-ranking properties (Patton 2011). Their use makes the comparison appropriate for this target, while the noise in a short test sample remains.
Code
evaluation_summary <- test_paths |>
group_by(model) |>
summarise(
observations = n(),
qlike = mean(log(fixed_origin) + squared_return / fixed_origin),
mse = mean((squared_return - fixed_origin)^2),
mean_forecast = mean(fixed_origin),
mean_proxy = mean(squared_return),
.groups = "drop"
)
policy_scores <- forecast_comparison |>
group_by(model, policy) |>
summarise(observations = n(),
qlike = mean(log(forecast) + squared_return / forecast),
mse = mean((squared_return - forecast)^2), .groups = "drop")
policy_scores |>
transmute(model, policy, qlike, mse_scaled = 1e8 * mse) |>
kable(caption = "Forecast losses across the three policies: the same twenty test dates for each estimator and policy.",
col.names = c("Estimator", "Forecasting policy", "QLIKE", "\\(\\mathrm{MSE}\\times10^8\\)"),
digits = c(0, 0, 6, 6), escape = FALSE, row.names = FALSE) |>
kable_styling(full_width = FALSE, latex_options = "HOLD_position")| Estimator | Forecasting policy | QLIKE | \(\mathrm{MSE}\times10^8\) |
|---|---|---|---|
| Full MLE | Original 20-day forecast | -8.494497 | 0.626010 |
| Full MLE | Daily update, fixed parameters | -8.640791 | 0.310189 |
| Full MLE | Daily re-estimation | -8.640958 | 0.309794 |
| Variance targeting | Original 20-day forecast | -8.483898 | 0.659523 |
| Variance targeting | Daily update, fixed parameters | -8.639352 | 0.312896 |
| Variance targeting | Daily re-estimation | -8.639243 | 0.312886 |
Code
qlike_best <- as.character(evaluation_summary$model[which.min(evaluation_summary$qlike)])
mse_best <- as.character(evaluation_summary$model[which.min(evaluation_summary$mse)])
loss_reading <- if (qlike_best == mse_best) {
paste(qlike_best, "has the lower loss under both QLIKE and MSE")
} else {
paste(qlike_best, "has the lower QLIKE, while", mse_best, "has the lower MSE")
}
update_gains <- policy_scores |>
select(model, policy, mse) |>
pivot_wider(names_from = policy, values_from = mse) |>
mutate(reduction = 100 * (1 - `Daily update, fixed parameters` / `Original 20-day forecast`))For the original twenty-day forecasts, Full MLE has the lower loss under both QLIKE and MSE. The QLIKE values are -8.494497 for Full MLE and -8.483898 for variance targeting; the corresponding scaled MSE values are 0.626010 and 0.659523. To understand this ordering, compare the risk forecasts with the size of the observed movements. The average squared return is 0.623 in squared-percentage-point units, compared with average variance forecasts of 1.227 and 1.251. Both models therefore assign more variance on average than the observed proxy in this window. Full MLE projects slightly less variance, and its forecasts yield smaller discrepancies under both loss functions. That is the empirical result for these dates; it gives no basis for claiming a statistically significant or lasting advantage over variance targeting.
The larger improvement in the table comes from updating the forecasts each day. Relative to the original twenty-day path, daily updating with fixed parameters reduces MSE by 50.4% for full MLE and 52.6% for variance targeting. QLIKE also falls for both estimators. In the middle column of the figure, the forecasts move downward as the quieter returns arrive, bringing them closer to the subsequent squared returns on average. Some large discrepancies remain because the forecast describes conditional variance, while each squared return records just one draw from that day’s return distribution. The improvement is measured over all twenty dates; it does not imply that every individual forecast becomes more accurate.
Daily re-estimation changes the picture very little beyond that update. The last two columns are almost indistinguishable, and their losses are very close. Full MLE improves slightly under both criteria when re-estimated; variance targeting has a marginally lower MSE but a marginally higher QLIKE. There is consequently little evidence here of an additional practical benefit from re-estimating every day. The original estimation sample already contains 1258 returns. Before the final test forecast, it has grown by only 19 returns, or 1.5% of its original size. The twentieth test return becomes available only after the last forecast has been issued. These few additions change the fitted coefficients little in this sample, whereas the daily variance recursion responds directly to each latest squared return. This explains why most of the adaptation comes from updating the variance state. A shorter estimation window or a different market episode could yield a more substantial effect from re-estimation.
The comparison therefore separates two decisions: whether forecasts can incorporate new observations, and whether the model is recalibrated when they do. Its interpretation must respect the information advantage of the daily policies. The original policy makes forecasts at horizons one through twenty from one date; the daily policies make twenty successive one-step forecasts. Lower daily-policy losses describe the benefit of that forecasting practice in this window, rather than a contest between procedures supplied with identical information. Within each policy, full MLE and variance targeting use the same information dates and can be compared directly.
A broader assessment would repeat this design over additional test windows and market conditions. Twenty daily outcomes provide limited evidence, especially with a noisy variance proxy. The favorable residual checks, the small in-sample likelihood difference, and the test losses answer related questions about remaining dependence, fit, and predictive accuracy. Together they give a fuller assessment than either visual agreement or a single loss ranking alone.
3.12 Common implementation mistakes
Several implementation errors can produce plausible-looking results while breaking the link between the model and the forecast. An especially easy one is to use \(u_t^2\) in the calculation of \(v_t\). That would let the model use the very return whose variance it is supposed to predict. Working through the first rows of the daily calculation table helps expose the mistake: the update must use \(u_{t-1}^2\) and \(v_{t-1}\), and the initialized value must appear before the first recursive update. The price-row numbers also need to remain distinct from the indices of the return vector after its first missing observation has been removed.
A similarly unobtrusive error is to mix returns expressed as percentages with variances calculated from decimal returns. Multiplying a return by 100 multiplies its variance by 10,000, so the mismatch changes the likelihood as well as the apparent magnitude of risk. Keeping decimal units throughout estimation and converting only for presentation makes it easier to check that the squared returns and conditional variances are on the same scale.
Even with correct timing and units, the optimizer’s output needs scrutiny. The parameter vector must satisfy \(\omega>0\), \(\alpha\ge0\), \(\beta\ge0\), and \(\alpha+\beta<1\) for the positive, finite-long-run-variance model used here. The coordinate bounds do not enforce the final joint restriction, which is why it is also checked inside the objective. Numerical convergence then needs to be interpreted alongside the attained likelihood. Trying several feasible starting vectors should give similar estimates in a stable problem, and scaling the parameters helps avoid search difficulties caused by their different magnitudes. An additional check is available in this application: the unrestricted likelihood should be at least as high as the likelihood at any admissible variance-targeting solution. A successful convergence code alone would not reveal a violation of that ordering.
Once a fit has passed these checks, a comparison can still be misleading if the reported statistic uses the wrong likelihood scale or observation set. The proportional score is sufficient for estimating parameters, but gaussian_loglik() supplies the complete value required by AIC and BIC. Both models must be evaluated on the same valid dates, with the treatment of the sample-based long-run variance made explicit in the parameter count. Otherwise an apparent preference can arise from the reporting convention rather than a substantive difference between models.
The interpretation needs the same care as the calculation. Information criteria describe a penalized fit and should be considered together with parameter plausibility and residual behavior. A claim about forecasting performance additionally requires later observations that were excluded from fitting. In that assessment, keeping the filtered GARCH path separate from the observed squared-return proxy prevents a comparison between model outputs from being mistaken for validation against data. These are the connections that allow the numerical results to answer the financial question posed at the start of the chapter.
3.13 Conclusion and extensions
We began with the observation that the size of daily S&P 500 movements changes over time. GARCH gives that observation a quantitative form by allowing the latest squared return to revise a variance estimate that also retains part of its preceding value. Following the recursion through the daily tables showed how this rule connects the data to the likelihood. Maximizing that likelihood then provided the coefficients needed to reconstruct historical volatility and project risk beyond the end of the sample.
The two estimators illustrate how a restriction can simplify this process without greatly changing its empirical description. Full MLE chooses \((\omega,\alpha,\beta)\) jointly, while variance targeting estimates the two weights and derives the intercept from a sample-based long-run variance. Both optimizations converge, and the freely estimated model attains the slightly higher Gaussian likelihood that its larger feasible set permits. Their persistent variance dynamics and similar standardized-residual diagnostics support broadly comparable readings of historical risk. The AIC/BIC preference depends on how the sample-based input is counted: variance targeting is favored under the conditional count \(k=2\), while counting \(V_L\) gives \(k=3\) and a narrow preference for Full MLE. Given the small difference, this comparison provides little evidence of a substantial gap in fit.
Forecasting carries the same recursion forward with a change in information. Historical variance responds to observed shocks; a forecast averages over shocks still to come, which explains its smoother approach to the long-run level. The simulated band around the variance-targeting forecast shows how future conditional volatility can nevertheless evolve along quite different paths with the fitted coefficients held fixed. The target throughout is the variability of future returns, so a useful risk forecast need not imply an ability to predict their direction.
The subsequent observations, from 2020-07-10 through 2020-08-06, let us distinguish that original prediction from the updates made possible by new returns. Both models revised their end-of-window risk assessment downward as the test unfolded. Comparing three forecasting policies showed that daily updating with fixed coefficients substantially reduced both loss criteria relative to the original twenty-day forecast. Re-estimating each day added little in this window, where the expanding sample gained at most 19 returns before the final forecast, alongside 1258 original training returns. Full MLE and variance targeting produced very similar paths under each policy. These findings illustrate the value of new information for daily risk assessment while keeping the distinction between one-step and multi-step forecasts explicit. They concern one short window with a noisy variance proxy; a general claim of forecasting superiority would require evaluation over additional windows and market conditions. Read together, the estimation results, residual checks, and test losses show both what this simple model can explain and the limits of the evidence used to assess it.
Further developments can be motivated by the assumptions that remain. An asymmetric response to gains and losses would lead toward EGARCH or GJR-GARCH, heavier conditional tails toward Student-t innovations, and joint risk across assets toward multivariate models such as DCC-GARCH. These extensions are discussed in Francq and Zakoian (2019), Tsay (2010), and McNeil et al. (2015). Each changes part of the model, but the route from an economic question to an explicit variance equation, an estimable criterion, and an honest forecast evaluation remains the basis for understanding and using the result.