那些年，我们追过的时间序列【2】随机波动率模型、QMLE ...

1、最简单的一种波动率模型： GARCH 模型

资产是存在波动率的，对于做衍生品交易尤其是期权交易尤其重要。因此，需要把资产的波动率考虑到模型中。Engle（1982）把波动率看成条件标准差，认为误差项和自变量之间存在着相关性，相关性用波动率表示，写为。后来，Bollerslev (1986) 扩展了ARCH模型提出了GARCH模型，形式如下：

GARCH 是比较简单的一种波动率方程，可以帮助我们判断一个资产的波动是多少。下面我们来看看怎么应用。我取了2021年至今英伟达公司的收盘价取对数收益率，如下图所示：

把对数收益放入GARCH（1，1）模型中，取条件分布为正太分布，得到下面的结果：
library(fGarch, quietly = TRUE)
mod1 <- garchFit(~ 1 + garch(1,1), data=ts(lgprice), trace=FALSE)
summary(mod1)

Title:
GARCH Modelling

Call:
garchFit(formula = ~1 + garch(1, 1), data = ts(lgprice), trace = FALSE)

Mean and Variance Equation:
data ~ 1 + garch(1, 1)
<environment: 0x0000019c553e5bd8>
[data = ts(lgprice)]

Conditional Distribution:
norm

Coefficient(s):
      mu    omega    alpha1    beta1
3.6943e-03  4.0443e-05  4.8032e-02  9.0033e-01

Std. Errors:
based on Hessian

Error Analysis:
      Estimate  Std. Error  t value Pr(>|t|)
mu    3.694e-03 1.721e-03 2.146 0.0319 *
omega  4.044e-05 4.833e-05 0.837 0.4027
alpha1 4.803e-02 3.451e-02 1.392 0.1640
beta1  9.003e-01 9.036e-02 9.964 <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log Likelihood:
546.7847 normalized:  2.169781

Standardised Residuals Tests:
                              Statistic p-Value
Jarque-Bera Test R Chi^2  22.53684  1.276991e-05
Shapiro-Wilk Test  R W    0.9860416 0.01469585
Ljung-Box Test    R Q(10)  6.850937  0.7394385
Ljung-Box Test    R Q(15)  8.478535  0.9031694
Ljung-Box Test    R Q(20)  16.38425  0.6925183
Ljung-Box Test    R^2  Q(10)  14.02068  0.1720506
Ljung-Box Test    R^2  Q(15)  18.39769  0.2423351
Ljung-Box Test    R^2  Q(20)  20.56725  0.4229852
LM Arch Test    R TR^2 15.52993  0.2137268

Information Criterion Statistics:
   AIC    BIC    SIC    HQIC
-4.307815 -4.251793 -4.308309 -4.285273 标准化残差和白噪声检验通过，把算出来的数据带回到GARCH中：

