28 状态空间模型

上一章的局部水平模型是线性高斯状态空间模型的一个简单特例。 本章给出状态空间模型， 举例说明这种模型能够表示的其它模型， 如ARIMA模型，结构时间序列模型， 时变回归模型，有自相关误差的回归模型， 随机波动率模型等， 并给出滤波、平滑、预报公式和参数估计方法。

参考：

28.1 状态空间模型的公式

28.1.1 线性高斯状态空间模型

状态空间模型有许多不同的表达形式， 按照(J. Durbin and Koopman 2012)的公式， 线性高斯模型为： \[\begin{align} \boldsymbol y_t =& Z_t \boldsymbol \alpha_t + \boldsymbol \varepsilon_t, \ \boldsymbol \varepsilon_t \sim \text{N}(0, H_t), \tag{28.1} \\ \boldsymbol \alpha_{t+1} =& T_t \boldsymbol \alpha_t + R_t \boldsymbol \eta_t, \ \boldsymbol \eta_t \sim \text{N}(0, Q_t), \tag{28.2} \end{align}\] 其中\(\boldsymbol y_t\)是\(t\)时刻的观测值， 为\(p \times 1\)向量； \(\boldsymbol \alpha_t\)时\(t\)时刻系统的状态， 是不可观测的\(m\times 1\)随机向量， 第一个方程称为观测方程，第二个方程称为状态方程。 \(\{\boldsymbol \varepsilon_t\}\)和\(\{\boldsymbol \eta_t \}\)相互独立， 都是独立同分布向量白噪声列， \(\boldsymbol \varepsilon_t\)为\(p\times 1\)随机向量， \(\boldsymbol \eta_t\)为\(r \times 1\)随机向量，\(r \leq m\)。 设各矩阵\(Z_t, T_t, R_t, H_t, Q_t\)已知， \(Z_t\)和\(T_{t-1}\)允许依赖于\(\boldsymbol y_1, \dots, \boldsymbol y_{t-1}\)， 初始状态\(\boldsymbol\alpha_1\)服从\(N(\boldsymbol a_1, P_1)\)， 设\(\boldsymbol a_1, P_1\)已知， \(\boldsymbol\alpha_1\)与\(\{\boldsymbol \varepsilon_t\}\)和\(\{\boldsymbol \eta_t \}\)独立。 当参数未知时， 设\(\boldsymbol\psi\)为未知参数， 矩阵\(Z_t, T_t, R_t, H_t, Q_t\)可以依赖于未知参数\(\boldsymbol\psi\)。

模型中的\(R_t\)常常是单位阵，\(r=m\)， 有些教材的模型就没有\(R_t\)这一项。 包含\(R_t\)的好处是， \(R_t\)常常是单位阵\(I_m\)的某些列组成的一个\(m \times r\)矩阵， 称为选择矩阵， 这允许某些状态分量对应的方程误差为0， 同时\(\boldsymbol \eta_t\)的方差阵\(Q_t\)还可以是满秩的\(r \times r\)正定阵， 如果没有\(R_t\)矩阵\(Q_t\)就可能不满秩。 如果\(R_t\)是一般的\(m \times r\)矩阵， 关于状态空间模型的大部分结论仍成立。

28.1.2 推广的状态空间模型

可以将线性高斯的状态空间模型， 推广到状态方程仍为线性高斯形式， 而观测方程的分布为非高斯分布， 或者观测方程中观测变量与状态变量的关系非线性， 更进一步可以推广到状态方程的关系也非线性， 分布为非高斯分布。

较一般的非线性、非高斯状态空间模型形式为： \[\begin{aligned} \boldsymbol y_t \sim& f_t(\boldsymbol\alpha_t; \boldsymbol\beta), \\ \boldsymbol\alpha_{t+1} \sim& g_t(\boldsymbol\alpha_t; \boldsymbol\theta), \end{aligned}\] 其中\(f_t(\cdot)\), \(g(\cdot)\)是密度函数或概率质量函数， \(\boldsymbol\beta\)和\(\boldsymbol\theta\)是超参数。 这样的模型一般需要用MCMC、序贯重要抽样等随机模拟方法进行滤波、平滑和估计。

28.2 状态空间模型示例

下面先给出多种时间序列模型的状态空间模型表示。

28.2.1 一元结构时间序列模型

前一章的局部水平模型是一元结构时间序列模型的特例。 在局部水平模型中增加一个斜率\(\nu_t\)项，变成 \[\begin{aligned} y_t =& \mu_t + e_t, \\ \mu_{t+1} =& \mu_t + \nu_t + \xi_t, \\ \nu_{t+1} =& \nu_t + \zeta_t, \end{aligned}\] 写成状态空间模型，为 \[\begin{aligned} y_t =& (1\ 0) \begin{pmatrix} \mu_t \\ \nu_t \end{pmatrix} + e_t, \\ \begin{pmatrix} \mu_{t+1} \\ \nu_{t+1} \end{pmatrix} =& \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} \mu_t \\ \nu_t \end{pmatrix} + \begin{pmatrix} \xi_t \\ \zeta_t \end{pmatrix} . \end{aligned}\]

再增加一个季节项\(\gamma_t\)， 季节项有多种模型可选， 比如\(\gamma_t = -\sum_{j=1}^{s-1} \gamma_{t-j} + \omega_t\)， 或者同一季度的随机游动， 或者三角多项式形式的表达式。 这里用 \[\begin{aligned} y_t =& \mu_t + \gamma_t + e_t, \\ \mu_{t+1} =& \mu_t + \nu_t + \xi_t, \\ \nu_{t+1} =& \nu_t + \zeta_t, \\ \gamma_{t+1} =& -\sum_{j=1}^{s-1} \gamma_{t+1-j} + \omega_t \end{aligned}\] 写成SSM形式， \[\begin{aligned} \boldsymbol\alpha_t =& (\mu_t, \nu_t, \gamma_t, \dots, \gamma_{t-s+2})^T, \\ y_t =& (1,0,1,0,\dots,0) \boldsymbol\alpha_t + e_t, \\ \boldsymbol\alpha_{t+1} =& T_t \boldsymbol\alpha_{t} + R_t \boldsymbol\eta_t . \end{aligned}\]

例如，当\(s=4\)（季度数据）时， \[\begin{aligned} \boldsymbol\alpha_t =& (\mu_t, \nu_t, \gamma_t, \gamma_{t-1}, \gamma_{t-2})^T, \\ y_t =& (1,0,1,0,0) \boldsymbol\alpha_t + e_t, \\ \begin{pmatrix} \mu_{t+1} \\ \nu_{t+1} \\ \gamma_{t+1} \\ \gamma_t \\ \gamma_{t-1} \end{pmatrix} =& \begin{pmatrix} 1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & -1 & -1 & -1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ \end{pmatrix} \begin{pmatrix} \mu_{t} \\ \nu_{t} \\ \gamma_{t} \\ \gamma_{t-1} \\ \gamma_{t-2} \end{pmatrix} \\ & + \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \xi_t \\ \zeta_t \\ \omega_t, \end{pmatrix}, \quad \begin{pmatrix} \xi_t \\ \zeta_t \\ \omega_t, \end{pmatrix} \sim \text{N} \left( \boldsymbol 0, \text{diag}(\sigma_{\xi}^2, \sigma_{\zeta}^2, \sigma_{\omega}^2) \right) . \end{aligned}\]

