2
$\begingroup$

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)
Improve this question
$\endgroup$
7
  • $\begingroup$ 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. ? $\endgroup$
    – Kermittfrog
    Apr 3, 2020 at 9:19
  • $\begingroup$ I am looking for what is typically done in those circumstances. $\endgroup$
    – Stéphane
    Apr 3, 2020 at 14:25
  • $\begingroup$ 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. $\endgroup$
    – Kermittfrog
    Apr 3, 2020 at 17:46
  • $\begingroup$ 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... $\endgroup$
    – Kermittfrog
    Apr 3, 2020 at 17:47
  • $\begingroup$ 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). $\endgroup$
    – Malick
    Apr 3, 2020 at 21:04

2 Answers 2

Reset to default
0
$\begingroup$

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)
Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Stéphane
    Apr 3, 2020 at 14:24
  • $\begingroup$ 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... $\endgroup$
    – Kermittfrog
    Apr 3, 2020 at 17:48
  • $\begingroup$ I appreciate the comment, @Kermittfrog, but I have to work this out from a frequentist perspective. It has to be done using maximum likelihood. $\endgroup$
    – Stéphane
    Apr 3, 2020 at 23:08
  • 1
    $\begingroup$ 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. $\endgroup$
    – Kermittfrog
    Apr 5, 2020 at 12:38
  • $\begingroup$ 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. $\endgroup$
    – Stéphane
    Apr 5, 2020 at 21:35
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.