时间序列波动率建模及预测

本文是对https://www.ricequant.com/community/topic/992/garch模型-对波动率建模
的修改和完善,追踪quantopian的原文,有的地方说的不清楚(不完善),所以给大家一个完整的时间序列波动率建模过程(如有错误请指出)。在时间序列中,$a_t$称为新信息或者扰动;在原文中的Garch(1,1)模型的新信息为$x_t$,我们的波动率模型是对新信息的$\sigma_t^2$建模,我们用蒙特卡洛方法模拟出一个GARCH(1,1)过程,变量的动态变化过程如原公式,此外原公式省略了$r_t=\mu+a_t$这个方程,并且假设$\mu$等于零,也是为了标准化收益率做铺垫,从而直接用模拟生成的$a_t$也就是$x_t$来寻求异常值

In [55]:
import cvxopt
from functools import partial
import math
import numpy as np
import scipy
from scipy import stats
import statsmodels as sm
from statsmodels import regression
from statsmodels.stats.stattools import jarque_bera

import matplotlib.pyplot as plt
In [56]:
# 参数初始化
a0 = 1.0
a1 = 0.1
b1 = 0.8
sigma1 = math.sqrt(a0 / (1 - a1 - b1))
In [57]:
def simulate_GARCH(T, a0, a1, b1, sigma1):
    
    # 初始化
    X = np.ndarray(T)
    sigma = np.ndarray(T)
    sigma[0] = sigma1
    
    for t in range(1, T):
        # 产生下一个Xt
        X[t - 1] = sigma[t - 1] * np.random.normal(0, 1)
        # 产生下一个sigma_t
        sigma[t] = math.sqrt(a0 + b1 * sigma[t - 1]**2 + a1 * X[t - 1]**2)
        
    X[T - 1] = sigma[T - 1] * np.random.normal(0, 1)    
    
    return X, sigma
In [58]:
X, _ = simulate_GARCH(10000, a0, a1, b1, sigma1) #生成模拟序列
X = X[1000:] # Drop burn in
X = X / np.std(X) # 标准化X

def compare_tails_to_normal(X):
    # Define matrix to store comparisons
    A = np.zeros((2,4))
    for k in range(4):
        A[0, k] = len(X[X > (k + 1)]) / float(len(X)) # 估计X的尾部,将X异常值得比率求出来
        A[1, k] = 1 - stats.norm.cdf(k + 1) # 正态分布的尾部的累计分布,与上面的值比较,上面的应该大于下面的表示异常值多
    return A

compare_tails_to_normal(X)#(而我们看到大部分确实是这样)
Out[58]:
array([[  1.54666667e-01,   2.53333333e-02,   1.44444444e-03,
          1.11111111e-04],
       [  1.58655254e-01,   2.27501319e-02,   1.34989803e-03,
          3.16712418e-05]])
In [59]:
X2 = np.random.normal(0, 1, 9000)
# 将两者存到一个矩阵中
both = np.matrix([X, X2])
#同时做出我们通过Garch(1,1)模拟生成的序列X及正态序列的收益时序图
plt.figure(figsize=(16, 9)) 
plt.plot(both.T, alpha=.7);
plt.axhline(X2.std(), color='yellow', linestyle='--')
plt.axhline(-X2.std(), color='yellow', linestyle='--')
plt.axhline(3*X2.std(), color='red', linestyle='--')
plt.axhline(-3*X2.std(), color='red', linestyle='--')
plt.xlabel('time')
plt.ylabel('sigma')
plt.legend(["GARCH","Normal"])
Out[59]:
<matplotlib.legend.Legend at 0x7f34d6dea0b8>

