Problem with the maximum likelihood for a GARCH-type of model

I'm currently working with the following GARCH process from Heston and Nandi (2000): \begin{align*} r_{t+1} - r_f &= \lambda h_{t+1} - \frac{h_{t+1}}{2} + \sqrt{h_{t+1}}z_{t+1} \\ h_{t+1} &= \omega + \beta h_t + \alpha \left( z_t - \gamma \sqrt{h_t} \right)^2 \end{align*} given $$z_{t+1} \sim N(0,1)$$, we can estimate the model parameters by maximum likelihood. I wrote some python code to simulate the process and, then, to compute the likelihood at the hypothesized parameter values. The density for each observation is given by:

$$$$f(r_{t+1} - r| h_{t+1}) = \frac{1}{\sqrt{2 \pi h_{t+1}}} \exp \left( \frac{-(r_{t+1} - r - \lambda h_{t+1} + \frac{h_{t+1}}{2})^2}{2 h_{t+1}} \right)$$$$

My problem is how should I work out the maximization? Obviously, $$|1 - \beta - \alpha \gamma^2 | < 1$$ ensures that the conditional variance process is covariance stationary. Moreover, $$(\omega + \alpha)/(1 - \beta - \alpha \gamma^2) > 0$$ ensures that unconditional variance is positive. So, as some might suspect, I am having trouble making sure the maximization algorithm can converge toward the real parameter values and I am looking for how I should approach this issue.

import numpy as np
from   numpy import sqrt, exp, log
from   matplotlib.pyplot import plot, hist
from   statistics import mean
from scipy.optimize import minimize
#%%

r       = 0.05/252

param   = [-9.765e-07, 2.194e-06, 0.8986, 205.15, 3.930]
omega, alpha, beta, gamma, Lambda = param

sigma2 = (omega+alpha)/(1-beta-alpha*gamma**2)
h0     = sigma2

T = 1000
z = np.random.normal(loc=0, scale=1, size=T)
R = np.zeros(shape=T)
h = h0*np.ones(shape=T)

for tt in range(0,T-1):
h[tt+1] = omega + beta*h[tt] + alpha*(z[tt] - gamma*sqrt(h[tt]))**2
R[tt+1] = r + Lambda*h[tt+1] - h[tt+1]/2 + sqrt(h[tt+1])*z[tt+1]

hh = h
Rt = R - r

def TS_Loglik_HN(Rt, h0, param):
'''
Author: Stéphane Surprenant, UQAM
Creation: 02/04/2020

Description: This function returns the value of the log-likelihood for the
Heston and Nandi (2000) process under the physical measure.

INPUTS       DESCRIPTION
Rt :         (float) Series of (log) returns minus the risk-free rate.
h0 :         (float) Initial value of the variance (Daily)
param:       (float) Parameters of the model
[omega, alpha, beta, gamma, Lambda] = param

OUTOUTS      DESRIPTION
loglik       (float) Log-likelihood value

Model:
Rt[tt+1] := R[tt+1] - r
= Lambda*h[tt+1] - h[tt+1]/2 + sqrt(h[tt+1])*z[tt+1]
h[tt+1]   = omega + beta*h[tt] + alpha*(z[tt] - gamma*sqrt(h[tt]))**2
'''
# Assign parameter values
omega, alpha, beta, gamma, Lambda = param

# Initialize matrices
T    = len(Rt)
h    = h0*np.ones(shape=T)
e    = np.zeros(shape=T)

# Filtering volatility
for tt in range(0,T-1):
e[tt]   = (Rt[tt] - Lambda*h[tt] + h[tt]/2)/sqrt(h[tt])
h[tt+1] = omega + beta*h[tt] + alpha*(e[tt] - gamma*sqrt(h[tt]))**2

e[T-1]   = (Rt[T-1] - Lambda*h[T-1] + h[T-1]/2)/sqrt(h[T-1])

# Compute Log-likelihood
l = -0.5*(log(2*np.pi) + log(h) + e**2)
loglik = sum(l)

return(loglik)

# Example:
f = lambda x: -TS_Loglik_HN(Rt, h0, x)
results = minimize(f, param)

• I have never done any kind of optimization with Python so far, but shouldn't you simply introduce your constraints to the optimiser (docs.scipy.org/doc/scipy/reference/generated/…), see section on constraints. ? – Kermittfrog Apr 3 at 9:19
• I am looking for what is typically done in those circumstances. – Stéphane Apr 3 at 14:25
• I just did that in R with a ‚standard‘ GARCH(1,1) today and all I had to do was to supply the usual linear constraints to the optimizer. – Kermittfrog Apr 3 at 17:46
• A VERY brute force approach would be to simply let your likelihood become VERY negative if a constraint is violated. Again, that is not the way to go, but sometimes quick&dirty suffices... – Kermittfrog Apr 3 at 17:47
• I'm not familiar with this model, do you know the necessary or sufficient conditions to obtain a positive cond. variance ? If yes, it is common to use an algorithm that allows imposing these conditions (non-linear constraints) during estimation (ex: see method='SLSQP' in scipy.optimize). – Malick Apr 3 at 21:04