由此，拟合出来的波动率图形如下：
vola <- volatility(mod1)
plot(ts(vola, start=start(ts(lgprice)), frequency=frequency(ts(lgprice))),
   xlab=&#34;day&#34;, ylab=&#34;Voltality&#34;)
abline(h=sd(ts(lgprice)), col=&#34;green&#34;)

用该模型的参数做预测，取两倍的方差作为范围，几乎可以覆盖95%的数据，结果如下：

GARCH还可以衍生出各种版本，比如 EGARCH、IGARCH、TGARCH、GARCH-M，等等。但是总的来说，都是对于波动率进行一些人为的限定以及靠过去的数据计算出来一些经验值，主流的金融论文都在用随机波动率模型来建模了（Stochastic Volatility Model）下面再来看看随机波动率模型是怎么做的
2、随机波动率模型（Stochastic Volatility Model）

当状态空间变为线性非高斯白噪声时，卡尔曼滤波器或者MLE就无法正常估计了。需要用到Quasi-MLE或者GMM估计。下面以随机波动率模型为例子（Stochastic Volatility Model）
假设为资产的收益率，是对应的波动率，是状态变量，服从iid，

2.1、Quasi Maximum Likelihood Estimation (QMLE) 估计：

按照 Ruiz （1994）的方法，对公式两边取 log，我们可以把上面的方程转化为下面的状态空间形式：

利用前一篇文章写的卡尔曼滤波器，就可以推导出参数了。（我认为在这里，QML和卡尔曼滤波+MLE是等价的，只是不同的叫法。）
代码如下：
y = log(resi^2)
num = length(y)
# Initial Parameters
phi0 = 0; phi1 =.95; sQ =.2; alpha = mean(y)
sR0 = 1; mu1 = -3; sR1 =2
init.par = c(phi0, phi1, sQ, alpha, sR0, mu1, sR1)
# Innovations Likelihood
Linn = function(para){
  phi0 = para[1]; phi1 = para[2]; sQ = para[3]; alpha = para[4]
  sR0 = para[5]; mu1 = para[6]; sR1 = para[7]
  sv = SVfilter(num, y, phi0, phi1, sQ, alpha, sR0, mu1, sR1)
  return(sv$like) }
# Estimation
(est = optim(init.par, Linn, NULL, method=&#39;BFGS&#39;, hessian=TRUE,
         control=list(trace=1,REPORT=1)))
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimates=est$par, SE)
rownames(u)=c(&#39;phi0&#39;,&#39;phi1&#39;,&#39;sQ&#39;,&#39;alpha&#39;,&#39;sigv0&#39;,&#39;mu1&#39;,&#39;sigv1&#39;); u
# Graphics (need filters at the estimated parameters)
phi0 = est$par[1]; phi1 = est$par[2]; sQ = est$par[3]; alpha = est$par[4]
sR0 = est$par[5]; mu1 = est$par[6]; sR1 = est$par[7]
sv = SVfilter(num,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1)参数结果如下

estimates	se
phi0	-1.082	1.486
phi1	0.517	0.414
sQ	0.519	0.417
alpha	1.832	1.893
sigv0	0.957	0.280
mu1	-2.362	0.524
sigv1	2.787	0.307

SVfilter = function(num,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1){
# Initialize
y=as.matrix(y)
Q=sQ^2
R0=sR0^2
R1=sR1^2
# initial values
xp=matrix(0,num,1)       # = h_t+1^t
Pp=matrix(0,num,1)       # = P_t+1^t
xp[1]= 0                # = h_1^0 | E(h_0) = mu 0 = 0    initial conditions
Pp[1]= phi1^2 + Q       # = P_1^0 | var(h_0) = Sigma0 = 1
pi1=.5                   # initial mix probs
pi0=.5
pit1=.5
pit0=.5
like=0                   # -log(likelihood)
#
for (i in 2:num){
  sig1 = Pp[i-1]+R1    #innov var
  sig0 = Pp[i-1]+R0
  k1 = phi1*Pp[i-1]/sig1
  k0 = phi1*Pp[i-1]/sig0
  e1 = y-xp[i-1]-mu1-alpha
  e0 = y-xp[i-1]-alpha
  den1 = (1/sqrt(sig1))*exp(-.5*e1^2/sig1)
  den0 = (1/sqrt(sig0))*exp(-.5*e0^2/sig0)
  denom = pi1*den1+pi0*den0
  pit1 = pi1*den1/denom
  pit0 = pi0*den0/denom
#
  xp = phi0 + phi1*xp[i-1] + pit0*k0*e0 + pit1*k1*e1
  Pp = (phi1^2)*Pp[i-1] + Q - pit0*(k0^2)*sig0 - pit1*(k1^2)*sig1
#
  like = like - 0.5*log(pi1*den1 + pi0*den0)
}
list(xp=xp,Pp=Pp,like=like)
}拟合的密度函数图像如下：
# densities plot (f is chi-sq, fm is fitted mixture)
x = seq(-15,6,by=.01)
f = exp(-.5*(exp(x)-x))/(sqrt(2*pi))
f0 = exp(-.5*(x^2)/sR0^2)/(sR0*sqrt(2*pi))
f1 = exp(-.5*(x-mu1)^2/sR1^2)/(sR1*sqrt(2*pi))
fm = (f0+f1)/2
plot(x, f, type=&#39;l&#39;); lines(x, fm, lty=2, lwd=2)

2.2、GMM 估计：(先鸽了，嘿嘿)

虽然 MLE可以被归纳到 GMM 框架下，但是当似然函数变得复杂之后，仿佛 GMM 也就是用来计算矩条件，之后还是要算似然函数的。暂时还没有参透他们之间的关系。要是有谁懂的话，麻烦留言告诉我其中的玄机。
3、一篇实证论文：加入随机波动率的利率期限结构模型

Joslin and Le (2021) 在利率期限结构模型上加入了随机波动率，构建了下面的模型：

他们以此为基础构建了一个以债券收益率为隐含因子的无风险套利模型：

通过QML估计，可以得到如下结果：

真实世界测度下的参数结果：

他们发现由于存在着随机波动率，导致最终的利率收益是由随机波动影响的而不是由隐含因子影响的。

Reference
Joslin, S., & Le, A. (2021). Interest rate volatility and no-arbitrage affine term structure models. Management Science.
Ruiz, E. (1994). Quasi-maximum likelihood estimation of stochastic volatility models. Journal of econometrics,63(1), 289-306.
Tsay, R. S. (2014). An introduction to analysis of financial data with R. John Wiley & Sons.