上面我们用蒙特卡洛模拟,模拟出了(收益率)序列$a_t$($r_t$),下来我们就要假装我们不知道一个序列是不是满足Garch(1,1),来进行拟合,如原文所提到的,先进行Arch检验:Arch效应的检验有两种方法,第一种是Ljung-Box统计量用于序列$a_t^2$,即检验残差是否相关。第二种是Engel(1982)提出的拉格朗日乘子检验。该检验等价于线性回归中用F统计量检验$\alpha_i=0(i=1,2…m)$,检验方程为$a_t^2= \alpha_0+ \alpha_1a_{t-1}^2+…+\alpha_ma_{t-m}^2+e_t,t=m+1,…T$ Quantopian原文中使用的应该是Ljung-Box检验的变种(因为我没找到完全吻合的- – !),P值越小越显著,Arch效应越显著,可以对波动率建模

In [60]:
X, _ = simulate_GARCH(1100, a0, a1, b1, sigma1)
X = X[100:] # Drop burn in

p = 20

# Drop the first 20 so we have a lag of p's
Y2 = (X**2)[p:]
X2 = np.ndarray((980, p))
for i in range(p, 1000):
    X2[i - p, :] = np.asarray((X**2)[i-p:i])[::-1]
con=np.empty((980,1))
for i in range(980):
    con[i,:]=1
X2=np.hstack((con,X2))
model = sm.regression.linear_model.OLS(Y2, X2)#原文中并没有拟合常数项,但是公式中显然用到了,到底拟不拟合计入F呢?
#大家还是按上面介绍的第二种方法做稳妥,这里我添加了,否则就跑偏了
model = model.fit()
theta = np.matrix(model.params)
omega = np.matrix(model.cov_HC0)
F = np.asscalar(theta * np.linalg.inv(omega) * theta.T)

#print np.asarray(theta.T).shape
print(theta)
plt.plot(range(21), np.asarray(theta.T))
plt.xlabel('Lag Amount')
plt.ylabel('Estimated Coefficient for Lagged Datapoint')

print ('F ='+ str(F))

chi2dist = scipy.stats.chi2(p)
pvalue = 1-chi2dist.cdf(F)
print ('p-value = ' + str(pvalue))#Pvalue is small,Arch condition is met.

# Finally let's look at the significance of each a_p as measured by the standard deviations away from 0
[[  5.00398970e+00   7.91909523e-02   8.16054764e-02   1.02422380e-01
    3.48122490e-02   8.52583813e-03  -3.48065747e-04   3.11362759e-02
    9.14522132e-02  -2.17775768e-02   2.22012631e-02  -3.66738314e-02
   -1.06271844e-02   6.83422118e-02   2.22335623e-02   8.14950773e-02
   -1.10637908e-03  -4.52856591e-03  -5.49474538e-04  -3.03768472e-02
    3.76003225e-02]]
F =572.2960640396009
p-value = 0.0

接下来就是用MLE和GMM方法来估计波动率方程$\sigma_t^2=\alpha_0+\alpha_1a_{t-1}^2+ b_1a_{t-2}^2$参数值,同理,我们先模拟出波动率

In [61]:
X, _ = simulate_GARCH(10000, a0, a1, b1, sigma1)
X = X[1000:] # Drop burn in
In [62]:
def compute_squared_sigmas(X, initial_sigma, theta):#生成sigma序列
    
    a0 = theta[0]
    a1 = theta[1]
    b1 = theta[2]
    
    T = len(X)
    sigma2 = np.ndarray(T)
    
    sigma2[0] = initial_sigma ** 2
    
    for t in range(1, T):
        # Here's where we apply the equation
        sigma2[t] = a0 + a1 * X[t-1]**2 + b1 * sigma2[t-1]
    
    return sigma2
In [63]:
plt.plot(range(len(X)), compute_squared_sigmas(X, np.sqrt(np.mean(X**2)), (1, 0.5, 0.5)))#这里的a0,a1,b1并不和我们初始一样,
#就是为了造成误差,否则就完全拟合了
plt.xlabel('Time')
plt.ylabel('Sigma');
In [64]:
def negative_log_likelihood(X, theta):
    
    T = len(X)
    
    # Estimate initial sigma squared
    initial_sigma = np.sqrt(np.mean(X ** 2))
    
    # Generate the squared sigma values
    sigma2 = compute_squared_sigmas(X, initial_sigma, theta)
    
    # Now actually compute
    return -sum(
        [-np.log(np.sqrt(2.0 * np.pi)) -
        (X[t] ** 2) / (2.0 * sigma2[t]) -
        0.5 * np.log(sigma2[t]) for
         t in range(T)]
    )