模型中可以加入商业周期项\(c_t\)， 可以增加若干个回归自变量，可以增加干预变量如： \(w_t = I_{\{t \geq \tau \} }\)。

例27.1 考虑Alcoa股票日现实波动率数据的局部水平模型， 时间期间为2003-01-02到2004-05-07， 共340个观测。

读入数据：

da <- readr::read_table2("aa-3rv.txt", 
  col_names=FALSE)
## Warning: `read_table2()` was deprecated in readr 2.0.0.
## Please use `read_table()` instead.
## 
## -- Column specification --------------------------------------------------------
## cols(
##   X1 = col_double(),
##   X2 = col_double(),
##   X3 = col_double()
## )
ts.alcoa <- ts(log(da[[2]]))

使用局部水平模型，状态空间模型为： \[\begin{aligned} y_t =& \mu_t + e_t, \\ \mu_{t+1} =& \mu_t + \eta_t . \end{aligned}\]

使用statespacer包进行估计， 这个包直接支持对结构时间序列模型进行简化的模型设定：

library(statespacer)
ssr1 <- statespacer(
  y = cbind(as.vector(ts.alcoa)),
  local_level_ind = TRUE,
  initial = rep(0.5*log(var(ts.alcoa)), 2),
  verbose = TRUE)
## Starting the optimisation procedure at: 2022-01-23 17:10:39
## initial  value 1.022439 
## iter  10 value 0.765458
## iter  20 value 0.764395
## final  value 0.764395 
## converged
## Finished the optimisation procedure at: 2022-01-23 17:10:39
## Time difference of 0.0313341617584229 secs
c("Var_obs"=ssr1$system_matrices$H$H, 
  "Var_level"=ssr1$system_matrices$Q$level)
##     Var_obs   Var_level 
## 0.230623632 0.005404681

超参数估计为 \[ \sigma_e^2 = 0.2306, \ \sigma_{\eta}^2 = 0.005405 . \]

做原始序列与平滑得到的趋势（水平值）的时间序列图， 包括95%预测区间：

smsd <- sqrt(ssr1$smoothed$V[1,1,])
plot(ts.alcoa, ylim=c(-1.5, 3))
lines(as.vector(time(ts.alcoa)), 
  ssr1$smoothed$level, col="green")
lines(as.vector(time(ts.alcoa)), 
  ssr1$smoothed$level - 1.96*smsd,
  lty=2, col="cyan")
lines(as.vector(time(ts.alcoa)), 
  ssr1$smoothed$level + 1.96*smsd,
  lty=2, col="cyan")
legend("bottomleft", lty=c(1,1,2), 
  col=c("black", "green", "cyan"),
  legend=c("Obs", "Smoothed", "95% CL"))

下面改用KFAS包进行估计。 这个包也支持对结构时间序列模型进行简化的模型设定。

library(KFAS)
## 模型设定
kfas01a <- SSModel(
  ts.alcoa ~ SSMtrend(1, Q=list(matrix(NA))), 
  H=matrix(NA))
## 超参数估计
kfas01b <- fitSSM(
  kfas01a, rep(log(var(ts.alcoa)), 2), method="BFGS")
c("Var_obs"=c(kfas01b$model$H), 
  "Var_level"=c(kfas01b$model$Q))
##     Var_obs   Var_level 
## 0.230651733 0.005404091
## 用估计的超参数进行滤波平滑
kfas01c <- KFS(kfas01b$model)

超参数估计结果与statespacer基本一致。 计算两个包的平滑结果的最大差距：

c("水平值平滑差距"=max(abs(c(ssr1$smoothed$level) - c(kfas01c$muhat))),
  "水平值平滑方差差距"=max(abs(c(ssr1$smoothed$V) - c(kfas01c$V_mu))))
##     水平值平滑差距 水平值平滑方差差距 
##       3.475598e-05       9.103878e-07

可见statespacer和KFAS两个包的平滑结果是一致的。

○○○○○

28.2.2 ARMA和ARIMA模型

考虑ARMA(\(p,q\))模型， 取\(r=\max(p,q+1)\)， 则模型可以写成 \[ y_t = \sum_{j=1}^r \phi_j y_{t-j} + \zeta_t + \sum_{j=1}^{r-q} \theta_j \zeta_{t-j}, \] 其中某些系数可以为0。 将此模型写成状态空间模型形式， 有多种不同形式， 这里给出一种形式比较复杂但计算方便的表达形式： \[\begin{aligned} \boldsymbol\alpha_t =& \begin{pmatrix} y_t & \\ \phi_2 y_{t-1} + \dots + \phi_r y_{t - (r-1)} & + \theta_1 \zeta_t + \dots + \theta_{r-1} \zeta_{t-(r-2)} \\ \phi_3 y_{t-1} + \dots + \phi_r y_{t - (r-2)} & + \theta_2 \zeta_t + \dots + \theta_{r-1} \zeta_{t-(r-3)} \\ \vdots & \vdots \\ \phi_r y_{t-1} & + \theta_{r-1} \zeta_t \end{pmatrix} \\ y_t =& (1, 0, 0, \dots, 0) \boldsymbol\alpha_t \\ \boldsymbol\alpha_{t+1} =& \begin{pmatrix} \phi_1 & 1 & 0 & \cdots & 0 \\ \phi_2 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & & \ddots & \vdots \\ \phi_{r-1} & 0 & 0 & \cdots & 1 \\ \phi_r & 0 & 0 & \cdots & 0 \end{pmatrix} \boldsymbol\alpha_t + \begin{pmatrix} 1 \\ \theta_1 \\ \theta_2 \\ \vdots \\ \theta_{r-1} \end{pmatrix} \zeta_{t+1} . \end{aligned}\]

对一个ARIMA(2,1,1)， 模型 \[\begin{aligned} y_t^* =& y_t - y_{t-1}, \\ y_t^* =& \phi_1 y_{t-1}^* + \phi_2 y_{t-2}^* + \zeta_t + \theta_1 \zeta_{t-1}, \end{aligned}\] 可以写成 \[\begin{aligned} \boldsymbol\alpha_t =& \begin{pmatrix} y_{t-1} \\ y_t^* \\ \phi_2 y_{t-1}^* + \theta_1 \zeta_t \end{pmatrix} \\ y_t =& (1,1,0) \boldsymbol\alpha_t \\ \boldsymbol\alpha_{t+1} =& \begin{pmatrix} 1 & 1 & 0 \\ 0 & \phi_1 & 1 \\ 0 & \phi_2 & 0 \end{pmatrix} \boldsymbol\alpha_t + \begin{pmatrix} 0 \\ 1 \\ \theta_1 \end{pmatrix} \zeta_{t+1} \end{aligned}\] 观测方程没有观测误差。

更一般的ARIMA以及带有乘性季节部分的ARIMA都可以类似表示。 将ARIMA模型表示成状态空间模型以后， 好处是状态空间模型的一系列工具都可以用在ARIMA的模型推断中， 比如SSM的精确最大似然估计和初始化方法都可以使用。 现在的估计ARIMA模型的软件中许多都是利用状态空间模型形式进行估计和推断。

