28 状态空间模型
上一章的局部水平模型是线性高斯状态空间模型的一个简单特例。 本章给出状态空间模型, 举例说明这种模型能够表示的其它模型, 如ARIMA模型,结构时间序列模型, 时变回归模型,有自相关误差的回归模型, 随机波动率模型等, 并给出滤波、平滑、预报公式和参数估计方法。
参考:
28.1 状态空间模型的公式
28.1.1 线性高斯状态空间模型
状态空间模型有许多不同的表达形式, 按照(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 一元结构时间序列模型
下面给出多种时间序列模型的状态空间模型表示和R程序应用。
前一章的局部水平模型是一元结构时间序列模型的特例。 在局部水平模型中增加一个斜率\(\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 \} }\)。
28.2.1 Alcoa现实波动率的局部水平模型
例28.1 考虑Alcoa股票日现实波动率数据的局部水平模型, 时间期间为2003-01-02到2004-05-07, 共340个观测。
读入数据:
<- readr::read_table("aa-3rv.txt",
da col_names=FALSE)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## X1 = col_double(),
## X2 = col_double(),
## X3 = col_double()
## )
<- ts(log(da[[2]])) ts.alcoa
使用局部水平模型,状态空间模型为: \[\begin{aligned} y_t =& \mu_t + e_t, \\ \mu_{t+1} =& \mu_t + \eta_t . \end{aligned}\]
使用statespacer包进行估计, 这个包直接支持对结构时间序列模型进行简化的模型设定:
library(statespacer)
<- statespacer(
ssr1 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: 2023-05-26 09:26:37
## initial value 1.022439
## iter 10 value 0.765458
## iter 20 value 0.764395
## final value 0.764395
## converged
## Finished the optimisation procedure at: 2023-05-26 09:26:37
## Time difference of 0.0184540748596191 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%预测区间:
<- sqrt(ssr1$smoothed$V[1,1,])
smsd plot(ts.alcoa, ylim=c(-1.5, 3))
lines(as.vector(time(ts.alcoa)),
$smoothed$level, col="green")
ssr1lines(as.vector(time(ts.alcoa)),
$smoothed$level - 1.96*smsd,
ssr1lty=2, col="cyan")
lines(as.vector(time(ts.alcoa)),
$smoothed$level + 1.96*smsd,
ssr1lty=2, col="cyan")
legend("bottomleft", lty=c(1,1,2),
col=c("black", "green", "cyan"),
legend=c("Obs", "Smoothed", "95% CL"))
下面改用KFAS包进行估计。 这个包也支持对结构时间序列模型进行简化的模型设定。
library(KFAS)
## 模型设定
<- SSModel(
kfas01a ~ SSMtrend(
ts.alcoa 1, Q=list(matrix(NA))),
H=matrix(NA))
## 超参数估计
<- fitSSM(
kfas01b rep(log(var(ts.alcoa)), 2), method="BFGS")
kfas01a, c("Var_obs"=c(kfas01b$model$H),
"Var_level"=c(kfas01b$model$Q))
## Var_obs Var_level
## 0.230651733 0.005404091
## 用估计的超参数进行滤波平滑
<- KFS(kfas01b$model) kfas01c
超参数估计结果与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两个包的平滑结果是一致的。
dlm包也可以实现这个模型。 在超参数\(\sigma_e^2\)和\(\sigma_{\eta}^2\)已知时, 可以指定模型如下:
library(dlm)
dlm(
FF = 1, # $Z_t$矩阵
V = 0.8, # $\sigma_e^2$
GG = 1, # $T_t$矩阵
W = 0.1, # $\sigma_{\eta}^2$
m0 = 0, # 与$a_1$作用相同,为$\mu_0$期望
C0 = 100 # 与$P_1$作用相同,为$\mu_0$方差
)
## $FF
## [,1]
## [1,] 1
##
## $V
## [,1]
## [1,] 0.8
##
## $GG
## [,1]
## [1,] 1
##
## $W
## [,1]
## [1,] 0.1
##
## $m0
## [1] 0
##
## $C0
## [,1]
## [1,] 100
也可以用dlmModPoly()
指定这样的局部水平模型:
dlmModPoly(
order = 1, # 局部水平模型
dV = 0.8, # $\sigma_e^2$
dW = 0.1, # $\sigma_{\eta}^2$
C0 = 100 # 与$P_1$作用相同,为$\mu_0$方差
)
## $FF
## [,1]
## [1,] 1
##
## $V
## [,1]
## [1,] 0.8
##
## $GG
## [,1]
## [1,] 1
##
## $W
## [,1]
## [1,] 0.1
##
## $m0
## [1] 0
##
## $C0
## [,1]
## [1,] 100
为了对超参数进行最大似然估计, 需要自定义一个给定超参数后输出dlm模型的函数:
<- function(para){
buildFun dlmModPoly(
order = 1,
dV = exp(para[[1]]),
dW = exp(para[[2]])
) }
注意模型中并没有输入数据。
用dlmMLE()
函数输入数据后对超参数进行最大似然估计:
<- dlmMLE(
fit_dlm # 输入的观测值序列
ts.alcoa, parm = c(0, 0), # 超参数对数值的初值
build = buildFun
)$conv fit_dlm
## [1] 0
获得估计的超参数对应的dlm模型:
<- buildFun(fit_dlm$par) dlm_al
估计的超参数值:
c("Var_obs"=V(dlm_al),
"Var_level"=W(dlm_al))
## Var_obs Var_level
## 0.23065283 0.00540346
估计结果与statespacer和KFAS结果基本相同。
为了获得滤波结果,
可以调用dlmFilter(dlm_al)
,
平滑结果使用dlmSmooth(dlm_al)
,
从结果中可以很容易地提取滤波估计与平滑估计,
但是条件分布方差阵估计使用了奇异值分解表示,
提取时要使用dlmSvd2var()
函数转换为方差阵。
○○○○○
28.2.2 尼罗河流量带变点的局部水平模型
例28.2 考虑尼罗河1871-1970年的年流量数据, 单位是\(10^8\)立方米, 在1898-1899年有明显的变点。 用局部水平模型建模并考虑变点的拟合。
基本R的datasets包中Nile变量为长度100的时间序列:
plot(Nile, main="Annual flow of Nile")
abline(v = c(1898, 1899), lty=3, col="yellow")
使用dlm包, 输入超参数建立dlm模型的函数:
library(dlm)
<- function(para){
buildFun <- dlmModPoly(
mod order = 1,
dV = exp(para[[1]]) # $\sigma_e^2$,非时变
)$JW <- matrix(1) # 示性函数,表示$\sigma_{\eta_t}^2$时变
mod# mod$X: 用来保存时变的成分
$X <- cbind(rep(exp(para[[2]]), length(Nile)))
mod# 给t=1898年的观测的$\eta_t$方差放大
<- which(time(Nile) == 1899)[1]
jj $X[jj,1] <- mod$X[jj,1] * (1 + exp(para[[3]]))
mod
mod }
在上面的模型中,
状态方差的扰动项\(\eta_t\)的方差是时变的,
但仅在1899年与其它年份不同。
为了容许模型中的矩阵时变,
在dlm中用JW
作为示性函数,
表示状态方程的扰动项方差阵那些元素是时变的,
然后将时变的内容保存在X
矩阵中,
矩阵的行数等于观测长度\(n\)。
用第三个钞参数表示1899年的状态方程扰动项方差放大倍数(对数尺度)。
三个超参数的最大似然估计:
<- dlmMLE(
fit_nile
Nile,parm = c(0, 0, 0),
build = buildFun)
$conv fit_nile
## [1] 0
获取估计的模型:
<- buildFun(fit_nile$par) dlm_nile
\(\sigma_e^2\)估计:
V(dlm_nile)
## [,1]
## [1,] 16300.33
除去变点以外的所有观测的\(\eta_t\)方差估计, 以及1899年的\(\eta_t\)方差估计:
$X[c(1, which(time(Nile) == 1899)[1]), 1] dlm_nile
## [1] 2.792224e-02 6.048379e+04
可以看出, 非变点处的\(\eta_t\)的标准差, 与观测\(y_t\)的取值范围相比, 基本上等于0。
计算滤波和平滑估计:
<- dlmFilter(Nile, dlm_nile)
filt_nile <- dlmSmooth(filt_nile) sm_nile
提取平滑估计:
<- dropFirst(sm_nile$s) sme_nile
提取平滑分布方差:
<- dlmSvd2var(sm_nile$U.S, sm_nile$D.S) |>
smv_nile unlist() |>
dropFirst()
作原始序列图, 添加平滑结果曲线(红色), 平滑区间(绿色):
plot(Nile, type="o", main="Smoothed Nile Flow")
lines(sme_nile, col="red")
lines(sme_nile - 1.96 * sqrt(smv_nile), lty=3, col="green")
lines(sme_nile + 1.96 * sqrt(smv_nile), lty=3, col="green")
abline(v=c(1898, 1899), lty=3, col="yellow")
这个模型中的\(\eta_t\)的方差是时变的, 水平在1898到1999之间有一个变点, 在变化之前的\(\eta_t\)的方差基本恒定, 水平也基本恒定, 变化之后也是如此。
○○○○○
28.2.3 英国燃气消耗量结构时间序列建模
基本R的datasets包的UKgas变量保存了英国1960年到1986年的季度燃气消耗时间序列。 使用结构时间序列建模, 需要有趋势部分和季节部分, dlm包支持这样的模型设定。
plot(UKgas, main="UK quarterly gas consumption")
因为趋势与季节部分呈现乘性的组合, 改用对数值:
plot(log(UKgas), main="UK quarterly gas consumption(log scale)")
指定模型:
library(dlm)
<- log(UKgas)
lGas <- dlmModPoly() + dlmModSeas(4) dlmGas
函数dlmModPoly()
会指定一个有局部水平和水平的增量(时变斜率)的模型,
dlmModSeas(4)
指定季节模型,
用加号表示这两部分的状态分别建模,
在观测方程中将两部分结果相加。
拟合模型:
<- function(para){
buildFun <- dlmModPoly() + dlmModSeas(4)
dlmGas # 仅指定了时变斜率和季节项的方差,
# 没有指定局部水平的方差
diag(W(dlmGas))[2:3] <- exp(para[1:2])
V(dlmGas) <- exp(para[3])
dlmGas
}<- dlmMLE(
fit_gas
lGas,parm = log(c(10^-5, 10^-3, 10^-3)),
build = buildFun )
$conv fit_gas
## [1] 0
估计得到的模型:
<- buildFun(fit_gas$par) dlmGas
观测方程方差:
drop(V(dlmGas))
## [1] 0.001822492
状态方程方差:
diag(W(dlmGas))[2:3]
## [1] 7.901263e-06 3.308590e-03
获取平滑得到的水平、季节项:
<- dlmSmooth(lGas, mod=dlmGas)
sm_gas <- cbind(
da_tmp
lGas,dropFirst(sm_gas$s[, c(1,3)]))
colnames(da_tmp) <- c("Gas", "Trend", "Seasonal")
plot(da_tmp, type="o",
main = "UK gas consumption structral model")
预报分布的均值和方差可以用dlmForecast()
获得。
下面的预测区间是\(90\%\)置信度的。
<- dlmFilter(lGas, mod = dlmGas)
filt_gas <- dlmForecast(filt_gas, nAhead = 20)
fore_gas # 预测
<- fore_gas$a[,1]
gas_pred20 # 预测标准误差
<- fore_gas$R |>
gas_predse sapply(\(x) sqrt(x[1,1]))
# 将末尾一段与预测结果连接
<- ts.union(
da_tmp window(lGas, start=c(1982,1)),
window(sm_gas$s[,1], start=c(1982,1)),
gas_pred20,- 1.645*gas_predse,
gas_pred20 + 1.645*gas_predse)
gas_pred20 plot(da_tmp, plot.type = "single",
type = "o",
pch = c(1, 0, 20, 3, 3),
col = c("darkgrey", "darkgrey",
"brown", "yellow", "yellow"),
ylab = "Gas consumption prediction"
)
使用dlm发现,
如果在dlmMLE()
中使用parm=rep(0, 3)
,
则将原始数据作自然对数变换的结果比较好,
状态方程扰动方差和观测方程扰动方差很小,
但改为常用对数变换时,
状态方程扰动项方差会增大一两个数量级,
使得上图中的预测区间会大得不可接受。
如果在dlmMLE()
中人为取parm
为较合适的初值,
就可以得到与自然对数变换相同的结果,
这说明dlm计算最大似然估计可能不太可靠,
很可能落入局部极值点。
dlmMLE()
的结果会依赖于参数的初始值parm
的选择。
28.3 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-1} \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.3 考虑太阳黑子数的ARMA建模。 数据为1700-1988的年数据。
数据:
library(datasets)
<- window(sunspot.year, start = 1770, end = 1869)
ts.ss plot(ts.ss)
可以看出,数据有周期性。
用R的arima函数建立AR(3)模型:
<- arima(
ss1 order=c(3,0,0))
ts.ss, 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)中特征多项式的根和复根对应的周期:
<- polyroot(c(1, -coef(ss1)[1:3])); rt1 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)
<- matrix(window(
mat.suns start = 1770, end = 1869))
sunspot.year, <- statespacer(
fit 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: 2023-05-26 09:26:40
## initial value 5.022561
## iter 10 value 4.114363
## final value 4.111163
## converged
## Finished the optimisation procedure at: 2023-05-26 09:26:40
## Time difference of 0.0616450309753418 secs
构造ARMA的结果:
<- rbind(
arma_coeff $system_matrices$AR$ARIMA1,
fit$standard_errors$AR$ARIMA1
fit
)<- cbind(
arma_coeff
arma_coeff,c(fit$smoothed$level[1],
sqrt(fit$system_matrices$Z_padded$level %*%
$smoothed$V[,,1] %*%
fitt(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.4 考虑航空乘客数的建模,1949-1961年的月度数据, 用ARIMA\((0,1,1)(0,1,1)_{12}\)模型。
plot(AirPassengers)
将数据取自然对数,
用arima()
建模:
<- arima(
ap01 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计算:
<- matrix(log(AirPassengers))
mat.ap
# 模型设定列表
<- list(list(
sarima_list s = c(12, 1), ar = c(0, 0), i = c(1, 1), ma = c(1, 1) ))
# 拟合模型
<- statespacer(
fit 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: 2023-05-26 09:26:40
## initial value -1.034434
## final value -1.616321
## converged
## Finished the optimisation procedure at: 2023-05-26 09:26:41
## Time difference of 0.461006164550781 secs
提前估计的参数:
<- rbind(
arma_coeff c(fit$system_matrices$SMA$SARIMA1$S1,
$system_matrices$SMA$SARIMA1$S12),
fitc(fit$standard_errors$SMA$SARIMA1$S1,
$standard_errors$SMA$SARIMA1$S12) )
fitrownames(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
<- rbind(
goodness_fit $system_matrices$Q$SARIMA1,
fit$diagnostics$loglik,
fit$diagnostics$AIC
fit
)rownames(goodness_fit) <- c("Variance", "Loglikelihood", "AIC")
goodness_fit
## [,1]
## Variance 0.001347882
## Loglikelihood 232.750284785
## AIC -3.010420622
结果与arima()
结果相同。
28.4 回归模型
28.4.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.4.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.4.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.4.4 Seatbelts数据建模
R的datasets包的Seatbelts数据是一个多元时间序列, 为月度数据, 保存了英国1969年到1984年与道路交通事故有关的数据。 分量含义:
- DriversKilled: 小汽车交通事故死亡人数。
- drivers: 小汽车交通事故死亡以及重伤人数。
- front: 前排驾乘人员死亡和重伤人数。
- rear: 后排乘客死亡和重伤人数。
- kms: 已行驶里程。
- PetroPrice: 汽油价格。
- VanKilled: 小型厢式载货汽车司机死亡人数。
- law: 强制安全带法律是否已实施的0-1变量,从1983-01-31起实行。
28.4.4.1 数据预处理
因为取正值的变量建模时比较受限制, 所以将这些变量取对数值。
library(statespacer)
library(datasets)
<- as.xts(Seatbelts) |> coredata()
da for(v in c("drivers", "front", "rear", "PetrolPrice", "kms"))
<- log(da[,v]) da[, v]
28.4.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}\]
<- cbind(da[,"drivers"])
y <- statespacer(
ssmr1 y = y, # 因变量
local_level_ind = TRUE, # 是否有局部水平项
# 设置局部水平项误差方差为零,从而非随机:
format_level = matrix(0),
BSM_vec = 12, # 季节项的周期
# 设置季节项误差项方差为零,从而非随机:
format_BSM_list = list(matrix(0)),
# 添加的外生自变量,每个因变量分量需要输入一个自变量矩阵,
# 所以即使是一元因变量也需要用列表
addvar_list = list(as.matrix(
c("PetrolPrice", "law")])),
da[, method = "BFGS",
# 唯一的方差参数$\sigma_e^2$的初值,
# 算法中使用其对数值的2倍以确保取正值
initial = 0.5 * log(var(c(y))),
verbose = TRUE)
## Starting the optimisation procedure at: 2023-05-26 09:26:42
## initial value -0.443442
## final value -0.735372
## converged
## Finished the optimisation procedure at: 2023-05-26 09:26:42
## Time difference of 0.0513920783996582 secs
观测误差方差\(\sigma_e^2\)估计:
$system_matrices$H$H[1,1] ssmr1
## [1] 0.007402481
非随机的局部水平\(\mu\)的平滑估计:
$smoothed$level[1,1] ssmr1
## [1] 6.401571
第一个自变量,即汽油价格对数值,的回归系数:
$smoothed$addvar_coeff[1, 1] ssmr1
## [1] -0.4521301
第二个自变量,即强制安全带法令是否实行,的回归系数, 和相应标准误差估计, t统计量值:
c(coef = ssmr1$smoothed$addvar_coeff[1, 2],
se = ssmr1$smoothed$addvar_coeff_se[1, 2],
tstat = ssmr1$smoothed$addvar_coeff[1, 2] /
$smoothed$addvar_coeff_se[1, 2]) ssmr1
## coef se tstat
## -0.19713947 0.02072792 -9.51081922
由此计算的法令影响的t统计量低于\(-2\), 在0.05水平下显著, 认为法令生效后水平伤亡水平显著降低。
下面将原始数据与平滑拟合结果(包括局部水平、季节项与自变量作用)作图:
plot(da[, c("drivers")], type = "l", ylim = c(6.95, 8.1),
xlab = "year", ylab = "logarithm of drivers")
lines(seq(tsp(Seatbelts)[1], tsp(Seatbelts)[2], 1/tsp(Seatbelts)[3]),
$smoothed$level
ssmr1+ 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.4.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\)相比, 程序只需将零方差改为非零方差:
<- cbind(da[,"drivers"])
y <- statespacer(
ssmr2 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: 2023-05-26 09:26:42
## 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: 2023-05-26 09:26:43
## Time difference of 0.770241975784302 secs
程序中还修改了原来的迭代初值, 原来的迭代初值导致迭代过程出错。 初值需要三个方差参数对应的初值, 这里只给了一个初值, 就重复利用。
估计结果:
## 观测误差方差估计
$system_matrices$H$H[1,1] ssmr2
## [1] 0.003786183
## 局部水平的系统方程的误差方差估计
$system_matrices$Q$level[1,1] ssmr2
## [1] 0.0002676875
## 季节项的系统方程的误差方差估计
$system_matrices$Q$BSM12[1,1] ssmr2
## [1] 1.161438e-06
## 关于汽油价格对数值这一自变量的回归系数估计
$smoothed$addvar_coeff[1, 1] ssmr2
## [1] -0.2913944
第二个自变量,即强制安全带法令是否实行,的回归系数, 和相应标准误差估计, t统计量值:
c(coef = ssmr2$smoothed$addvar_coeff[1, 2],
se = ssmr2$smoothed$addvar_coeff_se[1, 2],
tstat = ssmr2$smoothed$addvar_coeff[1, 2] /
$smoothed$addvar_coeff_se[1, 2]) ssmr2
## coef se tstat
## -0.23773735 0.04631717 -5.13281209
结果表明法令生效后伤亡水平显著降低。
原始数据与拟合结果(拟合结果是平滑结果, 利用了平滑的局部水平与自变量线性组合):
plot(da[, c("drivers")], type = "l", ylim = c(6.95, 8.1),
xlab = "year", ylab = "logarithm of drivers")
lines(seq(tsp(Seatbelts)[1], tsp(Seatbelts)[2], 1/tsp(Seatbelts)[3]),
$smoothed$level
ssmr2+ 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(Seatbelts)[1], tsp(Seatbelts)[2], 1/tsp(Seatbelts)[3]),
$smoothed$BSM12,
ssmr2type = "l", ylim = c(-0.2, 0.3),
xlab = "year", ylab = "stochastic seasonal")
abline(h = 0)
我们将随机的季节项改为非随机:
<- cbind(da[,"drivers"])
y <- statespacer(
ssmr3 y = y, # 因变量
local_level_ind = TRUE, # 是否有局部水平项
# 设置局部水平项误差方差为非零,从而随机:
format_level = matrix(1),
BSM_vec = 12, # 季节项的周期
# 设置季节项误差项方差为零,从而非随机:
format_BSM_list = list(matrix(0)),
# 添加的外生自变量,每个因变量分量需要输入一个自变量矩阵,
# 所以即使是一元因变量也需要用列表
addvar_list = list(as.matrix(
c("PetrolPrice", "law")])),
da[, method = "BFGS",
# 唯一的方差参数$\sigma_e^2$的初值,
# 算法中使用其对数值的2倍以确保取正值
initial = log(var(c(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: 2023-05-26 09:26:44
## initial value -0.577129
## iter 10 value -0.912856
## final value -0.912859
## converged
## Finished the optimisation procedure at: 2023-05-26 09:26:44
## Time difference of 0.247884035110474 secs
估计结果:
## 观测误差方差估计
$system_matrices$H$H[1,1] ssmr3
## [1] 0.004033516
## 局部水平的系统方程的误差方差估计
$system_matrices$Q$level[1,1] ssmr3
## [1] 0.0002681651
## 关于汽油价格对数值这一自变量的回归系数估计
$smoothed$addvar_coeff[1, 1] ssmr3
## [1] -0.2767301
第二个自变量,即强制安全带法令是否实行,的回归系数, 和相应标准误差估计, t统计量值:
c(coef = ssmr3$smoothed$addvar_coeff[1, 2],
se = ssmr3$smoothed$addvar_coeff_se[1, 2],
tstat = ssmr3$smoothed$addvar_coeff[1, 2] /
$smoothed$addvar_coeff_se[1, 2]) ssmr3
## coef se tstat
## -0.2375904 0.0464483 -5.1151578
结果表明法令生效后伤亡水平显著降低。
比较三个模型的AIC:
c(ssmr1$diagnostics$AIC,
$diagnostics$AIC,
ssmr2$diagnostics$AIC) ssmr3
## [1] -1.314494 -1.653950 -1.659052
第三个模型的AIC最小。
28.4.4.4 前后座驾乘人员同时建模的二元模型
考虑前座驾乘人员和后座乘客死亡、重伤人数对数的二元模型, 都以汽油价格对数值、行驶里程对数值和安全带法令为自变量。
<- da[,c("front", "rear")]
y <- da[,c("PetrolPrice", "kms", "law")]
Xmat <- statespacer(y = y, # 输入数据是两列矩阵,二元时间序列
ssmr4 # 观测方程误差项方差阵格式,这里是无限制的二阶方差阵
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: 2023-05-26 09:26:45
## initial value 0.389588
## iter 10 value -1.583525
## iter 20 value -1.676129
## iter 30 value -1.676530
## final value -1.676537
## converged
## Finished the optimisation procedure at: 2023-05-26 09:26:56
## Time difference of 11.4353930950165 secs
观测方程误差方差阵估计:
$system_matrices$H$H ssmr4
## [,1] [,2]
## [1,] 0.005402166 0.004449533
## [2,] 0.004449533 0.008566858
局部水平的系统方程的误差方差阵估计:
$system_matrices$Q$level ssmr4
## [,1] [,2]
## [1,] 0.0002556392 0.0002247045
## [2,] 0.0002247045 0.0002319561
6个回归系数及相应的t统计量:
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,]) |>
::kable(
knitrcol.names=c("", "Coef", "t stat"),
digits=4 )
Coef | t stat | |
---|---|---|
front PetrolPrice | -0.3076 | -2.8986 |
front kms | 0.1518 | 1.1678 |
front law | -0.3370 | -6.8466 |
rear PetrolPrice | -0.0858 | -0.7639 |
rear kms | 0.5504 | 3.7927 |
rear law | 0.0009 | 0.0171 |
与常识一致的是, 后排乘客的伤亡数并不受安全带法令的影响。 汽油价格上升对前排驾乘人员伤亡数有抑制作用。 行驶里程数对前排驾乘人员伤亡数没有显著影响, 但对后排乘客伤亡有正向的影响。
程序中还可以限制某个方差阵的秩,
比如取format_level
为
\[
\begin{pmatrix}
1 & 0 \\
1 & 0
\end{pmatrix}
\]
可以限制水平的系统方差的误差方差阵(\(2 \times 2\)方差阵)为秩等于1。
28.4.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(statespacer)
library(YieldCurve, quietly=TRUE, warn.conflicts=FALSE)
## Warning: package 'YieldCurve' was built under R version 4.2.3
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
<- FedYieldCurve["1984-12-31/2000-12-01"]
xts.yc <- index(xts.yc)
date.yc <- coredata(xts.yc)
y.yc plot(xts.yc, multi.panel=FALSE)
需要自己指定模型的各个矩阵和初值。
## 保存自定义设置用的列表
<- list()
spec ## 自己规定H矩阵(观测误差的方差阵)
$H_spec <- TRUE
spec## 状态向量维数,3个beta, 3个常数beta均值
$state_num <- 6
spec
## 待估参数个数20个,包括:
## 1个 - $\lambda$
## 1个 - $\sigma_e^2$
## 6个 - $\Sigma_{\boldsymbol\eta}\ 3 \times 3$对称阵
## 9个 - $\Phi_{3\times 3}$矩阵
## 3个 - $\boldsymbol\mu$
$param_num <- 20
spec
## 状态方程中$R$矩阵不起作用,方差结构直接编写在矩阵$Q$中
$R <- diag(1, 6, 6) # $I_6$
spec
## 初始状态$\boldsymbol\alpha_1$的方差,元素都取为0,因为平稳
$P_inf <- matrix(0, 6, 6)
spec
## 指定状态向量中不输出到观测的分量,
## 当压缩观测向量为标量时可提高计算效率
$state_only <- 4:6 spec
模型要估计的参数(超参数)为\(\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
元素设置为一个自定义函数:
$sys_mat_fun <- function(param) {
spec## 输入20个参数,具体意义见前面注释和说明
# 8个期限,即$\tau$值
<- c(3, 6, 12, 24, 36, 60, 84, 120)
maturity
# 从取了对数的参数值恢复$\lamabda$的值
<- exp(2 * param[1])
lambda
# \sigma_e^2的值
<- exp(2 * param[2])
sigma2
# H是观测方程的误差向量的协方差阵参数。
# 注意每个$t$对应的观测值是$8\times 1$向量。
<- sigma2 * diag(1, 8, 8)
H
# 观测方程的Z矩阵,一个$8 \times 6$矩阵,
# 前三列对应于$\beta_{1t}$, $\beta_{2t}$, $\beta_{3t}$
# 后三列对应于$\beta$的三个均值,系数为0
<- lambda * maturity
lambda_maturity <- exp(-lambda_maturity)
ze <- matrix(1, 8, 3)
Z 2] <- (1 - ze) / lambda_maturity
Z[, 3] <- Z[, 2] - ze
Z[,
# 状态方程的误差方差阵($6 \times 6$)的左上角$3 \times 3$部分
# 使用6个参数
<- Cholesky(
Q param = param[3:8],
decompositions = FALSE,
format = matrix(1, 3, 3))
# 从输入的9个矩阵元素生成满足平稳性条件的向量自回归系数矩阵$\Phi$
<- CoeffARMA(
Tmat A = array(param[9:17], dim = c(3, 3, 1)),
variance = Q,
ar = 1, ma = 0)$ar[,,1]
# 生成初始状态$\boldsymbol\alpha_1$的方差阵,确保其满足平稳性条件
<- kronecker(Tmat, Tmat)
T_kronecker <- solve(
Tinv diag(1, dim(T_kronecker)[1],
dim(T_kronecker)[2]) - T_kronecker)
<- matrix(Q)
vecQ <- Tinv %*% vecQ
vecPstar <- matrix(vecPstar, dim(Tmat)[1], dim(Tmat)[2])
P_star
## 添加对应于常数均值的部分
# 在观测方程的矩阵Z中增加对应于常数均值的3列
<- cbind(Z, matrix(0, 8, 3))
Z # 在状态方程的误差方差阵中增加对应于常数均值的部分,都是0
<- BlockMatrix(Q, matrix(0, 3, 3))
Q # 在6维的状态方程中,设置转移矩阵$T$
<- cbind(Tmat, diag(1, 3, 3) - Tmat)
Tmat <- rbind(Tmat, cbind(matrix(0, 3, 3), diag(1, 3, 3)))
Tmat # 初始状态的方差阵中对应于常数均值部分方差和协方差为0
<- BlockMatrix(P_star, matrix(0, 3, 3))
P_star
# $\beta$的均值是最后的3个待估参数,
# 也作为$t=1$时$\beta$的初始分布均值
<- matrix(param[18:20], 6, 1)
a1
# 函数返回模型所需的所有系统矩阵
return(list(
H = H, # 观测误差随机向量的方差阵
Z = Z, # 观测方程的观测矩阵
Tmat = Tmat, # 状态方程的状态转移矩阵
Q = Q, # 状态方程的误差随机向量的方差阵
a1 = a1, # 初始状态的均值
P_star = P_star)) # 初始状态的方差阵
}
为了从估计的参数(超参数)提取某些成分,
可以提供一个transform_fun
函数:
$transform_fun <- function(param) {
spec<- exp(2 * param[1])
lambda <- exp(2 * param[2])
sigma2 <- param[18:20]
means return(c(lambda, sigma2, means))
}
下面给20个参数设置算法迭代初值:
<- c(
initial -1, # $0.5 \ln\lambda$
-2, # $0.5 \ln\sigma_e^2$
# $\Sigma_{\boldsymbol\eta}$的LDL分解中D的对角元素的0.5倍对数值
-1, -1, -1,
# $\Sigma_{\boldsymbol\eta}$的LDL分解中L的元素值
0, 0, 0,
# VAR(1)模型系数矩阵
1, 0, 0,
0, 1, 0,
0, 0, 1,
# 三个beta的均值
0, 0, 0)
<- statespacer(
fit y = y.yc,
self_spec_list = spec,
collapse = TRUE,
initial = initial,
method = "BFGS",
control = list(maxit = 200),
verbose = TRUE)
## Starting the optimisation procedure at: 2023-05-26 09:27:02
## initial value 18.089728
## iter 10 value -6.160531
## iter 20 value -6.760946
## iter 30 value -6.953583
## iter 40 value -6.990606
## iter 50 value -6.996315
## iter 60 value -6.999004
## iter 70 value -7.002370
## iter 80 value -7.004278
## iter 90 value -7.005021
## iter 100 value -7.005174
## iter 110 value -7.005240
## iter 120 value -7.005375
## iter 130 value -7.005400
## iter 140 value -7.005430
## iter 150 value -7.005440
## iter 160 value -7.005462
## final value -7.005470
## converged
## Finished the optimisation procedure at: 2023-05-26 09:27:16
## Time difference of 13.7717099189758 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")
获取估计的参数:
<- data.frame(
parameters Parameter = c("lambda", "sigma2", "mu1", "mu2", "mu3"),
Value = fit$system_matrices$self_spec,
SE = fit$standard_errors$self_spec
)::kable(parameters, digits=4) knitr
Parameter | Value | SE |
---|---|---|
lambda | 0.0789 | 0.0017 |
sigma2 | 0.0035 | 0.0002 |
mu1 | 8.2277 | 2.0759 |
mu2 | -2.2748 | 1.4364 |
mu3 | -0.4256 | 0.7191 |
向量自回归系数矩阵:
$system_matrices$T$self_spec[1:3, 1:3] fit
## [,1] [,2] [,3]
## [1,] 0.98488006 -0.01776100 0.005741521
## [2,] -0.01980538 0.94678326 0.062497259
## [3,] -0.03354339 -0.04024349 0.964429517
状态方程误差方差阵:
$system_matrices$Q$self_spec[1:3, 1:3] fit
## [,1] [,2] [,3]
## [1,] 0.06723244 -0.05232934 0.06856474
## [2,] -0.05232934 0.07118337 -0.04047250
## [3,] 0.06856474 -0.04047250 0.24097759
28.4.6 时变系数CAPM模型
28.4.6.1 模型
考虑时变系数的资产定价模型(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.4.6.2 通用动力的例子
例28.5 考虑通用动力(GM)股票的月度超额收益率\(r_t\)的CAPM模型, 时间为1990年1月到2003年12月, 以标普500超额收益率为市场超额收益率\(r_{M,t}\)。
数据(收益率单位:百分之一):
<- readr::read_table(
da_gm "m-fac9003.txt",
col_types = cols(.default=col_double()))[,c("GM", "SP5")]
<- nrow(da_gm)
nda <- ts(as.matrix(da_gm),
ts.gm start=c(1990, 1), frequency=12)
plot(as.xts(ts.gm),
main="GM和标普500的月度超额收益率(%)")
28.4.6.3 滚动计算回归的简单做法
尝试滚动地进行线性回归, 取固定的拟合时间窗口为24个月, 获得动态的系数估计:
<- 24 # 滚动窗口
wh <- nrow(da_gm)
n <- matrix(NA, n, 2)
betamat for(t in wh:n){
<- (t-wh+1):t
tind <-
betamat[t,] lm(GM ~ SP5, data=da_gm[tind,]) |>
coef()
}<- ts(
ts.beta start=c(1990, 1), frequency=12)
betamat, plot(ts.beta)
这个结果可以作为基准。
28.4.6.4 使用KFAS包拟合时变系数回归模型
statespacer包没有提供时变系数回归模型, 而\(Z_t\)依赖于\(t\)时statespacer包没有提供自定义方法。
使用KFAS包。
这个包提供了SSMregression()
函数,
可以指定带有随机游动斜率的回归模型,
但不支持随机截距项,
尝试使用自定义的对应于截距项的自变量会使得滤波平滑失败。
我们使用SScustom()
直接指定各个矩阵。
library(KFAS)
# 将每一行是一个$Z_t$值的矩阵,
# 转换成$1 \times m \times n$的数组,
# $m$是状态空间维数,$n$是样本量
<- function(M){
dim2to3 <- t(M)
M dim(M) <- c(1, dim(M))
M
}
# Zt: 用三维数组表示的$Z_t, t=1,2,\dots,n$,
# 第三下标代表$t$,
# Zt[,,t]是$1 \times m$矩阵$Z_t$
<- dim2to3(cbind(1, da_gm[["SP5"]]))
Zt
<- SSModel(
mod_gm $GM ~ SSMcustom(
da_gmZ = Zt, # 观测方程的载荷阵
T = diag(2), # 状态方程的转移矩阵
# 状态的初始值
a1 = c(0, 1),
# Q是状态方程扰动项方差阵
Q = diag(NA, 2)),
# H是观测方程误差方差阵
H = NA)
<- fitSSM(
fit_gm
mod_gm, inits = c(0,0,0),
method="BFGS") #"Nelder-Mead")
$optim.out$value fit_gm
## [1] 589.0349
估计的随机截距和随机斜率的扰动标准差:
sqrt(diag(fit_gm$model$Q[,,1]))
## [1] 0.0043620203 0.0006551979
估计的回归误差项标准差:
sqrt(fit_gm$model$H[,,1])
## [1] 8.108269
计算平滑结果:
<- KFS(fit_gm$model) res_gm
GM原始值与拟合值(使用平滑)的曲线:
plot(ts.gm[,"GM"], main="GM and fitted")
lines(as.vector(time(ts.gm)),
c(fitted(res_gm, filtered=FALSE)),
col="red")
如果前面的程序中没有指定a1 = c(0, 1)
,
即初始的随机截距、随机斜率,
则拟合上图在左端的拟合效果很差。
估计的时变截距项和时变斜率项:
<- as.matrix(coef(res_gm))
states 2] <- states[,1] + states[,2]
states[,<- states[,2:3]
states plot(ts(states, start=c(1990, 1), frequency=12))
可以看出, 时变的截距项\(\beta_{0t}\)和时变斜率项\(\beta_{1t}\)的变化范围很小。 与滚动回归得到的结果不一致, 基本上没有起到时变的效果。
使用KFAS的SSMregression()
,
仅使用随机斜率项:
<- SSModel(
mod_gm $GM ~ SSMregression(
da_gm~ da_gm$SP5,
# Q是状态方程扰动项方差阵
Q = NA),
# H是观测方程误差方差阵
H = NA)
<- fitSSM(
fit_gm
mod_gm, inits = c(0, 0),
method="BFGS")
<- KFS(fit_gm$model) res_gm
拟合结果:
plot(ts.gm[,"GM"], main="GM and fitted with time-varying beta")
lines(as.vector(time(ts.gm)),
c(fitted(res_gm, filtered=FALSE)),
col="red")
这个拟合结果比较合理。
随机斜率项估计:
<- as.matrix(coef(res_gm))
states plot(ts(states[,2], start=c(1990, 1), frequency=12),
main="Time-varying Regression Slope")
这个结果与滚动窗口回归计算的斜率项还是不一致, 但是变化范围比较小, 有一定合理性。
28.4.6.5 使用dlm包拟合时变系数回归模型
dlm包也支持时变系数回归, 可以用dlm尝试估计时变截距和斜率。
输入超参数的建模函数:
library(dlm)
<- function(para){
buildFun dlmModReg(
cbind(da_gm[["SP5"]]), # 自变量值
dV = exp(para[1]), # 回归误差方差
dW= exp(para[2:3])) # 状态方差误差方差
}
超参数最大似然估计:
<- dlmMLE(
fit_dlm "GM"]],
da_gm[[parm = c(0,-2,-2),
build = buildFun)
$conv fit_dlm
## [1] 0
估计的数值算法收敛。 提取最优参数的模型:
<- buildFun(fit_dlm$par) dlm_gm
回归误差项标准差:
V(dlm_gm) |> sqrt()
## [,1]
## [1,] 8.130086
随机截距项和斜率项的标准差:
W(dlm_gm) |> diag() |> sqrt()
## [1] 9.177050e-04 1.271935e-06
估计得到的状态方程扰动项标准差很小, 比KFAS模型得到的结果小得多, 使得截距和斜率基本等于常数。
利用最大似然估计得到的超参数进行滤波和平滑:
<- dlmFilter(da_gm[["GM"]], dlm_gm)
filt_gm <- dlmSmooth(filt_gm) smth_gm
平滑得到的时变截距项和时变斜率项:
plot(ts(
$s[-1,1:2],
smth_gmstart=c(1990, 1),
frequency=12),
main="Time-varying intercept and slope")
与KFAS的结果不一致。 另外,实际上截距和斜率的变化范围很小, 基本上是常数。