原文链接：http://tecdat.cn/?p=22665
最近咱们被客户要求撰写对于状态空间模型的钻研报告，包含一些图形和统计输入。
状态空间建模是一种高效、灵便的办法，用于对大量的工夫序列和其余数据进行统计推断
摘要
本文介绍了状态空间建模，其观测值来自指数族，即高斯、泊松、二项、负二项和伽马散布。在介绍了高斯和非高斯状态空间模型的根本实践后，提供了一个泊松工夫序列预测的说明性例子。最初，介绍了与拟合非高斯工夫序列建模的其余办法的比拟。
绪论
状态空间模型为几种类型的工夫序列和其余数据的建模提供了一个对立的框架。结构性工夫序列、自回归综合挪动均匀模型（ARIMA）、简略回归、狭义线性混合模型和三次样条平滑模型只是一些能够示意为状态空间模型的统计模型的例子。最简略的一类状态空间模型是线性高斯状态空间模型（也被称为动静线性模型），常常被用于许多迷信畛域。
高斯状态空间模型
本节将介绍无关高斯状态空间模型实践的要害概念。因为卡尔曼滤波（Kalman filtering）背地的算法次要是基于 Durbin 和 Koopman（2012）以及同一作者的相干文章。对于具备间断状态和离散工夫距离的线性高斯状态空间模型 t =1, . . .，n，咱们有
其中 t∼N（0，Ht），ηt∼N（0，Qt）和 α1∼N（a1，P1）互相独立。咱们假如 yt 是一个 p×1，αt+ 1 是一个 m×1，ηt 是一个 k×1 的向量。α = (α > 1 , . . . , α> n) >，同样 y = (y > 1 , . . , y> n) >。
状态空间建模的次要指标是在给定观测值 y 的状况下取得潜状态 α 的常识。这能够通过两个递归算法实现，即卡尔曼滤波和平滑算法。从卡尔曼滤波算法中，咱们能够失去先行一步的预测后果和预测误差
和相干的协方差矩阵
利用卡尔曼滤波的后果，咱们建设了状态平滑方程，在工夫上向后运行，产生了
对于烦扰项 t 和 ηt，对于信号 θt = Ztαt，也能够计算相似的平滑预计。
高斯状态空间模型的例子
当初通过例子来阐明。咱们的工夫序列包含 1969-2007 年 40-49 岁年龄组每年每 10 万人中酒精相干的死亡人数（图 1）。数据取自统计局。对于观测值 y1, … . , yn，咱们假如在所有 t = 1, . . . , n，其中 μt 是一个随机游走的漂移过程
ηt∼N(0, σ2 η)。假如咱们没有对于初始状态 μ1 或斜率 ν 的先验信息。这个模型能够用状态空间的模式来写，定义为
在 KFAS 中，这个模型能够用以下代码来写。为了阐明问题，咱们手动定义所有的零碎矩阵，而不采纳默认值。
R> Zt <- matrix(c(1, 0), 1, 2)
R> model_gaussian <-Model(deaths / population ~ -1 +custom(Z = Zt)
第一个参数是定义观测值的公式（左侧 \~）和状态方程的构造（右侧）。这里死亡人数 / 人口是一个单变量工夫序列，状态方程是用矩阵来定义的，为了放弃模型的可识别性，截距项用 - 1 省略。观测程度方差通过参数 H 定义，NA 值代表未知方差参数 σ 2 和 σ 2 η。预计之后，进行过滤和平滑递归。
KF(fit_gaussian)
在这种状况下，最大似然估计值对于 σ 2 是 9.5，对于 σ 2 η 是 4.3。从卡尔曼滤波算法中，咱们失去了对状态的一步超前预测，at = (µt , νt)。请留神，即便斜率项 ν 在咱们的模型中被定义为工夫不变量（νt = ν），它也是由卡尔曼滤波算法递归预计的。因而，在每个工夫点 t，当新的观测值 yt 可用时，ν 的估计值被更新，以思考到 yt 所提供的新信息。在卡尔曼滤波完结时，an+ 1 给出了咱们对所有数据下恒定斜率项的最终预计。这里斜率项被预计为 0.84，标准误差为 0.34。对于 µt，卡尔曼滤波给出了一步的预测，然而因为状态是时变的，如果咱们对 t =1, …, n 的 µt 估计值感兴趣，咱们还须要运行平滑算法。n 的估计值。
图 1 显示了带有一步预测（红色）和平滑化（蓝色）的随机行走过程 µt 的估计值的察看后果。请留神典型的模型；在工夫 t，卡尔曼滤波器计算一步向前预测误差 vt = yt – µt，并应用它和先前的预测来修改下一个工夫点的预测。在这里，这在系列的开始阶段最容易看到，咱们的预测仿佛落后于观测值一个工夫步长。另一方面，平滑算法同时思考了每个工夫点的过来和将来的数值，从而产生了更平滑的潜过程的预计。
点击题目查阅往期内容
卡尔曼滤波器：用 R 语言中的 KFAS 建模工夫序列
左右滑动查看更多
01
02
03
04
非高斯状态空间模型的例子
与酒精无关的死亡也能够天然地被建模为泊松过程。当初咱们的观测值 yt 是第 t 年与酒精无关的死亡的理论计数，而变动的人口规模则由裸露项 ut 来思考。状态方程放弃不变，但察看方程当初的模式是 p(yt |µt) = Poisson(ute µt)。
R> Model(deaths ~ -1 +
+ distribution = "poisson")
与之前的高斯模型相比，咱们当初须要用参数 distribution（默认为 “ 高斯 ”）来定义观测数据的散布。咱们还通过参数 u 来定义裸露项，并应用 a1 和 P1 的默认值。在这个模型中，只有一个未知参数，即 σ 2 η。这个参数被预计为 0.0053，然而高斯模型和泊松模型之间 σ 2 η 的理论值不能间接比拟，因为不同模型对 µt 的解释不同。泊松模型的斜率项预计为 0.022，标准误差为 1.4×10-4，对应于死亡人数每年减少 2.3%。
图 2 显示了以高斯过程（蓝色）和泊松过程（红色）为模型（每 10 万人的死亡人数）的平滑预计。
任意的状态空间模型
通过联合后面的办法，能够绝对容易地构建大量的模型。对于这样做还不够的状况，能够通过间接定义零碎矩阵来构建任意状态空间模型。作为一个例子，咱们批改了之前的泊松模型，减少了一个额定的白噪声项，试图捕获数据的可能的适度离散。当初咱们的泊松强度模型是 ut exp(µt + t)，即
其中 ηt∼N(0, σ2 η)如前，t∼N(0, σ2)。这个模型能够用状态空间的模式来写，定义为
Model(deaths ~ trend(2, Q = list(NA, 0)) +
distribution = "poisson")
因为模型蕴含 P1 中的未知参数，咱们须要提供一个特定的模型更新函数。
R> update <- function(pars, model) {+ model[ "custom"] <- exp(pars)
+ }
<!—->
fit(model_poisson,method = "BFGS")
从图 3 中咱们看到，高斯构造工夫序列模型和带有额定白噪声的泊松构造工夫序列模型对平滑趋势 µt 的预计简直没有区别。这是因为泊松过程的强度绝对较高。
例子
我当初用一个比后面的例子更残缺的例子来阐明 KFAS 的应用。数据还是由酒精无关的死亡组成，但当初有四个年龄组，即 30-39 岁、40-49 岁、50-59 岁和 60-69 岁，被一起作为一个多变量泊松模型来建模。
1969-2012 年的死亡人数和相应年龄组的年人口规模都有，但作为阐明，咱们只应用 2007 年之前的数据，并对 2008-2013 年进行预测。图 4 显示了所有年龄组的每 10 万人的死亡人数。
ts.plot(window(data[, 1:4] / data[, 5:8], end = 2007)
这里我抉择了之前应用的泊松模型的一个多变量扩大。
这里 μt 是带有漂移成分的随机游走，νt 是一个恒定的斜率，t 是一个额定的白噪声成分，用于捕获序列的额定变动。我对程度和噪声成分的协方差构造不做限度。模型（4）能够用 KFAS 构建如下。
R> Model(Pred[, 1:4] ~
+ trend(2, Q = list(matrix(NA, 4, 4))
distribution = "poisson"
更新函数为
R> updatefn <- function(pars, model, ...){+ model[ etas] <- crossprod(Q)
+ crossprod(Q)
+ model
+ }
咱们能够先不通过模仿来预计模型参数，而后用这些估计值作为初始值再次运行重要性抽样的预计程序。在这种状况下，从重要性抽样步骤失去的后果实际上与从初始步骤失去的后果雷同。
fit(model, update,
+ method = "BFGS")
R> fit <- fit(model, updatefn = updatefn, inits =optimpar)
应用拟合模型的提取办法，咱们能够查看预计的协方差和相关矩阵。
R> varcordel["Q", "level"]
R> varcordel["Q", "custom"]
状态空间模型的参数估计通常工作量很大，因为似然面蕴含多个最大值，从而使优化问题高度依赖于初始值。通常状况下，未知参数与未察看到的潜在状态无关，如本例中的协方差矩阵，简直没有先验常识。
因而，要猜出好的初始值是很有挑战性的，特地是在更简单的环境中。因而，在能够正当地确定找到适当的最优值之前，倡议应用多种初始值配置，可能有几种不同类型的优化办法。这里咱们应用察看到的系列的协方差矩阵作为协方差构造的初始值。在非高斯模型的状况下，另一个问题是，似然计算是基于迭代程序的，该程序应用一些终止条件（如对数似然的绝对变动）进行，因而对数似然函数实际上蕴含一些噪声。这反过来又会影响 BFGS 等办法的梯度计算，在实践上能够失去不牢靠的后果。因而，有时倡议应用无导数的办法，如 Nelder-Mead。另一方面，BFGS 通常比 Nelder-Mead 快得多，因而我更违心先尝试 BFGS，至多在初步剖析中。咱们能够计算出状态的平滑预计。
R> out <- KF(model,)
咱们看到残差之间偶然有滞后的穿插相干，但总体上咱们能够对咱们的模型绝对称心。当初咱们能够用咱们预计的模型预测 2008-2013 年每个年龄组与酒精无关的死亡强度 e θt。因为咱们的模型是工夫变动的（u 变动），咱们须要通过 newdata 参数为将来的察看样本提供模型。
predict(model,
+ newdata +
+ interval = "confidence",)
for (i in 1:4) ts.plot(data[, i]/data[, 4 + i], trend[, i], pred[[i]]
图 7 显示了察看到的死亡人数，1969-2007 年的平滑趋势，以及 2008-2013 年的预测，还有 95% 的预测区间。当咱们将咱们的预测与实在的察看后果进行比拟时，咱们看到在事实中，最年长的年龄组（60-69 岁）的死亡人数略有减少，而在预测期间，另一个年龄组的死亡人数大幅降落。局部起因是在此期间酒精生产总量简直枯燥降落，而这又可能是因为 2008 年、2009 年和 2012 年酒精税的减少造成的。
探讨
状态空间模型提供了解决一大类统计问题的工具。在这里，我介绍了一个用于线性状态空间建模的办法。
点击文末 “浏览原文”
获取全文残缺代码数据资料。
本文选自《R 语言状态空间模型和卡尔曼滤波预测酒精死亡人数工夫序列》。
点击题目查阅往期内容
matlab 实现扩大卡尔曼滤波 (EKF) 进行故障检测 \
卡尔曼滤波器：用 R 语言中的 KFAS 建模工夫序列 \
状态空间模型：卡尔曼滤波器 KFAS 建模工夫序列 \
R 语言用 LOESS(部分加权回归)节令趋势合成（STL）进行工夫序列异样检测 \
应用 R 语言随机稳定模型 SV 解决工夫序列中的随机稳定率 \
PYTHON 用时变马尔可夫区制转换（MRS）自回归模型剖析经济工夫序列 \
R 语言无限混合模型 (FMM,finite mixture model)EM 算法聚类分析间歇泉喷发工夫 \
长短期记忆网络 LSTM 在工夫序列预测和文本分类中的利用 \
Python 随机稳定率 (SV) 模型对标普 500 指数工夫序列波动性预测 \
R 语言中 ARMA，ARIMA（Box-Jenkins），SARIMA 和 ARIMAX 模型用于预测工夫序列数据 \
R 语言应用 ARIMAX 预测失业率经济工夫序列数据 \
R 语言用 ARIMA 模型，ARIMAX 模型预测冰淇淋生产工夫序列数据 \
R 语言经济学：动静模型均匀 (DMA)、动静模型抉择(DMS) 预测原油工夫序列价格