反过来， 许多状态空间模型的观测值也服从ARIMA模型， 比如局部趋势模型， 带有斜率项的局部趋势模型等， 加入了季节项的结构时间序列模型等。 但是，结构时间序列模型转换成ARIMA模型会丢失一些有可解释含义的信息。

状态空间模型的优点：

  • 状态空间模型更灵活， 能够将随着时间增加的已知的机制变化添加到模型中， 而ARIMA则很难修改；
  • 容易处理缺失值；
  • 很容易增加额外的解释变量， 有回归自变量时允许回归系数为时变系数， 很容易进行日历调整；
  • 预测不需要单独的理论；
  • 不要求平稳或者差分后平稳。

另一方面， ARIMA模型则无法将趋势成分、季节项成分提取出来， 要求差分后平稳， 而经济和金融中许多数据是差分后也不平稳的。

例28.1 考虑太阳黑子数的ARMA建模。 数据为1700-1988的年数据。

数据：

library(datasets)
ts.ss <- window(sunspot.year, start = 1770, end = 1869)
plot(ts.ss)

可以看出，数据有周期性。

用R的arima函数建立AR(3)模型：

ss1 <- arima(
  ts.ss, order=c(3,0,0)
)
ss1
## 
## Call:
## arima(x = ts.ss, order = c(3, 0, 0))
## 
## Coefficients:
##          ar1      ar2     ar3  intercept
##       1.5471  -0.9915  0.2004    48.5119
## s.e.  0.0983   0.1545  0.0990     6.0560
## 
## sigma^2 estimated as 220.2:  log likelihood = -412.94,  aic = 835.88

估计的模型为 \[ (1 - 1.5471B + 0.9915B^2 - 0.2004B^3)(y_t - 48.5119) = e_t, \ e_t \sim \text{WN}(0, 220.2) . \] 模型诊断：

tsdiag(ss1)

从模型诊断看是合适的。

估计的AR(3)中特征多项式的根和复根对应的周期：

rt1 <- polyroot(c(1, -coef(ss1)[1:3])); rt1
## [1] 1.045502+0.808659i 1.045502-0.808659i 2.855749-0.000000i
2*pi/Mod(rt1[1])
## [1] 4.753712

求出的周期与实际观察到的大约11年周期不符。

用statespacer包将AR(3)估计成状态空间模型， 此包直接支持恢复ARIMA参数：

library(datasets)
library(statespacer)
mat.suns <- matrix(window(
  sunspot.year, start = 1770, end = 1869))
fit <- statespacer(
  y = mat.suns,
  H_format = matrix(0), 
  # H是观测方程的误差方差阵
  local_level_ind = TRUE,
  arima_list = list(c(3,0,0)),
  format_level = matrix(0),
  # format_level是结构时间序列模型中水平值$\mu_t$的方程的误差方差阵
  initial = c(0.5*log(var(mat.suns)), 0, 0, 0),
  verbose=TRUE)
## Starting the optimisation procedure at: 2022-01-25 12:00:35
## initial  value 5.022561 
## iter  10 value 4.114363
## final  value 4.111163 
## converged
## Finished the optimisation procedure at: 2022-01-25 12:00:35
## Time difference of 0.0751800537109375 secs

构造ARMA的结果：

arma_coeff <- rbind(
   fit$system_matrices$AR$ARIMA1,
   fit$standard_errors$AR$ARIMA1
)
arma_coeff <- cbind(
   arma_coeff,
   c(fit$smoothed$level[1],
     sqrt(fit$system_matrices$Z_padded$level %*%
          fit$smoothed$V[,,1] %*%
          t(fit$system_matrices$Z_padded$level))
   )
)
rownames(arma_coeff) <- c("coefficient", "std_error")
colnames(arma_coeff) <- c("ar1", "ar2", "ar3", "intercept")
arma_coeff
##                    ar1       ar2       ar3 intercept
## coefficient 1.55976415 -1.005462 0.2129622 48.605905
## std_error   0.09962468  0.155982 0.1003591  6.358039

与arima()函数结果相近但不完全相同。

例28.2 考虑航空乘客数的建模，1949-1961年的月度数据， 用ARIMA\((0,1,1)(0,1,1)_{12}\)模型。 
plot(AirPassengers)

将数据取自然对数， 用arima()建模：

ap01 <- arima(
  log(AirPassengers),
  order=c(0,1,1),
  seasonal=list(order=c(0,1,1), frequency=12))
ap01
## 
## Call:
## arima(x = log(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0, 
##     1, 1), frequency = 12))
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 estimated as 0.001348:  log likelihood = 244.7,  aic = -483.4

模型为： \[ (1 - B) (1 - B^{12}) \ln y_t = (1 - 0.4018B)(1 - 0.5569B^{12}) e_t, \ e_t \sim \text{WN}(0, 0.001348) . \]

利用statespacer计算：

mat.ap <- matrix(log(AirPassengers))

# 模型设定列表
sarima_list <- list(list(
  s = c(12, 1), ar = c(0, 0), i = c(1, 1), ma = c(1, 1) ))

# 拟合模型
fit <- statespacer(
  y = mat.ap,
  H_format = matrix(0),
  sarima_list = sarima_list,
  initial = c(0.5*log(var(diff(mat.ap))), 0, 0),
  verbose = TRUE)
## Starting the optimisation procedure at: 2022-01-23 17:10:40
## initial  value -1.034434 
## final  value -1.616321 
## converged
## Finished the optimisation procedure at: 2022-01-23 17:10:41
## Time difference of 0.584496021270752 secs

提前估计的参数：

arma_coeff <- rbind(
   c(fit$system_matrices$SMA$SARIMA1$S1, fit$system_matrices$SMA$SARIMA1$S12),
   c(fit$standard_errors$SMA$SARIMA1$S1, fit$standard_errors$SMA$SARIMA1$S12)
)
rownames(arma_coeff) <- c("coefficient", "std_error")
colnames(arma_coeff) <- c("ma1 s = 1", "ma1 s = 12")
arma_coeff
##               ma1 s = 1  ma1 s = 12
## coefficient -0.40188859 -0.55694248
## std_error    0.08963614  0.07309788
goodness_fit <- rbind(
   fit$system_matrices$Q$SARIMA1,
   fit$diagnostics$loglik,
   fit$diagnostics$AIC
)
rownames(goodness_fit) <- c("Variance", "Loglikelihood", "AIC")
goodness_fit
##                        [,1]
## Variance        0.001347882
## Loglikelihood 232.750284785
## AIC            -3.010420622

结果与arima()结果相同。

28.2.3 回归模型

28.2.3.1 常系数回归模型

考虑回归模型 \[\begin{aligned} y_t = \boldsymbol X_t \boldsymbol\beta + e_t, \ e_t \sim \text{ iid N}(0, H_t), \ t=1,2,\dots,n , \end{aligned}\] 其中\(\boldsymbol X_t\)为\(1 \times k\)非随机的已知值， \(\boldsymbol\beta\)是\(k \times 1\)未知的回归系数向量。 可以写成状态空间模型： \[\begin{aligned} \boldsymbol \alpha_t =& \boldsymbol\beta, \\ \boldsymbol\alpha_{t+1} =& I \boldsymbol\alpha_t, \\ y_t =& \boldsymbol X_t \boldsymbol\alpha_t + e_t . \end{aligned}\] 状态方程没有误差。

