Problema de máxima verosimilitud para un modelo de tipo GARCH

Question

Actualmente estoy trabajando con el siguiente proceso GARCH de Heston y 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*}$$ dado $z_{t+1} \sim N(0,1)$ podemos estimar los parámetros del modelo por máxima verosimilitud. Escribí un código python para simular el proceso y, a continuación, calcular la verosimilitud en los valores hipotéticos de los parámetros. La densidad para cada observación viene dada por:

$$\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}$$

Mi problema es ¿cómo debo calcular la maximización? Obviamente, $|1 - \beta - \alpha \gamma^2 | < 1$ asegura que el proceso de varianza condicional es estacionario en la covarianza. Además, $(\omega + \alpha)/(1 - \beta - \alpha \gamma^2) > 0$ garantiza que la varianza incondicional sea positiva. Así que, como algunos podrían sospechar, estoy teniendo problemas para asegurarme de que el algoritmo de maximización pueda converger hacia los valores reales de los parámetros y estoy buscando cómo debería abordar esta cuestión.

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 1

Si no me equivoco, como ya has dicho tienes la relación a largo plazo

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

Le sugiero que imponga las siguientes restricciones que deberían garantizar $h_t$ para mantenerse positivo:

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

Sospecho que no tiene que imponer ninguna restricción a $\gamma$ per se. Por desgracia, estas restricciones no pueden escribirse en términos de restricciones de (in)igualdad lineal, pero eso no debería ser un gran problema, en realidad. HTH

P.D.: En la práctica, por supuesto, se fijan los límites en algo así como 1E-6 o 1-1E-6.

Answer 2

Si tienes la probabilidad de la ruta, puedes intentar escribir esa función y optimizarla directamente. Podrías tener algunos problemas con la parte de la varianza. Esto se parece mucho a la inferencia de parámetros para SDE, asimilación de datos, etc.

Creo que si escribes una función de verosimilitud adecuada con priores para todos los parámetros y la misma a través de algún MCMC o MC (Gibbs) que está garantizado para trabajar para usted.

También se puede intentar un enfoque de inferencia variacional y simplemente optimizar para el MLE de los parámetros.

Si escribes la Probabilidad arriba (en látex) podría ser más fácil discutir y notar cualquier problema de estabilidad.

ACTUALIZACIÓN:

Por lo tanto, para el enfoque MLE puro puede tratar de optimizar la log-verosimilitud como lo está haciendo. Si no converge, tal vez intente hacer un análisis de estabilidad. Una prueba rápida de cordura es si empiezas cerca de los valores reales (los conoces en este caso ya que los has generado) y ves si converge. El cálculo de la hessiana podría dar alguna idea también, pero esto es básicamente el análisis de estabilidad. Otra forma de depuración es tratar de ajustar un parámetro a la vez con todos los demás parámetros dados correctamente o al menos cerca de los valores correctos. Yo estaría un poco preocupado por $h$ siendo casi nulo, pero no he comprendido del todo el proceso, así que tal vez esté bien.

He empezado a trastear con el código y, o bien he introducido un error y luego lo he arreglado, o bien tienes un error de fuera a dentro. De cualquier manera, es posible que desee añadir las mismas comprobaciones. Básicamente estoy comprobando que puedo retroceder $h$ y $z$ (su e[tt]) correctamente.

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)