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)