随机波动率模型的有效矩量法（EMM）

问：

我们正试图通过EMM根据历史股价回报校准Heston模型的参数。然而，我们首先尝试使用EMM的简单随机波动率模型。

我们在SAS中遇到了一个EMM过程的例子，它与Anderson论文中的一个模型相对应。

链接到安徒生论文：https://www.uh.edu/~bsorense/Emmpublished.pdf

链接到 SAS 示例：https://support.sas.com/rnd/app/ets/examples/emmweb/index.htm#:~:text=The%20Efficient%20Method%20of%20Moments,Method%20of%20Moments%20（GMM.）

我们试图在 Python 中复制这个模型，因为我们不熟悉 SAS。我们能够正确地拟合 GARCH（1,1） 参数。我们从GARCH模型的第一个时刻（分数）正确计算协方差矩阵。但是，在尝试执行 SAS 实现的矩模拟法 （SMM） 组件时，我们无法对齐 a、b 和 s 的校准参数。我们怀疑 Python 实现中的最小化算法存在问题。任何关于我们哪里出错的见解将不胜感激。

**注意： TEST.csv 文件包含在 SAS 中模拟的返回[y]（被视为历史返回）。EMM_u_z.csv 文件包含模拟正态随机变量 [u 和 z]，用于模拟 SAS（唯一种子值）中的矩**

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)