Efficient Method of Moments(EMM) for Stochastic volatility model

Question

We are attempting to calibrate the parameters of the Heston model via EMM on historical stock price returns. However, we are first trying a simple stochastic volatility model using EMM.

We have come across an example of the EMM procedure in SAS, which corresponds to one of the models in the paper by Anderson.

Link to Andersen paper: https://www.uh.edu/~bsorense/Emmpublished.pdf

Link to SAS example: https://support.sas.com/rnd/app/ets/examples/emmweb/index.htm#:~:text=The%20Efficient%20Method%20of%20Moments,Method%20of%20Moments%20(GMM.)

We are attempting to replicate this model in Python as we are unfamiliar with SAS. We are able to fit the GARCH(1,1) parameters correctly. We compute the covariance matrix from the first moments (scores) of the GARCH model correctly. However, when attempting to do the Simulated Method of Moments (SMM) component of the SAS implementation, we are unable to align on the calibrated parameters for a, b and s. We suspect that there is an issue in the minimization algorithm in our Python implementation. Any insights into where we are going wrong would be greatly appreciated.

** NOTE: The TEST.csv file contains the returns[y] simulated in SAS(seen as historical returns). The EMM_u_z.csv file contains the simulated normal random variables[u and z] for the simulated method of moments from SAS(unique seed values)**

import pandas as pd

from sci..py.optimize import minimize

from arch import arch_model

directory = r'C:\Users\Heston'
filename = '\\TEST.csv'
filename2 = '\\EMM_u_z.csv'

data = pd.read_csv(directory+filename)
u_z = pd.read_csv(directory + filename2)

nobs = 1000
nsim = 20000

y = data['y'].values

print(np.mean(y))
print(np.std(y))

# Estimate GARCH(1,1) model
model = arch_model(y, vol='Garch', lags = 1, p=1, q=1, )

result = model.fit()
print(result.summary)

gsigmasq = result.conditional_volatility**2

score_1 = np.zeros(nobs)
score_2 = np.zeros(nobs)
score_3 = np.zeros(nobs)

for i in range(nobs):
   if i == 0:
       score_1[i] = (-1 + y[i]**2/gsigmasq[i])/gsigmasq[i]
       score_2[i] = 0
       score_3[i] = 0
   else:
       score_1[i] = (-1 + y[i]**2/gsigmasq[i])/gsigmasq[i]
       score_2[i] = (-1 + y[i]**2/gsigmasq[i])*gsigmasq[i-1]/gsigmasq[i]
       score_3[i] = (-1+y[i]**2/gsigmasq[i])*y[i-1]**2/gsigmasq[i]
       
score = np.vstack([score_1,score_2,score_3]).T

MSE = np.mean(result.resid**2)

#Volatility of GARCH model
params_fitted = result.params
omega = params_fitted['omega']
alpha = params_fitted['alpha[1]']
beta = params_fitted['beta[1]']
mu = params_fitted['mu']

# Parameter initialization
a_init, b_init, s_init = -0.736, 0.9, 0.363
params_init = np.array([a_init, b_init, s_init])
u_z['_ah_0'] = omega
u_z['_gh_1'] = beta
u_z['_ah_1'] = alpha
u_z['_mse_'] = MSE

# Bounds for optimization
bounds = [(-np.inf, np.inf), (-1, 1), (0, np.inf)]  # (a, b, s)

# GMM objective function
def gmm_objective(params, data, scores):
   a, b, s = params
   lsigmasq = np.exp(a)  # Replace with appropriate lag calculation
   lnsigmasq = a + b * np.log(lsigmasq) + s * data['u']
   sigmasq = np.exp(lnsigmasq)
   ysim = np.sqrt(sigmasq) * data['z']
   
   gsigmasq = data['_ah_0'] + data['_gh_1'] * data['_mse_'].shift(1) + data['_ah_1'] * (ysim ** 2)

   m1 = (-1 + (ysim ** 2) / gsigmasq) / gsigmasq
   m2 = (-1 + (ysim ** 2) / gsigmasq) * gsigmasq.shift(1) / gsigmasq
   m3 = (-1 + (ysim ** 2) / gsigmasq) * (ysim ** 2).shift(1) / gsigmasq

   moments = np.column_stack((m1, m2, m3))
   
   # GMM moment conditions
   gmm_conditions = np.nanmean(moments, axis = 0) # Use sample moments
   
   # Compute the variance-covariance matrix
   vhat = np.cov(scores, rowvar=False)

   gmm_objective_value = gmm_conditions @ np.linalg.inv(vhat) @ gmm_conditions.T
   
   return gmm_objective_value  # GMM objective function

# GMM estimation
result = minimize(gmm_objective, params_init, args=(u_z,score), method='BFGS', bounds = bounds, tol = 10e-12)

# Extract estimated parameters
params_estimated = result.x

# Display results
print("Estimated Parameters:")
print("a:", params_estimated[0])
print("b:", params_estimated[1])
print("s:", params_estimated[2])
print("\nMinimization Result:")
print(result)