这个回归问题的加权最小二乘解为 \[\begin{aligned} \hat{\boldsymbol\beta} = (\sum_{i=1}^n \boldsymbol X_t^T H_t^{-1} \boldsymbol X_t)^{-1} \sum_{i=1}^n \boldsymbol X_t^T H_t^{-1} y_t . \end{aligned}\] 写成状态空间模型后， 用滤波算法对\(\boldsymbol\alpha_t = \boldsymbol\beta\)进行估计， 等价于对\(\boldsymbol\beta\)进行递推改进估计， 相当于递推的最小二乘估计。

28.2.3.2 变系数回归模型

考虑变系数的回归模型 \[\begin{aligned} y_t = \boldsymbol X_t \boldsymbol\beta_t + e_t, \ e_t \sim \text{ iid N}(0, H_t), \ t=1,2,\dots,n , \end{aligned}\] 其中的随时间而变化的回归系数可以用如下的随机游动建模： \[\begin{aligned} \boldsymbol\beta_{t+1} =& \boldsymbol\beta_t + \boldsymbol\eta_t, \end{aligned}\] 这样，以\(\beta_t\)为状态向量， 就可以将变系数回归模型写成简单的状态空间形式： \[\begin{aligned} y_t =& \boldsymbol X_t \boldsymbol\beta_t + e_t, \\ \boldsymbol\beta_{t+1} =& I \boldsymbol\beta_t + I \boldsymbol\eta_t . \end{aligned}\] 这样的模型可以用标准的状态空间模型工具进行估计和推断。

28.2.3.3 带有ARMA误差的回归模型

设 \[\begin{aligned} y_t = \boldsymbol X_t \boldsymbol\beta + \xi_t, \end{aligned}\] \(\xi_t\)服从ARMA模型， 可以写成状态空间形式， 状态方程为 \[ \boldsymbol\alpha_{t+1} = T \boldsymbol\alpha_t + R \boldsymbol\eta_t, \] 则\(y_t\)可以写成状态空间模型形式， 状态变量为\(\boldsymbol\alpha_t^*\)， 模型为 \[\begin{aligned} \boldsymbol\alpha_t^* =& \begin{pmatrix} \boldsymbol\beta \\ \boldsymbol\alpha_t \end{pmatrix}, \\ y_t =& (\boldsymbol X_t, 1, 0, \dots, 0) \boldsymbol\alpha_t^*, \\ \boldsymbol\alpha_{t+1}^* =& \begin{pmatrix} I_k & 0 \\ 0 & T \end{pmatrix} \boldsymbol\alpha_t^* + \begin{pmatrix} 0 \\ R \end{pmatrix} \eta_t, \end{aligned}\] 可以用状态空间模型的工具进行估计和推断。

28.2.4 Seatbelts数据建模

R的datasets包的Seatbelts数据是一个多元时间序列， 为月度数据， 保存了英国1969年到1984年与道路交通事故有关的数据。 分量含义：

  • DriversKilled: 小汽车交通事故死亡人数。
  • drivers: 小汽车交通事故死亡以及重伤人数。
  • front: 前排驾乘人员死亡和重伤人数。
  • rear: 后排乘客死亡和重伤人数。
  • kms: 已行驶里程。
  • PetroPrice: 汽油价格。
  • VanKilled: 小型厢式载货汽车司机死亡人数。
  • law: 强制安全带法律是否已实施的0-1变量，从1983-01-31起实行。

28.2.4.1 数据预处理

因为取正值的变量建模时比较受限制， 所以将这些变量取对数值。

library(statespacer)
library(datasets)
da <- Seatbelts
for(v in c("drivers", "front", "rear", "PetrolPrice", "kms"))
  da[, v] <- log(da[,v])

28.2.4.2 确定性水平和季节项、有回归自变量的一元模型

以小汽车交通事故死亡、重伤人数（对数值）为因变量。 以汽油价格（对数值）和是否实行了安全带法令为自变量。 模型为 \[\begin{aligned} y_t =& \mu_t + \gamma_t + \beta_1 x_{1t} + \beta_2 x_{2t} + e_t, \\ \mu_{t+1} =& \mu_{t}, \\ \gamma_{t+1} =& -\sum_{j=1}^{11} \gamma_{t+1-j} . \end{aligned}\]

y <- as.matrix(da[,"drivers"])
ssmr1 <- statespacer(
  y = y, # 因变量
  local_level_ind = TRUE, # 是否有局部水平项
  # 设置局部水平项误差方差为零，从而非随机:
  format_level = matrix(0), 
  BSM_vec = 12, # 季节项的周期
  # 设置季节项误差项方差为零，从而非随机:
  format_BSM_list = list(matrix(0)), 
  # 添加的外生自变量，每个因变量分量需要输入一个自变量矩阵，
  # 所以即使是一元因变量也需要用列表
  addvar_list = list(as.matrix(da[, c("PetrolPrice", "law")])),
  method = "BFGS",
  # 唯一的方差参数$\sigma_e^2$的初值，
  # 算法中使用其对数值的2倍以确保取正值
  initial = 0.5 * log(var(as.vector(y))), 
  verbose = TRUE)
## Starting the optimisation procedure at: 2022-01-25 12:00:35
## initial  value -0.443442 
## final  value -0.735372 
## converged
## Finished the optimisation procedure at: 2022-01-25 12:00:35
## Time difference of 0.0653839111328125 secs

观测误差方差\(\sigma_e^2\)估计：

ssmr1$system_matrices$H$H[1,1]
## [1] 0.007402481

非随机的局部水平\(\mu\)的平滑估计：

ssmr1$smoothed$level[1,1]
## [1] 6.401571

第一个自变量，即汽油价格对数值，的回归系数：

ssmr1$smoothed$addvar_coeff[1, 1]
## [1] -0.4521301

第二个自变量，即强制安全带法令是否实行，的回归系数：

ssmr1$smoothed$addvar_coeff[1, 2]
## [1] -0.1971395

下面将原始数据与平滑拟合结果（包括局部水平、季节项与自变量作用）作图：

plot(da[, c("drivers")], type = "l", ylim = c(6.95, 8.1),
     xlab = "year", ylab = "logarithm of drivers")
lines(seq(tsp(da)[1], tsp(da)[2], 1/tsp(da)[3]), 
      ssmr1$smoothed$level
      + ssmr1$smoothed$BSM12
      + ssmr1$smoothed$addvar, 
      type = 'l', col = "red")
legend("topright",
  c("log(drivers)", "Smoothed fit"), 
  lty = c(1,1), lwd=c(2.5, 2.5), col = c("black", "red"))

28.2.4.3 随机水平和季节项、有回归自变量的一元模型

模型为 \[\begin{aligned} y_t =& \mu_t + \gamma_t + \beta_1 x_{1t} + \beta_2 x_{2t} + e_t, \\ \mu_{t+1} =& \mu_{t} + \eta_t, \\ \gamma_{t+1} =& -\sum_{j=1}^{11} \gamma_{t+1-j} + \zeta_t . \end{aligned}\]

和确定性的\(\mu\), \(\gamma_t\)相比， 程序只需将零方差改为非零方差：