If you have the path likelihood, you can try just writing that function and optimizing it directly. You might have some issues with the variance piece. This looks a lot like parameter inference for SDE, data-assimilation etc.

I think if you write a proper likelihood function with priors for all parameters and same via some MCMC or MC (Gibbs) that is guaranteed to work for you.

You can also try a variational inference approach and just optimize for MLE of params.

If you write out the Likelihood above (in latex) it might be easier to discuss and notice any stability issues.

UPDATE:

So for pure MLE approach you can just try to optimize the log-likelihood as you are doing. If it's not converging maybe try to do stability analysis. A quick sanity test is if you start near the real values (you know them in this case since you generated them) and see if it converges. Calculating the hessian might give some insight too but this is basically stability analysis. Another debug is to just try to fit one parameter at a time with all the other parameters given correctly or at least nearby the correct values. I would be a bit worried about $$h$$ being close to zero but I haven't fully grasped the process so maybe that is ok.

Started messing around with the code and either I introduced a bug and then fixed it or you have an off-by-one error. Either way, you might want to add the same checks. Basically I'm just checking that I can back out $$h$$ and $$z$$ (your e[tt]) properly.

from statistics import mean

import numpy as np
from numpy import exp, log, sqrt
from pylab import *
from scipy.optimize import minimize

r = 0.05 / 252
param = [9.765e-07, 2.194e-06, 0.8986, 205.15, 3.930]
omega, alpha, beta, gamma, Lambda = param

def get_h0(param):
omega, alpha, beta, gamma, Lambda = param
sigma2 = (omega + alpha) / (1 - beta - alpha * gamma ** 2)
h0 = sigma2
return h0

h0 = get_h0(param)

def rhs_h(param, h, z):
omega, alpha, beta, gamma, Lambda = param
return omega + beta * h + alpha * (z - gamma * sqrt(h)) ** 2

def rhs_R(param, h, z):
omega, alpha, beta, gamma, Lambda = param
return Lambda * h - h / 2 + sqrt(h) * z

def get_paths(param):
omega, alpha, beta, gamma, Lambda = param
assert omega > 0
assert alpha > 0
assert beta > 0
assert beta + alpha * gamma ** 2 < e
np.random.seed(10)
T = 10
z = np.random.normal(loc=0, scale=1, size=T)
R = np.zeros(shape=T - 1)
h = h0 * np.ones(shape=T)
for i in range(0, T - 1):
h[i + 1] = rhs_h(param, h[i], z[i])
R[i] = r + rhs_R(param, h[i], z[i])
return R, h, z

def get_h_z_from_R(Rt, h0, param):
omega, alpha, beta, gamma, Lambda = param
T = len(Rt)
h = np.empty(shape=T)
h[0] = h0
z = np.zeros(shape=T)
for i in range(0, T - 1):
z[i] = (Rt[i] - Lambda * h[i] + h[i] / 2) / sqrt(h[i])
h[i + 1] = omega + beta * h[i] + alpha * (z[i] - gamma * sqrt(h[i])) ** 2
z[T - 1] = (Rt[T - 1] - Lambda * h[T - 1] + h[T - 1] / 2) / sqrt(h[T - 1])
return h, z

R, h, z = get_paths(param)
Rt = R - r

h_check, z_check = get_h_z_from_R(Rt, h0, param)
assert np.allclose(z[:-1], z_check)
assert np.allclose(h[:-1], h_check)

• This function here is the log likelihood. My problem is what would be the best way to make sure the optimization algorithm doesn't input nonsensical values that does things like assigning a negative value to conditional variance at some point in the sample. – Stéphane Apr 3 at 14:24
• The proper prior on the parameters would directly encapsulate the constraints. You’d swap speed with memory... but an MCMC like estimation is of course always nice... – Kermittfrog Apr 3 at 17:48
• I appreciate the comment, @Kermittfrog, but I have to work this out from a frequentist perspective. It has to be done using maximum likelihood. – Stéphane Apr 3 at 23:08
• I think @Stéphane might need help around the constraints piece, if I am not mistaken. But your stepwise approach is still a valuable hint. – Kermittfrog Apr 5 at 12:38
• Thanks for the update. $h_t$ is the daily variance, so the yearly volatility is $\sqrt{h_t*252}$. For 18% yearly volatility, you get a daily value of 0,000128571, so, yes, it's always pretty small. – Stéphane Apr 5 at 21:35

If I am not mistaken, as you have already stated you have the long run relationship

$$h\left(1-\beta-\alpha\gamma^2\right)=\omega + \alpha$$

I suggest you impose the following restrictions that should ensure $$h_t$$ to stay positive:

\begin{align} \omega&>0\\ \alpha&>0\\ \beta &>0\\ \beta+\alpha\gamma^2&<1\\ \end{align}

I suspect you do not have to impose any restrictions on $$\gamma$$ per se. Unfortunately, these restrictions cannot be written in terms of linear (in)equality restrictions, but that should not be that much of a problem, really. HTH

PS: In practice, of course, you set your boundaries to something like 1E-6 or 1-1E-6.