Heston Nandi Garch Implementation Problem for Python

Question

I have a coded my own Garch class in order to implement the Heston-Nandi Garch model.

import numpy as np
import pandas as pd
import math
from scipy.optimize import *

Ticker=['^GSPC']
df=pd.DataFrame()
for i in Ticker:
      df[i]=wb.DataReader(i,'yahoo', start='01-09-2012',end = '01-09-2014')['Adj Close']

df2 = np.log(df/df.shift(1)).dropna() * 100

This code extracts the data and rescale them by multiplying the rate of returns by 100 for stability purpose(otherwise we end up with a messy log likelihood function to optimize).The following class compute our Garch.

class NGARCH :
    def __init__(self,df,r = 0):
        
        self.df = np.array(df)
        self.r = r
        self.NGarch_Log_Like()
        
    def NGarch(self,params):
        w,Beta,alpha,y,lambda_  = params
        h = np.zeros(len(self.df))
        e = np.zeros(len(self.df))
        for i in range(0,len(self.df)):
            if i == 0:
                h[0] = (w + alpha) / ( 1 - alpha * y**2 - Beta ) #Unconditional Variance 
                if h[0] < 0 or 1 - alpha * y**2 - Beta < 1e-3 :#or w > 1e-8 or alpha > 1e-05 or w < 0 or alpha < 0 or Beta > 0.90 :
                    return 1e50
            else :
                h[i] = w + Beta * h[i-1] + alpha * (e[i-1] - y * np.sqrt(h[i-1]))**2
                e[i] = (self.df[i] - self.r - lambda_ * h[i] + h[i] / 2)/np.sqrt(h[i])
        LogL = - .5 * np.sum(-np.log( 2 * np.pi) - np.log(h[:]) - (e[:])**2)
        return LogL
    
    def NGarch_Model(self):
        h = np.zeros(len(self.df))
        e = np.zeros(len(self.df))
        for i in range(0,len(self.df)):
            if i == 0:
                h[0] = (self.w + self.alpha) / (1 - self.alpha * self.y**2 - self.Beta)
            else :
                h[i] =  self.w + self.Beta * h[i-1] + self.alpha * (e[i-1]- self.y * np.sqrt(h[i-1]))**2
                e[i] = (self.df[i] - self.r - self.lambda_ * h[i] + h[i] / 2)/np.sqrt(h[i])
        return h,e
   
    def NGarch_Log_Like(self):
        #w,Beta,alpha,y,lambda_ = 9.765e-10, 0.90,2.194e-06,100.15,10
        #x0 = np.array([w,Beta,alpha,y,lambda_])
        #cons = ({'type': 'ineq', 'fun': lambda x: 1-x[1]-x[2]*(x[3]**2)})#{'type': 'ineq', 'fun': lambda x: np.array(x)})
        bounds = ((1e-8,1e-120),(1e-6,1),(1e-6,1e-12),(1e-10,1e3),(1e-10,100))
        #sol = minimize(self.NGarch,x0,bounds = bounds,constraints = cons,tol = 1e-9,method = 'Nelder-Mead')
        sol = differential_evolution(self.NGarch,bounds = bounds,tol = 1e-9,polish = False,workers = -1)
        self.w,self.Beta,self.alpha,self.y,self.lambda_ = sol.x
        print("Log Likelihood : {} ".format(-sol.fun))
        print(100 * "-")
        print("Stationarity Condition distance : ", 1-(self.Beta + self.alpha * self.y**2 ))
        print(100 * "-")
        print(" Parameters : ",sol.x)
        print(100 * "-")

However if you insert:

f = NGARCH(df2,0)

The result is such :

Log Likelihood : -9949.302700236756 
----------------------------------------------------------------------------------------------------
Stationarity Condition distance :  0.0010000008123174364
----------------------------------------------------------------------------------------------------
 Parameters :  [9.99632246e-09 9.17942150e-01 9.99999894e-07 2.84706616e+02
 1.13789977e+01]
----------------------------------------------------------------------------------------------------

The log likelihood is terrible for a dataframe with less than 700 data and it seems like the differential evolution algorithm is maximizing the stationary condition’s constraint and therefore our generated variance is always barely stationary.Better log likelihood appears if we are using the Local optimisor(Nelder-Mead).However the process stays barely stationary.

Could someone help because I was stuck with that during the whole day.

Thank you.

PS : The whole math behind that is available on the following link: Problem with the maximum likelihood for a GARCH-type of model