Estimating Heston–Nandi GARCH Models with Constrained Likelihood

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


# 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:

\begin{equation} 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) \end{equation}

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)
```

## Answer by safetyduck (score 0, accepted)

https://quant.stackexchange.com/a/53019

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)
```

## Answer by Kermittfrog (score 1)

https://quant.stackexchange.com/a/53060

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.