In [65]:
# 进行MLE估计,求出估计值
objective = partial(negative_log_likelihood, X)
# Define the constraints for our minimizer
def constraint1(theta):
    return np.array([1 - (theta[1] + theta[2])])

def constraint2(theta):
    return np.array([theta[1]])

def constraint3(theta):
    return np.array([theta[2]])

cons = ({'type': 'ineq', 'fun': constraint1},
        {'type': 'ineq', 'fun': constraint2},
        {'type': 'ineq', 'fun': constraint3})

# Actually do the minimization
result = scipy.optimize.minimize(objective, (1, 0.5, 0.5),
                        method='SLSQP',
                        constraints = cons)
theta_mle = result.x
print ('theta MLE: ' + str(theta_mle))
theta MLE: [ 1.08330365  0.09764981  0.79638653]

估计值已经出来,那估计的效果怎么样呢,原文用两个效果衡量,并且都是绝对比较:
1.How fat are the tails of the residuals 也就是收益率$x_t$,或者$a_t$.
2.How normal are the residuals under the Jarque-Bera normality test(从下面的P-value很大可以看出,无法拒绝正太性原假设)
但是绝对比较并不符合统计要求,完整的做法是对$\{a_t\}$做Ljung-Box,对$\{a_t^2\}$做Ljung-Box以及保证估计参数的Pvalue很小

In [66]:
def check_theta_estimate(X, theta_estimate):
    initial_sigma = np.sqrt(np.mean(X ** 2))
    sigma = np.sqrt(compute_squared_sigmas(X, initial_sigma, theta_estimate))
    epsilon = X / sigma
    print ('Tails table')
    print (compare_tails_to_normal(epsilon / np.std(epsilon)))
    print ( '')
    
    _, pvalue, _, _ = jarque_bera(epsilon)
    print ('Jarque-Bera probability normal: ' + str(pvalue))
    
check_theta_estimate(X, theta_mle)
Tails table
[[  1.54888889e-01   2.22222222e-02   2.00000000e-03   1.11111111e-04]
 [  1.58655254e-01   2.27501319e-02   1.34989803e-03   3.16712418e-05]]

Jarque-Bera probability normal: 0.218225625935

除了用MLE估计,还可用GMM估计,具体原理见原文(思想是迭代)

In [67]:
#定义求阶矩公式
def standardized_moment(x, mu, sigma, n):
    return ((x - mu) ** n) / (sigma ** n)
In [68]:
def gmm_objective(X, W, theta):
    # Compute the residuals for X and theta
    initial_sigma = np.sqrt(np.mean(X ** 2))
    sigma = np.sqrt(compute_squared_sigmas(X, initial_sigma, theta))
    e = X / sigma
    
    # Compute the mean moments
    m1 = np.mean(e)
    m2 = np.mean(e ** 2) - 1
    m3 = np.mean(standardized_moment(e, np.mean(e), np.std(e), 3))
    m4 = np.mean(standardized_moment(e, np.mean(e), np.std(e), 4) - 3)
    
    G = np.matrix([m1, m2, m3, m4]).T
    
    return np.asscalar(G.T * W * G)

def gmm_variance(X, theta):
    # Compute the residuals for X and theta    
    initial_sigma = np.sqrt(np.mean(X ** 2))
    sigma = np.sqrt(compute_squared_sigmas(X, initial_sigma, theta))
    e = X / sigma

    # Compute the squared moments
    m1 = e ** 2
    m2 = (e ** 2 - 1) ** 2
    m3 = standardized_moment(e, np.mean(e), np.std(e), 3) ** 2
    m4 = (standardized_moment(e, np.mean(e), np.std(e), 4) - 3) ** 2
    
    # Compute the covariance matrix g * g'
    T = len(X)
    s = np.ndarray((4, 1))
    for t in range(T):
        G = np.matrix([m1[t], m2[t], m3[t], m4[t]]).T
        s = s + G * G.T
    
    return s / T