y <- as.matrix(da[,"drivers"])
ssmr2 <- statespacer(
  y = y, # 因变量
  local_level_ind = TRUE, # 是否有局部水平项
  # 设置局部水平项误差方差为非零，从而随机:
  format_level = matrix(1), 
  BSM_vec = 12, # 季节项的周期
  # 设置季节项误差项方差为零，从而随机:
  format_BSM_list = list(matrix(1)), 
  # 添加的外生自变量，每个因变量分量需要输入一个自变量矩阵，
  # 所以即使是一元因变量也需要用列表
  addvar_list = list(as.matrix(da[, c("PetrolPrice", "law")])),
  method = "BFGS",
  # 唯一的方差参数$\sigma_e^2$的初值，
  # 算法中使用其对数值的2倍以确保取正值
  initial = log(var(as.vector(y))), 
  verbose = TRUE)
## Warning: Number of initial parameters is less than the required amount of
## parameters (3), recycling the initial parameters the required amount of times.
## Starting the optimisation procedure at: 2022-01-25 14:48:52
## initial  value -0.172423 
## iter  10 value -0.893962
## iter  20 value -0.915493
## iter  30 value -0.915515
## final  value -0.915517 
## converged
## Finished the optimisation procedure at: 2022-01-25 14:48:52
## Time difference of 0.846170902252197 secs

程序中还修改了原来的迭代初值， 原来的迭代初值导致迭代过程出错。 初值需要三个方差参数对应的初值， 这里只给了一个初值， 就重复利用。

估计结果：

## 观测误差方差估计
ssmr2$system_matrices$H$H[1,1]
## [1] 0.003786183
## 局部水平的系统方程的误差方差估计
ssmr2$system_matrices$Q$level[1,1]
## [1] 0.0002676875
## 季节项的系统方程的误差方差估计
ssmr2$system_matrices$Q$BSM12[1,1]
## [1] 1.161438e-06
## 关于汽油价格对数值这一自变量的回归系数估计
ssmr2$smoothed$addvar_coeff[1, 1]
## [1] -0.2913944
## 关于安全带强制法令这一自变量的回归系数估计
ssmr2$smoothed$addvar_coeff[1, 2]
## [1] -0.2377374

原始数据与拟合结果（拟合结果是平滑结果， 利用了平滑的局部水平与自变量线性组合）：

plot(da[, c("drivers")], type = "l", ylim = c(6.95, 8.1),
     xlab = "year", ylab = "logarithm of drivers")
lines(seq(tsp(da)[1], tsp(da)[2], 1/tsp(da)[3]), 
      ssmr2$smoothed$level 
      + ssmr2$smoothed$BSM12
      + ssmr2$smoothed$addvar, 
      type = 'l', col = "red")
legend("topright",
  c("log(drivers)", "smoothed fit"), 
  lty = c(1,1), lwd=c(2.5, 2.5), col = c("black", "red"))

与非随机的水平和季节项的拟合结果比较， 随机的水平和季节项使得拟合结果更贴近原始数据。

实际上，从季节项的系统方程的误差方差\(10^{-6}\)级别可以看出季节项的随机性很小， 下面作平滑估计的季节项图形：

plot(seq(tsp(da)[1], tsp(da)[2], 1/tsp(da)[3]), 
     ssmr2$smoothed$BSM12,
     type = "l", ylim = c(-0.2, 0.3),
     xlab = "year", ylab = "stochastic seasonal")
abline(h = 0)

我们将随机的季节项改为非随机：

y <- as.matrix(da[,"drivers"])
ssmr3 <- statespacer(
  y = y, # 因变量
  local_level_ind = TRUE, # 是否有局部水平项
  # 设置局部水平项误差方差为非零，从而随机:
  format_level = matrix(1), 
  BSM_vec = 12, # 季节项的周期
  # 设置季节项误差项方差为零，从而非随机:
  format_BSM_list = list(matrix(0)), 
  # 添加的外生自变量，每个因变量分量需要输入一个自变量矩阵，
  # 所以即使是一元因变量也需要用列表
  addvar_list = list(as.matrix(da[, c("PetrolPrice", "law")])),
  method = "BFGS",
  # 唯一的方差参数$\sigma_e^2$的初值，
  # 算法中使用其对数值的2倍以确保取正值
  initial = log(var(as.vector(y))), 
  verbose = TRUE)
## Warning: Number of initial parameters is less than the required amount of
## parameters (2), recycling the initial parameters the required amount of times.
## Starting the optimisation procedure at: 2022-01-25 14:48:53
## initial  value -0.577129 
## iter  10 value -0.912856
## final  value -0.912859 
## converged
## Finished the optimisation procedure at: 2022-01-25 14:48:53
## Time difference of 0.255355834960938 secs

估计结果：

## 观测误差方差估计
ssmr3$system_matrices$H$H[1,1]
## [1] 0.004033516
## 局部水平的系统方程的误差方差估计
ssmr3$system_matrices$Q$level[1,1]
## [1] 0.0002681651
## 关于汽油价格对数值这一自变量的回归系数估计
ssmr3$smoothed$addvar_coeff[1, 1]
## [1] -0.2767301
## 关于安全带强制法令这一自变量的回归系数估计
ssmr3$smoothed$addvar_coeff[1, 2]
## [1] -0.2375904

比较三个模型的AIC：

c(ssmr1$diagnostics$AIC, 
  ssmr2$diagnostics$AIC, 
  ssmr3$diagnostics$AIC)
## [1] -1.314494 -1.653950 -1.659052

第三个模型的AIC最小。

28.2.4.4 前后座驾乘人员同时建模的二元模型

考虑前座驾乘人员和后座乘客死亡、重伤人数对数的二元模型， 都以汽油价格对数值、行驶里程对数值和安全带法令为自变量。

y <- log(da[,c("front", "rear")])
Xmat <- as.matrix(da[,c("PetrolPrice", "kms", "law")])
ssmr4 <- statespacer(y = y, # 输入数据是两列矩阵，二元时间序列
  # 观测方程误差项方差阵格式，这里是无限制的二阶方差阵
  H_format = matrix(1,2,2),
  local_level_ind = TRUE,
  # 局部水平的系统方程的误差项方差设置，无限制：
  format_level = matrix(1, 2, 2), 
  BSM_vec = 12,
  # 季节项的系统方程的误差项方差设置，设置为零，表示非随机季节项
  format_BSM_list = list(matrix(0, 2, 2)),
  # 每个因变量所需的自变量矩阵：
  addvar_list = list(Xmat, Xmat),
  method = "BFGS",
  initial = 0.5 * log(diag(var(y))),
  verbose = TRUE)
## Warning: Number of initial parameters is less than the required amount of
## parameters (6), recycling the initial parameters the required amount of times.
## Starting the optimisation procedure at: 2022-01-25 15:32:39
## initial  value -2.682940 
## iter  10 value -4.398713
## iter  20 value -4.978251
## iter  30 value -5.058729
## iter  40 value -5.058840
## iter  40 value -5.058840
## iter  40 value -5.058840
## final  value -5.058840 
## converged
## Finished the optimisation procedure at: 2022-01-25 15:33:01
## Time difference of 21.8911490440369 secs

观测方程误差方差阵估计：

