source("common/functions.r", chdir = TRUE)
library(nloptr)21 Manual estimation
Specialised libraries made estimation easy in Chapter 18. They also hid the machinery. This chapter opens that box.
Two pieces come out of it, the optimisation algorithm and the likelihood function. Knowing how each works is what lets us diagnose an estimation that fails and adapt a model the packages do not cover, which in quantitative finance is most of them.
21.1 Data and libraries
An optimiser in each language, since we are no longer using the packages’ own estimation routines.
import numpy as np
from scipy.optimize import minimize
from scipy.special import gammaln
import sys
sys.path.insert(0, 'common')
from functions import ProcessRawDatausing CSV, DataFrames, Dates, Statistics
using Optim
using SpecialFunctions
using Printf
include("common/functions.jl");data = ProcessRawData()
y_vec = na.omit(data$sp500$y)
y_vec = y_vec - mean(y_vec)data = ProcessRawData()
y = data['sp500']['y'].dropna().values
y = y - np.mean(y)data_jl = ProcessRawData();
y_jl = collect(skipmissing(data_jl["sp500"].y));
y_jl = y_jl .- mean(y_jl);21.2 Optimisers
Maximum likelihood estimation searches for the parameter values that maximise the likelihood function. The algorithm that conducts that search is called an optimiser.
21.2.1 How optimisers search
A classical algorithm is Newton-Raphson. It starts from an initial guess for the parameters, uses the gradient of the log-likelihood to find the direction of improvement and the curvature, the second derivative, to size each step, and repeats until the parameters stop changing.
The curvature at the optimum also carries information about precision, since a sharply peaked likelihood pins the parameters down tightly and a flat one leaves them uncertain.
Problems are rare for a small model such as GARCH(1,1), but grow more likely as the number of parameters increases, and the search can settle on a local maximum rather than the global one.
21.2.2 Optimisers minimise
Optimisers minimise functions by default. Since we need to maximise the likelihood, we return the negative log-likelihood.
21.2.3 Ensuring positive parameters
The parameters that go into the likelihood function must stay non-negative, with \(\GARCHconst\) strictly positive so that the variance recursion cannot collapse to zero (the degrees of freedom \(\DOF\) carries its own constraint, \(\DOF > 2\), once the GARCH-t likelihood arrives). There are at least two ways to ensure non-negativity:
- Take the absolute value of the parameters inside the likelihood function
- Impose constraints when calling the optimising function (bounds)
Only R’s implementation sets explicit bounds, through nloptr. Python and Julia use unbounded Nelder-Mead and rely on the absolute-value reflection inside the likelihood, together with a large penalty value whenever the computed \(\Vol_t^2 \leq 0\).
The reflection has a cost. abs() puts a kink in the objective at zero and gives every sign-flipped parameter vector the same likelihood, so the optimiser searches a space of mirrored duplicate optima. That is harmless for a derivative-free method such as Nelder-Mead. R keeps the reflection alongside its bounds as a defensive redundancy, though bounded Nelder-Mead never proposes a point outside lb = c(0, 0, 0), so abs() never has anything to correct there.
R provides the built-in optim() function, which uses Nelder-Mead by default. For more advanced optimisation, the nloptr package provides access to the NLopt library with many algorithms.
# Built-in optim with Nelder-Mead
res = optim(c(0.000001, 0.1, 0.85), fn = likelihood_function)
# NLopt with bounds
opts = list("algorithm" = "NLOPT_LN_NELDERMEAD", "xtol_rel" = 1.0e-8, maxeval = 1e4)
res = nloptr(x0 = c(0.000001, 0.1, 0.85),
eval_f = likelihood_function,
lb = c(0, 0, 0),
ub = c(1, 1, 1),
opts = opts)Python’s scipy.optimize module provides various optimisation algorithms. The minimize function is the main interface.
# Nelder-Mead (no bounds)
res = minimize(likelihood_function, x0=[0.000001, 0.1, 0.85], method='Nelder-Mead')
# L-BFGS-B (with bounds)
res = minimize(likelihood_function, x0=[0.000001, 0.1, 0.85], method='L-BFGS-B',
bounds=[(1e-10, 1), (1e-10, 1), (1e-10, 1)])Julia’s Optim.jl package provides optimisation algorithms. The optimize function is the main interface.
# Nelder-Mead (no bounds)
res = optimize(likelihood_function, [0.000001, 0.1, 0.85], NelderMead())
# L-BFGS with box constraints
res = optimize(likelihood_function, [0.0, 0.0, 0.0], [1.0, 1.0, 1.0],
[0.000001, 0.1, 0.85], Fminbox(LBFGS()))21.3 Coding the likelihood
Coding the likelihood by hand is what lets us change the model — a custom distribution or an extra term in the variance equation needs only a new likelihood function.
21.3.1 Likelihood issues
- Value of \(\Vol_1^2\): We need to assign a value to \(\Vol_1^2\), the initial value of the conditional variance. In general, we choose the unconditional variance of the data for this.
- Initial values of the parameters: Optimisation works iteratively. We assign initial values for the parameters, making sure they satisfy the function’s restrictions and the parameter bounds. This will have a small effect on the parameter values but can have a large effect on computing time.
We also need to pick an optimisation method and algorithm.
21.3.2 Normal GARCH likelihood
The log-likelihood function for a normal GARCH(1,1) model is:
\[ \begin{split} \log\lik &= \underbrace{-\frac{\SampleSize-1}{2}\log(2\pi)}_{\text{Constant}} \\ &\quad - \frac{1}{2}\sum_{t=2}^{\SampleSize}\left(\log(\GARCHconst + \ARCHcoeff \CompoundReturns^2_{t-1} + \GARCHcoeff \Vol^2_{t-1}) + \frac{\CompoundReturns^2_t}{\GARCHconst + \ARCHcoeff \CompoundReturns^2_{t-1} + \GARCHcoeff \Vol^2_{t-1}}\right) \end{split} \]
21.3.3 GARCH-t likelihood
The Student-t density is given by \[ \frac{\Gamma\left(\frac{\DOF+1}{2}\right)}{\sqrt{\DOF\pi}\Gamma\left(\frac{\DOF}{2}\right)} \left(1+\frac{x^2}{\DOF}\right)^{-\frac{\DOF+1}{2}} \]
Start with the Student-t density. Then standardise it and express it with time-varying volatility:
\[ \frac{1}{\Vol_t} \frac{\Gamma\left(\frac{\DOF+1}{2}\right)}{\sqrt{(\DOF-2)\pi}\Gamma\left(\frac{\DOF}{2}\right)} \left(1+\frac{\CompoundReturns_t^2}{\Vol_t^2(\DOF-2)}\right)^{-\frac{\DOF+1}{2}} \]
Taking logs, the constant part is: \[ \log\Gamma\left(\frac{\DOF+1}{2}\right) - \frac{1}{2}\log((\DOF-2)\pi) - \log\Gamma\left(\frac{\DOF}{2}\right) \]
The time-varying part is: \[ -\log\Vol_t - \frac{\DOF+1}{2}\left(\log\left(\Vol_t^2(\DOF-2)+\CompoundReturns_t^2\right) - \log\left(\Vol_t^2(\DOF-2)\right)\right) \]
Summing over \(t = 2, \ldots, \SampleSize\) gives the full log-likelihood: \[ \log\lik = \underbrace{(\SampleSize-1)\left(\log\Gamma\left(\frac{\DOF+1}{2}\right) - \frac{1}{2}\log((\DOF-2)\pi) - \log\Gamma\left(\frac{\DOF}{2}\right)\right)}_{\text{Constant}} + \sum_{t=2}^{\SampleSize}\left(-\log\Vol_t - \frac{\DOF+1}{2}\left(\log\left(\Vol_t^2(\DOF-2)+\CompoundReturns_t^2\right) - \log\left(\Vol_t^2(\DOF-2)\right)\right)\right) \]
21.4 Implementation
This is where the formulas become an estimator. We combine the likelihood with the optimiser and deal with the parts that usually cause trouble.
21.4.1 Practical considerations
One practical issue arises before estimation starts. The optimiser calls the likelihood function many times, so we pre-calculate squared returns, sample size and variance outside the loop.
Parameter positivity is handled as described in Section 21.2.3.
21.4.2 Normal GARCH estimation
y2 = y_vec^2
sigma2_init = var(y_vec)
T_len = length(y2)
LL_GARCH = function(par) {
par = abs(par)
omega = par[1]
alpha = par[2]
beta = par[3]
sigma2 = sigma2_init
loglikelihood = 0
for (i in 2:T_len) {
sigma2 = omega + alpha * y2[i-1] + beta * sigma2
if (sigma2 <= 0) return(1e10)
loglikelihood = loglikelihood + log(sigma2) + y2[i] / sigma2
}
loglikelihood = -0.5 * loglikelihood - (T_len - 1) / 2 * log(2 * pi)
return(-loglikelihood)
}
opts = list("algorithm" = "NLOPT_LN_NELDERMEAD", "xtol_rel" = 1.0e-8, maxeval = 1e5)
res_garch = nloptr(x0 = c(0.000001, 0.1, 0.85),
eval_f = LL_GARCH,
lb = c(0, 0, 0),
ub = c(1, 1, 1),
opts = opts)
if (res_garch$status < 0) warning("Normal GARCH nloptr did not converge: status ", res_garch$status)
cat("Normal GARCH(1,1) Parameters:\n")
cat("Omega:", sprintf("%.6f", res_garch$solution[1]), "\n")
cat("Alpha:", sprintf("%.6f", res_garch$solution[2]), "\n")
cat("Beta: ", sprintf("%.6f", res_garch$solution[3]), "\n")
cat("Log-likelihood:", sprintf("%.2f", -res_garch$objective), "\n")
r_loglik_garch = -res_garch$objectiveNormal GARCH(1,1) Parameters:
Omega: 0.000002
Alpha: 0.105267
Beta: 0.879347
Log-likelihood: 29234.41y2 = y**2
sigma2_init = np.var(y, ddof=1)
T_len = len(y2)
def LL_GARCH(par):
par = np.abs(par)
omega, alpha, beta = par
sigma2 = sigma2_init
loglikelihood = 0
for i in range(1, T_len):
sigma2 = omega + alpha * y2[i-1] + beta * sigma2
if sigma2 <= 0:
return 1e10
loglikelihood += np.log(sigma2) + y2[i] / sigma2
loglikelihood = -0.5 * loglikelihood - (T_len - 1) / 2 * np.log(2 * np.pi)
return -loglikelihood
res_garch = minimize(LL_GARCH, x0=[0.000001, 0.1, 0.85], method='Nelder-Mead',
options={'xatol': 1e-8, 'fatol': 1e-8, 'maxiter': 10000})
if not res_garch.success:
print(f"Warning: Normal GARCH did not converge: {res_garch.message}")
print("Normal GARCH(1,1) Parameters:")
print(f"Omega: {np.abs(res_garch.x[0]):.6f}")
print(f"Alpha: {np.abs(res_garch.x[1]):.6f}")
print(f"Beta: {np.abs(res_garch.x[2]):.6f}")
print(f"Log-likelihood: {-res_garch.fun:.2f}")
py_loglik_garch = -res_garch.funNormal GARCH(1,1) Parameters:
Omega: 0.000002
Alpha: 0.105267
Beta: 0.879347
Log-likelihood: 29234.41using Statistics, Optim, Printf
y2_jl = y_jl.^2;
h_jl = var(y_jl);
N_jl = length(y2_jl);
function gLikelihoodA(theta, y2, h)
o = abs(theta[1])
a = abs(theta[2])
b = abs(theta[3])
N = length(y2)
lik = 0.0
u = h
for i in 2:N
u = o + a * y2[i-1] + b * u
if u <= 0 || !isfinite(u)
return 1e10
end
lik += log(u) + y2[i] / u
end
nll = -(-0.5 * lik - (N - 1) * 0.5 * log(2 * pi))
return isfinite(nll) ? nll : 1e10
end
lg(b) = gLikelihoodA(b, y2_jl, h_jl);
res_garch_jl = optimize(lg, [0.000001, 0.1, 0.85], NelderMead(),
Optim.Options(iterations = 100000));
if !Optim.converged(res_garch_jl)
println("Warning: Normal GARCH did not converge")
end
params_garch = abs.(Optim.minimizer(res_garch_jl));
jl_loglik_garch = -Optim.minimum(res_garch_jl);
println("Normal GARCH(1,1) Parameters:")
@printf("Omega: %.6f\n", params_garch[1])
@printf("Alpha: %.6f\n", params_garch[2])
@printf("Beta: %.6f\n", params_garch[3])
@printf("Log-likelihood: %.2f\n", jl_loglik_garch)Normal GARCH(1,1) Parameters:
Omega: 0.000002
Alpha: 0.105267
Beta: 0.879347
Log-likelihood: 29234.4121.4.3 GARCH-t estimation
LL_tGARCH = function(par) {
par = abs(par)
omega = par[1]
alpha = par[2]
beta = par[3]
nu = par[4]
n1 = nu + 1.0
n2 = nu - 2.0
sigma2 = sigma2_init
loglikelihood = 0
for (i in 2:T_len) {
sigma2 = omega + alpha * y2[i-1] + beta * sigma2
if (sigma2 <= 0) return(1e10)
loglikelihood = loglikelihood + log(sigma2) +
n1 * (log(n2 * sigma2 + y2[i]) - log(n2 * sigma2))
}
loglikelihood = -loglikelihood / 2.0
loglikelihood = loglikelihood + (T_len - 1) * (lgamma(n1/2) -
lgamma(nu/2) - 0.5 * log(pi * n2))
return(-loglikelihood)
}
opts = list("algorithm" = "NLOPT_LN_NELDERMEAD", "xtol_rel" = 1.0e-8, maxeval = 1e5)
res_tgarch = nloptr(x0 = c(0.000001, 0.1, 0.85, 6),
eval_f = LL_tGARCH,
lb = c(0, 0, 0, 2.1),
ub = c(1, 1, 1, 100),
opts = opts)
if (res_tgarch$status < 0) warning("tGARCH nloptr did not converge: status ", res_tgarch$status)
cat("GARCH-t(1,1) Parameters:\n")
cat("Omega:", sprintf("%.6f", res_tgarch$solution[1]), "\n")
cat("Alpha:", sprintf("%.6f", res_tgarch$solution[2]), "\n")
cat("Beta: ", sprintf("%.6f", res_tgarch$solution[3]), "\n")
cat("Nu: ", sprintf("%.6f", res_tgarch$solution[4]), "\n")
cat("Log-likelihood:", sprintf("%.2f", -res_tgarch$objective), "\n")
r_loglik_tgarch = -res_tgarch$objectiveGARCH-t(1,1) Parameters:
Omega: 0.000001
Alpha: 0.099790
Beta: 0.894967
Nu: 6.304794
Log-likelihood: 29443.61def LL_tGARCH(par):
par = np.abs(par)
omega, alpha, beta, nu = par
if nu < 2.1 or nu > 100:
return 1e10
n1 = nu + 1.0
n2 = nu - 2.0
sigma2 = sigma2_init
loglikelihood = 0
for i in range(1, T_len):
sigma2 = omega + alpha * y2[i-1] + beta * sigma2
if sigma2 <= 0 or not np.isfinite(sigma2):
return 1e10
loglikelihood += np.log(sigma2) + n1 * (np.log(n2 * sigma2 + y2[i]) - np.log(n2 * sigma2))
loglikelihood = -loglikelihood / 2.0
loglikelihood += (T_len - 1) * (gammaln(n1/2) - gammaln(nu/2) - 0.5 * np.log(np.pi * n2))
nll = -loglikelihood
return nll if np.isfinite(nll) else 1e10
res_tgarch = minimize(LL_tGARCH, x0=[0.000001, 0.1, 0.85, 6], method='Nelder-Mead',
options={'xatol': 1e-8, 'fatol': 1e-8, 'maxiter': 10000})
if not res_tgarch.success:
print(f"Warning: tGARCH did not converge: {res_tgarch.message}")
print("GARCH-t(1,1) Parameters:")
print(f"Omega: {np.abs(res_tgarch.x[0]):.6f}")
print(f"Alpha: {np.abs(res_tgarch.x[1]):.6f}")
print(f"Beta: {np.abs(res_tgarch.x[2]):.6f}")
print(f"Nu: {np.abs(res_tgarch.x[3]):.6f}")
print(f"Log-likelihood: {-res_tgarch.fun:.2f}")
py_loglik_tgarch = -res_tgarch.funGARCH-t(1,1) Parameters:
Omega: 0.000001
Alpha: 0.099790
Beta: 0.894967
Nu: 6.304796
Log-likelihood: 29443.61using Optim, SpecialFunctions, Printf
function LL_tGARCH_jl(theta, y2, h, N)
o = abs(theta[1])
a = abs(theta[2])
b = abs(theta[3])
nu = abs(theta[4])
# Constraint: nu in [2.1, 100], matching R's nloptr bounds
if nu < 2.1 || nu > 100
return 1e10
end
n1 = nu + 1.0
n2 = nu - 2.0
u = h
loglik = 0.0
for i in 2:N
u = o + a * y2[i-1] + b * u
if u <= 0 || !isfinite(u)
return 1e10
end
loglik += log(u) + n1 * (log(n2 * u + y2[i]) - log(n2 * u))
end
loglik = -loglik / 2.0
loglik += (N - 1) * (loggamma((nu + 1) / 2) - loggamma(nu / 2) - 0.5 * log(pi * (nu - 2)))
nll = -loglik # Negative log-likelihood for minimisation
return isfinite(nll) ? nll : 1e10
end
# Use NelderMead
ll(b) = LL_tGARCH_jl(b, y2_jl, h_jl, N_jl);
res_tgarch_jl = optimize(ll, [0.000001, 0.1, 0.85, 6.0], NelderMead(),
Optim.Options(iterations = 100000));
if !Optim.converged(res_tgarch_jl)
println("Warning: tGARCH did not converge")
end
params_tgarch = abs.(Optim.minimizer(res_tgarch_jl));
jl_loglik_tgarch = -Optim.minimum(res_tgarch_jl);
println("GARCH-t(1,1) Parameters:")
@printf("Omega: %.6f\n", params_tgarch[1])
@printf("Alpha: %.6f\n", params_tgarch[2])
@printf("Beta: %.6f\n", params_tgarch[3])
@printf("Nu: %.6f\n", params_tgarch[4])
@printf("Log-likelihood: %.2f\n", jl_loglik_tgarch)GARCH-t(1,1) Parameters:
Omega: 0.000001
Alpha: 0.099788
Beta: 0.894968
Nu: 6.304814
Log-likelihood: 29443.6121.5 Comparison
The table below compares the log-likelihood values across the three implementations. R uses box-bounded optimisation, so every parameter stays inside explicit lower and upper limits. Python and Julia instead use unbounded Nelder-Mead with absolute-value reflection and penalty terms in the likelihood. This keeps parameters non-negative and the degrees of freedom within range.
None of the three enforces covariance stationarity, \(\ARCHcoeff + \GARCHcoeff < 1\). R’s bounds constrain each parameter separately, not their sum, so the two coefficients can still sum above one. R, Python and Julia each report successful convergence for this sample, and the log-likelihoods agree to within rounding in the table below.
|Model | R| Python| Julia|
|:-----------------|--------:|--------:|--------:|
|Normal GARCH(1,1) | 29234.41| 29234.41| 29234.41|
|GARCH-t(1,1) | 29443.61| 29443.61| 29443.61|The Student-t model achieves a higher log-likelihood than the normal model. A more flexible model always fits at least as well as one it nests, and the normal GARCH is the \(\DOF \to \infty\) limit of the GARCH-t, so part of the gain is mechanical. The improvement here is large enough to point to fat tails. Financial returns produce extreme outcomes to which the normal distribution assigns too little probability. A likelihood ratio test or an information criterion, not the log-likelihood alone, is what would settle the question formally.
21.6 Exercise
Implement a manual GJR-GARCH(1,1) estimation. Modify the likelihood function to include the leverage parameter \(\GJRleverage\), and compare your results across the three languages. Discuss the additional optimisation challenges introduced by the asymmetric term.
The GJR-GARCH model is: \[ \Vol_t^2 = \GARCHconst + \ARCHcoeff \CompoundReturns_{t-1}^2 + \GJRleverage \CompoundReturns_{t-1}^2 \Indicator_{\CompoundReturns_{t-1}<0} + \GARCHcoeff \Vol_{t-1}^2 \]
where \(\Indicator_{\CompoundReturns_{t-1}<0}\) is an indicator function that equals 1 when returns are negative.