Efficient Method of Moments(EMM) for Stochastic volatility model

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.

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 insight into where we are going wrong would be greatly appreciated.

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 scipy.optimize import minimize

from arch import arch_model

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

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)
$$$$
`
• I’ve never used EMM, but in my experience with calibrating SV models is that they’re very sensitive to your initial guesses and scaling/damping factors. I haven’t checked your code, but if you have implemented it correctly, you are most likely experiencing the drawbacks of the Heston - it’s incredibly unstable and lacks tractability. Heston has a book on Heston calibration with matlab and C code. Nov 14, 2023 at 11:22