ssmr4$system_matrices$H$H
##              [,1]         [,2]
## [1,] 0.0001216484 0.0001122534
## [2,] 0.0001122534 0.0002452852

局部水平的系统方程的误差方差阵估计：

ssmr4$system_matrices$Q$level
##              [,1]         [,2]
## [1,] 5.827611e-06 5.711507e-06
## [2,] 5.711507e-06 6.543960e-06

6个回归系数及相应的t统计量：

knitr::kable(data.frame(
  x1 = c("front PetrolPrice", "front kms", "front law", 
    "rear PetrolPrice", "rear kms", "rear law"),
  x2 = ssmr4$smoothed$addvar_coeff[1,],
  x3 = ssmr4$smoothed$addvar_coeff[1,] 
  / ssmr4$smoothed$addvar_coeff_se[1,]),
  col.names=c("", "Coef", "t stat"),
  digits=4  )
Coef t stat
front PetrolPrice -0.0451 -2.8250
front kms 0.0204 1.0428
front law -0.0525 -7.0897
rear PetrolPrice -0.0132 -0.6996
rear kms 0.0898 3.6799
rear law -0.0001 -0.0165

与常识一致的是， 后排乘客的伤亡数并不受安全带法令的影响。 汽油价格上升对前排驾乘人员伤亡数有抑制作用。 行驶里程数对前排驾乘人员伤亡数没有显著影响， 但对后排乘客伤亡有正向的影响。

程序中还可以限制某个方差阵的秩， 比如取format_level为 \[ \begin{pmatrix} 1 & 0 \\ 1 & 0 \end{pmatrix} \] 可以限制水平的系统方差的误差方差阵（\(2 \times 2\)方差阵）为秩等于1。

28.2.5 动态Nelson-Siegel模型

考虑如下的关于利率期限结构的动态Nelson-Siegel模型： \[\begin{aligned} y_t(\tau) =& \beta_{1t} + \beta_{2t} \left( \frac{1 - e^{-\lambda \tau}}{\lambda \tau} \right) + \beta_{3t} \left( \frac{1 - e^{-\lambda \tau}}{\lambda \tau} - e^{-\lambda \tau} \right) + e_t(\tau), \\ & \tau=3, 6, 12, 24, 36, 60, 84, 120, \ t=1,2,\dots, n . \end{aligned}\] 其中\(\lambda\)是未知参数， \(\tau\)是贷款期限， \(y_t(\tau)\)是\(t\)时刻期限为\(\tau\)的贷款利息。 设各\(e_t(\tau)\)相互独立， 服从\(\text{N}(0, \sigma_e^2)\)； \(\beta_{1t}\), \(\beta_{2t}\), \(\beta_{3t}\)是时变的回归系数， 设其服从一阶向量自回归模型： \[ \boldsymbol\beta_{t+1} =(I - \Phi) \boldsymbol\mu + \Phi\boldsymbol\beta_t + \boldsymbol\eta_t, \] 其中\(\boldsymbol\beta_t = (\beta_{1t}, \beta_{2t}, \beta_{3t})^T\), \(\boldsymbol\mu\)为\(\boldsymbol\beta_t\)的均值， \(\Phi\)为\(3 \times 3\)的回归系数矩阵， 设\(\boldsymbol\eta_t \sim \text{N}(\boldsymbol 0, \Sigma_{\eta})\), \(\boldsymbol\beta_1 \sim \text{N}(\boldsymbol\mu, P_{\boldsymbol\beta})\)， 初始分布中的\(P_{\boldsymbol\beta}\)也满足平稳性条件 \[ P_{\boldsymbol\beta} - \Phi P_{\boldsymbol\beta} \Phi^T = \Sigma_{\boldsymbol\eta} . \]

因为对每个\(t\)有8个不同的\(\tau\)对应的观测值\(y_t(\tau)\)， 将这8个观测值写成一个\(8 \times 1\)观测值向量\(\boldsymbol y_t\)。 对包含期望值\(\boldsymbol\mu\)的VAR(1)模型， 可转换成如下的状态方程： \[\begin{aligned} \boldsymbol\alpha_{t} =& \begin{pmatrix} \boldsymbol\beta_t \\ \boldsymbol\mu \end{pmatrix}, \\ \begin{pmatrix} \boldsymbol\beta_{t+1} \\ \boldsymbol\mu \end{pmatrix} =& \begin{pmatrix} \Phi & I - \Phi \\ \boldsymbol 0 & I \end{pmatrix} \begin{pmatrix} \boldsymbol\beta_t \\ \boldsymbol\mu \end{pmatrix} + \begin{pmatrix} I \\ \boldsymbol 0 \end{pmatrix} \boldsymbol\eta_t , \ \boldsymbol\eta_t \sim \text{N}(\boldsymbol 0, \Sigma_{\boldsymbol\eta}) . \end{aligned}\]

令\(\boldsymbol\tau = (3, 6, 12, 24, 36, 60, 84, 120)^T\), 观测方程可以写成 \[\begin{aligned} \boldsymbol y_t = \begin{pmatrix} Z^{(1)} & \boldsymbol 0_{8 \times 3} \end{pmatrix} \boldsymbol\alpha_t + \boldsymbol e_t, \ \boldsymbol e_t \sim \text{N}(\boldsymbol 0, \sigma_e^2 I_8), \end{aligned}\] 其中\(Z^{(1)}\)是\(8 \times 3\)矩阵， 第一列元素都等于1， 第二列元素为各\(\frac{1 - e^{-\lambda \tau}}{\lambda \tau}\)值， 第三列元素为各\(\frac{1 - e^{-\lambda \tau}}{\lambda \tau} - e^{-\lambda \tau}\)值。

这样将DNS(动态Nelson-Siegel)模型写成了状态空间形式， 可以用statespacer包进行估计和平滑。 如果不用状态空间模型方法， 可以先对每个\(t\)的8个观测回归得到3个系数在\(t\)时刻的估计值， 然后再对估计的\(\beta_{jt}\)建模， 这样的做法不能充分利用模型结构。

statespacer包对于ARIMA、结构时间序列模型都有一些方便的模型参数规定方式和输出方式， 对于更一般的模型就需要用户指定需要各个矩阵和初始值。

示例所用的数据是YieldCurve扩展包中的FedYieldCurve数据集。 取1984-12-31到2000-12-01的数据子集。 数据是\(192 \times 8\)矩阵， 有192天(对应\(t\))，8个期限（对应\(\tau\)）。

library(YieldCurve, quietly=TRUE, warn.conflicts=FALSE)
data(FedYieldCurve)
str(FedYieldCurve)
## An 'xts' object on 1981-12-31/2012-11-30 containing:
##   Data: num [1:372, 1:8] 12.9 14.3 13.3 13.3 12.7 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:8] "R_3M" "R_6M" "R_1Y" "R_2Y" ...
##   Indexed by objects of class: [Date] TZ: 
##   xts Attributes:  
##  NULL
xts.yc <- FedYieldCurve["1984-12-31/2000-12-01"]
date.yc <- index(xts.yc)
y.yc <- coredata(xts.yc)
plot(xts.yc, multi.panel=FALSE)

需要自己指定模型的各个矩阵和初值。

