source("common/functions.r", chdir = TRUE)
library(tsgarch)
library(xts)20 When things go wrong
The estimation in the previous chapter went through without any problems, but we are not always so lucky. GARCH estimation maximises a non-linear likelihood, and when it fails it usually fails with an incomprehensible error message. The cause is one of four things — the data, the specification, the optimiser, or our own code — and the message rarely says which.
These problems can arise because of the trade-off between safe code and speed. If we want to ensure that nothing goes wrong, we need many checks and robust algorithms, which are slow, even very slow. Consequently, we opt for algorithms that are fast and usually work but sometimes fail.
20.1 Data and libraries
The GARCH packages again. What differs here is the data we feed them.
import numpy as np
import pandas as pd
import sys
sys.path.insert(0, 'common')
from functions import ProcessRawData
from arch import arch_modelusing CSV, DataFrames, Dates, Statistics
using ARCHModels
using Random
include("common/functions.jl");data = ProcessRawData()
dates = as.Date(as.character(data$sp500$date), format = "%Y%m%d")
y = xts(data$sp500$y, order.by = dates)
y = na.omit(y)
y = y - mean(y)data = ProcessRawData()
y = data['sp500']['y'].dropna().values
y = y - np.mean(y)data = ProcessRawData();
y = collect(skipmissing(data["sp500"].y));
y = y .- mean(y);20.2 Numerical precision problems
Most calculations are done with what is known as double precision floating point, called float64, a number that uses 64 bits in computer memory. A floating point number such as \(0.0234=2.34\times 10^{-2}\) is written in code as 2.34e-2.
The largest representable number is about \(1.8\times10^{308}\) (approximately \(2^{1024}\)), and the smallest normal positive number is \(2^{-1022}\) (subnormal numbers reach as low as \(2^{-1074}\)). A float64 is laid out in three fields — one sign bit, eleven exponent bits and fifty-two stored fraction bits. The exponent fixes the magnitude and the fraction carries the precision. Because normal numbers have an implicit leading one that is not stored, those fifty-two stored bits deliver fifty-three bits of significand precision.
Floating-point arithmetic loses accuracy when an operation combines values of very different magnitudes, because the smaller one is rounded away against the larger. Summing the same numbers in a different order therefore gives a different answer, and adding the smallest first loses least. Taking the logarithm of a product can overflow or underflow where summing the logarithms does not, as Chapter 26 shows.
When returns are very small or very large, squaring them amplifies the differences. The covariance matrix of the parameter estimates is obtained from the inverse of the negative Hessian of the log-likelihood, the matrix of second derivatives at the maximum. It is tricky to estimate this Hessian, especially if we want it done quickly. That is why the scale of the data matters here. Badly scaled data leave the Hessian ill-conditioned, so the standard errors can be unreliable or fail outright even when the point estimates themselves converge.
20.3 A failure to start with
Nothing has gone wrong yet, so let us make it go wrong. Daily returns are of the order of 0.01, and squaring them takes us to 0.0001. We fit a Gaussian GARCH(1,1) to the first 2,000 returns at that raw scale, with no mean equation, and see what each package makes of it.
window = y[1:2000]
spec_raw = garch_modelspec(window, model = "garch", order = c(1, 1), distribution = "norm")
fit_raw = estimate(spec_raw)
round(coef(fit_raw), 8) omega alpha1 beta1
0.00000023 0.03202469 0.96556317import warnings
window_py = y[:2000]
with warnings.catch_warnings(record = True) as caught:
warnings.simplefilter("always")
model_raw = arch_model(window_py, vol = 'Garch', p = 1, q = 1, mean = 'Zero', dist = 'normal')
fit_raw = model_raw.fit(disp = 'off')
for w in caught:
print(f"{w.category.__name__}: {str(w.message).splitlines()[0]}")
print(fit_raw.params.round(8))DataScaleWarning: y is poorly scaled, which may affect convergence of the optimizer when
omega 0.000001
alpha[1] 0.050000
beta[1] 0.930000
Name: params, dtype: float64window_jl = y[1:2000];
fit_raw = fit(GARCH{1,1}, window_jl; meanspec=NoIntercept, dist=StdNormal);
round.(coef(fit_raw), digits=8)R and Julia estimate this without complaint. Python does not. arch warns that the series is poorly scaled and then returns \(\ARCHcoeff=0.05\) and \(\GARCHcoeff=0.93\), which are not estimates at all — they are the values the optimiser started from, and it never moved away from them. Note that the fit reports success. Nothing raised an exception, and only the warning and the suspiciously round numbers reveal that the answer is worthless.
An estimation that fails loudly is easy to deal with. One that quietly returns its starting values is not.
The rest of this section fixes it.
20.4 Rescaling to fix problems
Rescaling the data resolves numerical precision problems by putting the numbers that enter the likelihood on a similar magnitude. Two adjustments are available, and they do different things. De-meaning shifts the location of the data and leaves its spread alone. Dividing by the standard deviation changes the magnitude and leaves the location alone.
20.4.1 De-meaning
The simplest adjustment is to remove the mean from the estimation window rather than the entire dataset. When we imported the data, we removed the mean from all the observations. If, however, we de-mean a specific estimation window, the calculation is more likely to succeed.
window_data = y[720:2719]
window_data = window_data - as.numeric(mean(window_data))
spec = garch_modelspec(window_data, model = "garch", order = c(1, 1), distribution = "std")
Res = estimate(spec)
coef(Res) omega alpha1 beta1 shape
3.732500e-07 5.471533e-02 9.437558e-01 5.773846e+00window_data = y[719:2719] # Python is 0-indexed
window_data = window_data - np.mean(window_data)
model = arch_model(window_data, vol='Garch', p=1, q=1, dist='t', mean='Zero')
res = model.fit(disp='off')
print(res.params)/Users/jond/.cache/uv/archive-v0/ir4pFzJQ83o6Q8vWzfwj6/lib/python3.14/site-packages/arch/univariate/base.py:694: DataScaleWarning: y is poorly scaled, which may affect convergence of the optimizer when
estimating the model parameters. The scale of y is 9.011e-05. Parameter
estimation work better when this value is between 1 and 1000. The recommended
rescaling is 100 * y.
This warning can be disabled by either rescaling y before initializing the
model or by setting rescale=False.
self._check_scale(resids)omega 0.000002
alpha[1] 0.102114
beta[1] 0.882570
nu 6.534431
Name: params, dtype: float64using ARCHModels, Statistics
window_data = y[720:2719];
window_data = window_data .- mean(window_data);
res_jl = fit(GARCH{1,1}, window_data; dist=StdT, meanspec=NoIntercept);
coef(res_jl)De-meaning alone does not resolve the scale problem. The returns remain tiny in absolute terms, and the Python fit still warns about poorly scaled data. The next step addresses the scale itself.
20.4.2 Normalising variance
Normalising the variance to be one also works, but then we need to rescale \(\GARCHconst\) back since \[
\Vol^2 = \frac{\GARCHconst}{1-\ARCHcoeff-\GARCHcoeff}
\] and so if we do data=data/sd(data) we recover the original scale with \[
\hat{\GARCHconst}\times \Var(\CompoundReturns)
\]
window_data = y[720:2719]
scale_factor = as.numeric(sd(window_data))
scaled_data = window_data / scale_factor
spec = garch_modelspec(scaled_data, model = "garch", order = c(1, 1), distribution = "std")
Res = estimate(spec)
coef(Res)
# Rescale omega back
omega_rescaled = coef(Res)["omega"] * scale_factor^2
omega_rescaled omega alpha1 beta1 shape
0.003959947 0.053036774 0.945527321 5.776259351
omega
3.57022e-07window_data = y[719:2719]
scale_factor = np.std(window_data, ddof=1)
scaled_data = window_data / scale_factor
model = arch_model(scaled_data, vol='Garch', p=1, q=1, dist='t', mean='Zero')
res = model.fit(disp='off')
print(res.params)
# Rescale omega back
omega_rescaled = res.params['omega'] * scale_factor**2
print(f"Rescaled omega: {omega_rescaled}")omega 0.003661
alpha[1] 0.046374
beta[1] 0.951598
nu 6.009776
Name: params, dtype: float64
Rescaled omega: 3.300653209916117e-07using ARCHModels, Statistics
window_data = y[720:2719];
scale_factor_jl = std(window_data);
scaled_data = window_data ./ scale_factor_jl;
res_jl = fit(GARCH{1,1}, scaled_data; dist=StdT, meanspec=NoIntercept);
coef(res_jl)
# Rescale omega back
omega_rescaled_jl = coef(res_jl)[1] * scale_factor_jl^2It is usually best to do both, de-meaning and normalising the data.
20.4.3 Why does rescaling work?
De-meaning and normalisation work because they make most of the numbers close to one. When we square them, the gap between the big and the small numbers does not become excessively large, so we do not get the loss of precision we discussed above.
20.5 Detecting failures in automated procedures
Suppose rescaling does not work and estimation still fails. If doing one estimation, it is easy to spot the problem and do something else. But what if we are doing backtesting and want to generate a large number of one-day risk forecasts? Then, it becomes annoying to have the code fail.
In R, the try() function catches errors so the code can continue even when individual operations fail.
# This will cause an error
log("string")Error: non-numeric argument to mathematical functionWith try():
res1 = try(log(1))
res2 = try(log("string"), silent = TRUE)Both ran without stopping the script. We can check the results:
res1
class(res1)[1] 0
[1] "numeric"res2
class(res2)[1] "Error in log(\"string\") : non-numeric argument to mathematical function\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
<simpleError in log("string"): non-numeric argument to mathematical function>
[1] "try-error"Detect failures by:
class(res2) == "try-error"[1] TRUEAfter running try() we check the output, and if its class is "try-error", we know the calculation failed and should deal with that.
For GARCH estimation in a loop:
for (i in 1:n_windows) {
res = try(estimate(spec), silent = TRUE)
if (class(res) == "try-error") {
# Handle failure: use previous estimate, skip, or try alternative
next
}
# Process successful result
}In Python, we use try/except blocks to catch exceptions:
import math
# This will cause an error
try:
result = math.log("string")
except (TypeError, ValueError) as e:
print(f"Error caught: {e}")
result = NoneError caught: must be real number, not strFor GARCH estimation:
for i in range(n_windows):
try:
model = arch_model(window_data, vol='Garch', p=1, q=1)
res = model.fit(disp='off')
# Process successful result
except Exception as e:
print(f"Estimation failed for window {i}: {e}")
# Handle failure: use previous estimate, skip, or try alternative
continueMore specific exception types are also available:
from arch.univariate.base import ConvergenceWarning
import warnings
with warnings.catch_warnings():
warnings.filterwarnings('error', category=ConvergenceWarning)
try:
res = model.fit(disp='off')
except ConvergenceWarning:
print("Model did not converge")In Julia, we use try/catch blocks:
# This will cause an error
result = try
log("string")
catch e
println("Error caught: ", e)
nothing
endError caught: MethodError(log, ("string",), 0x00000000000097f7)For GARCH estimation:
for i in 1:n_windows
result = try
fit(GARCH{1,1}, window_data)
catch e
println("Estimation failed for window $i: $e")
nothing
end
if result === nothing
# Handle failure: use previous estimate, skip, or try alternative
continue
end
# Process successful result
end20.6 A checklist for the other failures
Scale is the failure this chapter works through, because it is the one that bites most often. Several others turn up in practice, and the remedies belong in one place.
The parameters are constrained. We need \(\GARCHconst > 0\), \(\ARCHcoeff_i \geq 0\) and \(\GARCHcoeff_j \geq 0\) for the variance to stay positive, and often \(\sum_i \ARCHcoeff_i + \sum_j \GARCHcoeff_j < 1\) so that the unconditional variance is finite. That last condition is covariance stationarity, though the slides and much of the literature call it stationarity for short. When the optimiser is pushed onto one of these boundaries the estimate can be unreliable even if the fit reports success, and variance targeting, tighter bounds or a simpler specification usually help. A binding stationarity constraint is often a symptom of a structural break rather than a numerical problem.
Starting values matter, as the Python failure above showed in the worst way. If the optimiser is landing on local maxima, fit a simpler model first and use its estimates to initialise the larger one, or start it from several plausible points and keep the best likelihood.
Sometimes the model is simply wrong for the data. Begin with a normal GARCH before adding complexity, read the residual diagnostics from the simpler fit, and use likelihood ratio tests to compare nested models and information criteria such as AIC and BIC for models that are not nested. If asymmetry or long memory is the problem, GJR-GARCH and FIGARCH are the usual next steps, and a persistent structural break may call for a regime-switching model.
Data problems are the last category, and the one where caution is most needed. Check for recording errors such as decimal point shifts and sign errors, confirm that the sampling frequency is what we think it is, and decide deliberately whether to interpolate missing values or exclude them.
Be careful with outliers. In a volatility model the extreme returns are the signal, not noise, and they are precisely what the model is fitted to capture. Do not alter an observation without evidence that it is a recording error. Winsorising or trimming can be run as a documented sensitivity check, to see how much the conclusions depend on a handful of days, but not as a default cleaning step.
20.7 Solver sensitivity
How much do the results depend on the optimiser rather than the data? We fit the same GARCH(1,1) model in all three languages and vary the solver settings within each one.
To make the comparison mean anything, every tab below fits the same model to the same data — the first 2,000 returns, divided by their standard deviation, with no mean equation and normal innovations. Only the solver settings change.
tsgarch passes its control list straight to nloptr, so the entries have to be ones nloptr recognises. The safe way to build one is to start from nloptr_fast_options() and change a single field. Note that estimate() writes the fitted values back into the specification object, so each fit needs a fresh spec or it starts where the previous one stopped.
scaled = y[1:2000] / as.numeric(sd(y[1:2000]))
fresh_spec = function() {
garch_modelspec(scaled, model = "garch", order = c(1, 1), distribution = "norm")
}
ctrl = nloptr_fast_options()
cat("algorithm:", ctrl$algorithm, " maxeval:", ctrl$maxeval, " xtol_rel:", ctrl$xtol_rel, "\n")
res1 = estimate(fresh_spec(), control = nloptr_fast_options())
ctrl_short = nloptr_fast_options()
ctrl_short$maxeval = 5
res2 = estimate(fresh_spec(), control = ctrl_short)
ctrl_loose = nloptr_fast_options()
ctrl_loose$xtol_rel = 1e-2
res3 = estimate(fresh_spec(), control = ctrl_loose)
round(rbind(default = coef(res1), maxeval_5 = coef(res2), xtol_1e2 = coef(res3)), 6)Warning messages:
1: In sqrt(diag(cov)) : NaNs produced
2: In sqrt(diag(object$cov.fixed)) : NaNs producedalgorithm: NLOPT_LD_SLSQP maxeval: 1000 xtol_rel: 1e-14
omega alpha1 beta1
default 0.003657 0.032025 0.965563
maxeval_5 0.018327 0.061887 0.921144
xtol_1e2 0.003861 0.030701 0.966461The truncated run stops well short of the optimum. tsgarch records the optimality conditions, and they show it:
c(default = res1$conditions$kkt2, maxeval_5 = res2$conditions$kkt2) default maxeval_5
TRUE FALSEarch uses scipy’s SLSQP and does not let us choose a different optimiser. The options argument is forwarded to SLSQP, which accepts entries such as maxiter and ftol, and passing anything else through it is silently ignored. What we can vary is tol, options and starting_values.
scaled_py = y[:2000] / np.std(y[:2000], ddof = 1)
base = arch_model(scaled_py, vol = 'Garch', p = 1, q = 1, mean = 'Zero', dist = 'normal')
res1 = base.fit(disp = 'off')
res2 = base.fit(disp = 'off', options = {'maxiter': 3})
res3 = base.fit(disp = 'off', starting_values = np.array([0.2, 0.2, 0.6]))
for label, r in [("default", res1), ("maxiter=3", res2), ("start values", res3)]:
print(f"{label:14s} {np.round(r.params.values, 6)} loglik {r.loglikelihood:.4f} flag {r.convergence_flag}")/Users/jond/github/notebook/_prebake/Volatility/univariate.vol.problems/prebake.py:87: ConvergenceWarning: The optimizer returned code 9. The message is:
Iteration limit reached
See scipy.optimize.fmin_slsqp for code meaning.
res2 = base.fit(disp = 'off', options = {'maxiter': 3})default [0.004038 0.032841 0.964215] loglik -2690.8134 flag 0
maxiter=3 [0.008861 0.042972 0.948299] loglik -2692.2256 flag 9
start values [0.004038 0.032842 0.964214] loglik -2690.8134 flag 0A non-zero convergence_flag is the warning sign. Code 9 means the iteration limit was reached.
ARCHModels forwards its algorithm keyword to Optim, so here we can genuinely change the optimiser, which is the one thing the Python tab cannot do. The numeric Optim controls are not reachable this way — passing iterations or g_tol to fit raises a MethodError — so we vary the algorithm instead.
using Optim
scaled_jl = y[1:2000] ./ std(y[1:2000]);
res1 = fit(GARCH{1,1}, scaled_jl; meanspec=NoIntercept, dist=StdNormal);
res2 = fit(GARCH{1,1}, scaled_jl; meanspec=NoIntercept, dist=StdNormal, algorithm=LBFGS());
res3 = fit(GARCH{1,1}, scaled_jl; meanspec=NoIntercept, dist=StdNormal, algorithm=NelderMead());
for (label, r) in [("BFGS", res1), ("LBFGS", res2), ("NelderMead", res3)]
println(rpad(label, 12), round.(coef(r), digits=6), " loglik ", round(loglikelihood(r), digits=4))
endBFGS [0.003659, 0.965567, 0.032013] loglik -2691.6829
LBFGS [0.003659, 0.965567, 0.032013] loglik -2691.6829
NelderMead [0.05, 0.9, 0.05] loglik -2718.276ARCHModels reports coefficients in the order \(\GARCHconst\), \(\GARCHcoeff\), \(\ARCHcoeff\).
Three things come out of this.
Within a language, a solver setting that is tight enough will stop the optimiser before it arrives. In R a maxeval of five leaves the estimates a long way from the optimum and the second-order optimality condition fails. In Python maxiter of three returns code 9. Neither raises an error, so the only way to know is to read the convergence information the package provides rather than the coefficients alone.
Different starting values, by contrast, mostly do not matter here. Python starting from \((0.2, 0.2, 0.6)\) lands on the same optimum as the default, which tells us the likelihood is well behaved once the data are scaled. That is a property of this problem, not a general guarantee.
Across languages the remaining differences are not caused by the optimiser. Both tsgarch and arch default to SLSQP, and ARCHModels uses BFGS from Optim, but BFGS and LBFGS give the same answer here while NelderMead stalls at its starting values. What separates the packages is how they initialise the variance recursion and where they stop, which is also why their log-likelihoods are not directly comparable.