In [69]:
# Initialize GMM parameters
W = np.identity(4)
gmm_iterations = 10

# First guess
theta_gmm_estimate = theta_mle

# Perform iterated GMM
for i in range(gmm_iterations):
    # Estimate new theta
    objective = partial(gmm_objective, X, W)
    result = scipy.optimize.minimize(objective, theta_gmm_estimate, constraints=cons)
    theta_gmm_estimate = result.x
    print ('Iteration ' + str(i) + ' theta: ' + str(theta_gmm_estimate))
    
    # Recompute W
    W = np.linalg.inv(gmm_variance(X, theta_gmm_estimate))
    

check_theta_estimate(X, theta_gmm_estimate)
Iteration 0 theta: [ 0.89961561  0.08336233  0.82815649]
Iteration 1 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 2 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 3 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 4 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 5 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 6 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 7 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 8 theta: [ 0.89962901  0.08348014  0.828286  ]
Iteration 9 theta: [ 0.89962901  0.08348014  0.828286  ]
Tails table
[[  1.55111111e-01   2.25555556e-02   2.22222222e-03   1.11111111e-04]
 [  1.58655254e-01   2.27501319e-02   1.34989803e-03   3.16712418e-05]]

Jarque-Bera probability normal: 0.231900700185

既然我们已经估计出了模型参数,那么开始预测(装逼)吧!

In [70]:
sigma_hats = np.sqrt(compute_squared_sigmas(X, np.sqrt(np.mean(X**2)), theta_mle))
initial_sigma = sigma_hats[-1]
initial_sigma#拿出最后一个值,要从此值开始我们的预测之旅
Out[70]:
2.6194141066428118
In [71]:
a0_estimate = theta_gmm_estimate[0]
a1_estimate = theta_gmm_estimate[1]
b1_estimate = theta_gmm_estimate[2]
#这里用了GMM的估计值,差不了多少
X_forecast, sigma_forecast = simulate_GARCH(100, a0_estimate, a1_estimate, b1_estimate, initial_sigma)#这样就求出了预计收益率和预计波动率
In [72]:
plt.plot(range(-100, 0), X[-100:], 'b-')
plt.plot(range(-100, 0), sigma_hats[-100:], 'r-')
plt.plot(range(0, 100), X_forecast, 'b--')
plt.plot(range(0, 100), sigma_forecast, 'r--')
plt.xlabel('Time')
plt.legend(['X', 'sigma']);
In [73]:
plt.plot(range(-100, 0), X[-100:], 'b-')
plt.plot(range(-100, 0), sigma_hats[-100:], 'r-')
plt.xlabel('Time')
plt.legend(['X', 'sigma'])
Out[73]:
<matplotlib.legend.Legend at 0x7f34de445240>

因为只是模拟了一次,显然不具备代表性,所以文章模拟了100次,并作出每次的拟合值及极值序列

In [74]:
max_X = [-np.inf]
min_X = [np.inf]
for i in range(100):
    X_forecast, sigma_forecast = simulate_GARCH(100, a0_estimate, a1_estimate, b1_estimate, initial_sigma)
    if max(X_forecast) > max(max_X):
        max_X = X_forecast
    elif min(X_forecast) < min(max_X):
        min_X = X_forecast
    plt.plot(range(0, 100), X_forecast, 'b--', alpha=0.05)
    plt.plot(range(0, 100), sigma_forecast, 'r--', alpha=0.05)

# Draw the most extreme X values specially
plt.plot(range(0, 100), max_X, 'g--', alpha=1.0)
plt.plot(range(0, 100), min_X, 'g--', alpha=1.0);

总结:时间序列波动率建模的顺序就是:
1.建立原文省略的均值方程(甚至可以用AR,MA,ARMA建立)
2.检验ARCH效应,如果显著,则对均值方程(收益率方程)和波动率方程建模
3.利用MLE(GMM)拟合
4.检查拟合效果