## 保存自定义设置用的列表
spec <- list()
## 自己规定H矩阵（观测误差的方差阵）
spec$H_spec <- TRUE
## 状态向量维数，3个beta, 3个常数beta均值
spec$state_num <- 6

## 待估参数个数20个，包括：
##   1个 - $\lambda$
##   1个 - $\sigma_e^2$
##   6个 - $\Sigma_{\boldsymbol\eta}$
##   9个 - $\Phi_{3\times 3}$矩阵
##   3个 - $\boldsymbol\mu$
spec$param_num <- 20

## 状态方程中$R$矩阵不起作用，方差结构直接编写在矩阵$Q$中
spec$R <- diag(1, 6, 6) # $I_6$

## 初始状态$\boldsymbol\alpha_1$的方差，元素都取为0，因为平稳？
spec$P_inf <- matrix(0, 6, 6)

## 指定状态向量中不输出到观测的分量，
## 当压缩观测向量为标量时可提高计算效率
spec$state_only <- 4:6

模型要估计的参数（超参数）为\(\lambda\), \(\sigma_e^2\), \((\Sigma_{\boldsymbol\eta})_{3\times 3}\), \(\Phi_{3\times 3}\), \(\boldsymbol\mu_{3\times 1}\)。

对于必须取正值的参数， 将其迭代计算的参数取自然对数后除以2。

对于要求正定的\(\Sigma_{\boldsymbol\eta}\)， 将其做Cholesky分解\(LDL^T\)， 其中\(D\)是对角元素为正值的对角阵， \(L\)是对角元素都等于1的下三角阵， 用\(D\)和\(L\)的6个元素来作为待估参数。 statespacer提供了Cholesky函数将这样编码的参数转换成正定矩阵。

对向量自回归系数矩阵\(\Phi\)， 要求其满足平稳性条件， statespacer包的CoeffARMA()函数用来改造输入的系数矩阵使其满足条件。

因为模型的各个重要矩阵、初始值等都需要编程计算而不是常数值， 所以将自定义设定的sys_mat_fun元素设置为一个自定义函数：

spec$sys_mat_fun <- function(param) {
  ## 输入20个参数，具体意义见前面注释和说明
  
  # 8个期限，即$\tau$值
  maturity <- c(3, 6, 12, 24, 36, 60, 84, 120)
  
  # 从取了对数的参数值回复$\lamabda$的值
  lambda <- exp(2 * param[1])
  
  # \sigma_e^2的值
  sigma2 <- exp(2 * param[2])
  # H是观测方程的误差向量的协方差阵参数。
  # 注意每个$t$对应的观测值是$8\times 1$向量。
  H <- sigma2 * diag(1, 8, 8)

  # 观测方程的Z矩阵，一个$8 \times 6$矩阵，
  # 前三列对应于$\beta_{1t}$, $\beta_{2t}$, $\beta_{3t}$
  # 后三列对应于$\beta$的三个均值，系数为0
  lambda_maturity <- lambda * maturity
  ze <- exp(-lambda_maturity)
  Z <- matrix(1, 8, 3)
  Z[, 2] <- (1 - ze) / lambda_maturity
  Z[, 3] <- Z[, 2] - ze
  # 在观测方程的矩阵Z中增加对应于常数均值的3列
  Z <- cbind(Z, matrix(0, 8, 3)) 

  # 状态方程的误差方差阵（$6 \times 6$）的左上角$3 \times 3$部分
  # 使用6个参数
  Q <- Cholesky(
    param = param[3:8], 
    decompositions = FALSE, 
    format = matrix(1, 3, 3))
  # 在状态方程的误差方差阵中增加对应于常数均值的部分，都是0
  Q <- BlockMatrix(Q, matrix(0, 3, 3)) 
  
  # 从输入的9个矩阵元素生成满足平稳性条件的向量自回归系数矩阵$\Phi$
  Tmat <- CoeffARMA(
    A = array(param[9:17], dim = c(3, 3, 1)),
    variance = Q,
    ar = 1, ma = 0)$ar[,,1]

  # 生成初始状态$\boldsymbol\alpha_1$的方差阵，确保其满足平稳性条件
  T_kronecker <- kronecker(Tmat, Tmat)
  Tinv <- solve(diag(1, dim(T_kronecker)[1], dim(T_kronecker)[2]) - T_kronecker)
  vecQ <- matrix(Q)
  vecPstar <- Tinv %*% vecQ
  P_star <- matrix(vecPstar, dim(Tmat)[1], dim(Tmat)[2])

  
  # $\beta$的均值是最后的3个待估参数，也作为$t=1$时$\beta$的初始分布均值
  a1 <- matrix(param[18:20], 6, 1) 
  
  ## 在6维的状态方程中，设置转移矩阵$T$
  Tmat <- cbind(Tmat, diag(1, 3, 3) - Tmat)
  Tmat <- rbind(Tmat, cbind(matrix(0, 3, 3), diag(1, 3, 3)))
  
  ## 初始状态的方差阵中对应于常数均值部分方差和协方差为0
  P_star <- BlockMatrix(P_star, matrix(0, 3, 3))

  # 函数返回模型所需的所有系统矩阵
  return(list(
    H = H,   # 观测误差随机向量的方差阵
    Z = Z,   # 观测方程的观测矩阵
    Tmat = Tmat, # 状态方程的状态转移矩阵
    Q = Q,   # 状态方程的误差随机向量的方差阵
    a1 = a1, # 初始状态的均值
    P_star = P_star)) # 初始状态的方差阵
}

为了从估计的参数（超参数）提取某些成分， 可以提供一个transform_fun函数：

spec$transform_fun <- function(param) {
  lambda <- exp(2 * param[1])
  sigma2 <- exp(2 * param[2])
  means <- param[18:20]
  return(c(lambda, sigma2, means))
}

下面给20个参数设置算法迭代初值：

initial <- c(
  -1, # $0.5 \ln\lambda$
  -2, # $0.5 \ln\sigma_e^2$
  0, 0, 0, # $\Sigma_{\boldsymbol\eta}$的LDL分解中D的对角元素的0.5倍对数值
  0, 0, 0, # $\Sigma_{\boldsymbol\eta}$的LDL分解中L的元素值
  4, 0, 0, 
  0, 3, 0, 
  0, 0, 2, # VAR(1)模型系数矩阵
  0, 0, 0) # 三个beta的均值
fit <- statespacer(
  y = y.yc,
  self_spec_list = spec,
  collapse = TRUE,
  initial = initial,
  method = "BFGS",
  verbose = TRUE)
## Starting the optimisation procedure at: 2022-01-23 17:10:42
## initial  value -1.899501 
## iter  10 value -6.401173
## iter  20 value -6.930366
## iter  30 value -6.971825
## iter  40 value -6.989413
## iter  50 value -6.991768
## iter  60 value -6.997905
## iter  70 value -7.001514
## iter  80 value -7.003293
## iter  90 value -7.003989
## iter 100 value -7.004314
## final  value -7.004314 
## stopped after 100 iterations
## Finished the optimisation procedure at: 2022-01-23 17:10:51
## Time difference of 8.87524509429932 secs

算法有可能在中间因矩阵不可逆而中断， 这时可以修改初值。 模型估计一次需要大约半分钟。 将计算系统矩阵的函数改写为C++版本可能会加快速度。

平滑得到的\(\beta_{1t}\)的时间序列图形：

plot(date.yc, fit$smoothed$a[, 1], type = 'l', 
  xlab = "year", ylab = "Level of yield curve")

平滑得到的\(\beta_{2t}\)的时间序列图形：

plot(date.yc, fit$smoothed$a[, 2], type = 'l', 
  xlab = "year", ylab = "Slope of yield curve")

平滑得到的\(\beta_{3t}\)的时间序列图形：

plot(date.yc, fit$smoothed$a[, 3], type = 'l', 
  xlab = "year", ylab = "Shape of yield curve")

获取估计的参数：

parameters <- data.frame(
  Parameter = c("lambda", "sigma2", "mu1", "mu2", "mu3"), 
  Value = fit$system_matrices$self_spec,
  SE = fit$standard_errors$self_spec
)
knitr::kable(parameters, digits=4)
Parameter Value SE
lambda 0.0789 0.0017
sigma2 0.0035 0.0002
mu1 7.3857 4.0530
mu2 -1.5138 2.3778
mu3 -0.1008 0.8667

向量自回归系数矩阵：

fit$system_matrices$T$self_spec[1:3, 1:3]
##             [,1]        [,2]        [,3]
## [1,]  0.99097905 -0.01235063 0.008124358
## [2,] -0.02412752  0.94346721 0.061020992
## [3,] -0.02714121 -0.03808673 0.973355560

状态方程误差方差阵：

fit$system_matrices$Q$self_spec[1:3, 1:3]
##             [,1]        [,2]        [,3]
## [1,]  0.06728484 -0.05240711  0.07024729
## [2,] -0.05240711  0.07102692 -0.04152358
## [3,]  0.07024729 -0.04152358  0.24260575

28.2.6 时变系数CAPM模型

考虑时变系数的资产定价模型(CAPM)： \[\begin{aligned} r_t =& \beta_{0t} + \beta_{1t} r_{M,t} + e_t, \ e_t \sim \text{ iid N}(0, \sigma_e^2), \\ \beta_{0,t+1} =& \beta_{0,t} + u_t, \ u_t \sim \text{ iid N}(0, \sigma_{u}^2), \\ \beta_{1,t+1} =& \beta_{1,t} + v_t, \ v_t \sim \text{ iid N}(0, \sigma_{v}^2), \end{aligned}\] 其中\(\{e_t\}\), \(\{u_t\}\), \(\{v_t\}\)相互独立， \(r_t\)是某金融资产的超额收益率， \(r_{M,t}\)是市场的超额收益率。 可以写成状态空间模型 \[\begin{aligned} r_t =& (1, r_{M,t}) \begin{pmatrix} \beta_{0t} \\ \beta_{1t} \end{pmatrix} + e_t, \\ \begin{pmatrix} \beta_{0,t+1} \\ \beta_{1,t+1} \end{pmatrix} =& \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} \beta_{0,t} \\ \beta_{1,t} \end{pmatrix} + \begin{pmatrix} u_t \\ v_t \end{pmatrix} . \end{aligned}\] 令状态变量\(\boldsymbol \alpha_t = (\beta_{0t}, \beta_{1t})^T\), 观测变量为\(r_t\)，且 \[\begin{aligned} Z_t =& (1, r_{M,t}) , \quad H_t = \sigma_e^2, \\ T_t =& I_2, \quad R_t = I, \quad \boldsymbol\eta_t = (u_t, v_t)^T, \\ Q_t =& \begin{pmatrix} \sigma_u^2 & 0 \\ 0 & \sigma_v^2 \end{pmatrix} , \end{aligned}\] 则时变系数CAPM对应于状态空间模型 \[\begin{aligned} r_t =& Z_t \boldsymbol\alpha_t + e_t, \ e_t \sim \text{N}(0, H_t), \\ \boldsymbol\alpha_{t+1} =& T_t \boldsymbol\alpha_t + R_t \boldsymbol\eta_t, \ \boldsymbol\eta_t \sim \text{N}(0, Q_t) . \end{aligned}\]

例28.3 考虑通用动力(GM)股票的月度超额收益率\(r_t\)的CAPM模型， 时间为1990年1月到2003年12月， 以标普500超额收益率为市场超额收益率\(r_{M,t}\)。

数据（收益率单位：百分之一）：

da <- readr::read_table2(
  "m-fac9003.txt",
  col_types = cols(.default=col_double()))[,c("GM", "SP5")]
## Warning: `read_table2()` was deprecated in readr 2.0.0.
## Please use `read_table()` instead.
nda <- nrow(da)
ts.gm <- ts(as.matrix(da),
  start=c(1990, 1), frequency=12)
plot(as.xts(ts.gm),
  main="GM和标普500的月度超额收益率(%)")

library(statespacer)

## 保存自定义设置用的列表
spec <- list()

## 自己规定H矩阵（观测误差的方差阵）
spec$H_spec <- TRUE

## 状态向量维数，2个beta
spec$state_num <- 2

## 待估参数个数3个，包括：
##   $\sigma_e^2$, $\sigma_u^2$, $\sigma_v^2$
## 采用发散先验
spec$param_num <- 3

# 观测方程的Z矩阵，一个$1 \times 2$矩阵，
# $(1, r_{M,t})$，输入为dim=c(1,2,N)数组
spec$Z <- array(c(rbind(1, da[["SP5"]])), dim=c(1,2,nda))

## 状态转移矩阵
spec$Tmat <- diag(1, 2, 2)

## 状态方程中$R$矩阵不起作用
spec$R <- diag(1, 2, 2) # $I_2$

## 初始状态
spec$a1 <- matrix(0, 2, 1)

## 初始状态$\boldsymbol\alpha_1$的方差，用发散先验。
spec$P_inf <- matrix(1E6, 2, 2)

spec$sys_mat_fun <- function(param) {
  ## 输入3个参数，具体意义见前面注释和说明
  
  # \sigma_e^2的值
  sigma2 <- exp(2 * param[1])
  # H是观测方程的误差向量的协方差阵参数。
  H <- sigma2 * diag(1, 1, 1)

  # 状态方程的误差方差阵
  sigmau2 <- exp(2 * param[2])
  sigmav2 <- exp(2 * param[3])
  Q <- diag(c(sigmau2, sigmav2))
  
  # 函数返回模型所需的所有依赖于参数的系统矩阵
  return(list(
    H = H,   # 观测误差随机向量的方差阵
    Q = Q))   # 状态方程的误差随机向量的方差阵
}

因为涉及到依赖于时间的系统矩阵\(Z_t\)， statespacer包没有给出具体做法， 估计程序运行时出错。

fit <- statespacer(
  y = matrix(da[["GM"]], ncol=1),
  self_spec_list = spec,
  initial = c(2,2,2),
  method = "BFGS",
  verbose = TRUE)

B 参考文献

Beijers, Dylan. 2020. statespacer: State Space Modelling in r. https://dylanb95.github.io/statespacer/.
Durbin, J., and S. J. Koopman. 2012. Time Series Analysis by State Space Methods. Oxford University Press.
Tsay, Ruey S. 2010. Analysis of Financial Time Series. 3rd Ed. John Wiley & Sons, Inc.