NO.39.tip： 结构时间序列模型

背景

我们将讨论时间序列预测（time series forecasting），这是在给定到目前为止数据的情况下，计算未来观测的预测分布问题，即计算 $p({\mathbf{y}}_{t+h}\mid {\mathbf{y}}_{1:t})$。该模型还可以选择性地以已知的未来输入为条件，以获得 $p({\mathbf{y}}_{t+h}\mid {\mathbf{y}}_{1:t},{\mathbf{u}}_{1:t+h})$。目前有很多方法可以求解这个问题我们将重点关注结构时间序列（Structural Time Series, STS）模型，该模型是根据线性高斯状态空间模型定义的。

许多经典的时间序列方法，例如 ARMA（AutoRegressive Moving Average，自回归移动平均）方法，可以表示为结构时间序列模型。然而，结构时间序列方法具有更大的灵活性。例如，我们可以创建非线性、非高斯，甚至是层次扩展模型，如下文所述。

概述

结构时间序列模型的基本思想是将观测到的标量时间序列表示为 $C$ 个单独分量（component）的总和：

$$f(t)=f_1(t)+f_2(t)+\cdots +f_C(t)+{ϵ}_t\text{\hspace{1em}(1)}$$

其中，${ϵ}_t\sim \mathcal{N}(0,{\sigma }^2)$。例如，我们可能有一个季节性分量，它导致观测值上下波动，还有一个增长分量，它导致观察值随着时间的推移而变大。每个潜在过程 $f_c(t)$ 都由一个线性高斯状态空间模型建模，该模型也称为动态线性模型。由于这些模型均是线性的，我们可以将这些模型组合成一个单独的线性高斯状态空间模型。特别是，在标量观测的情况下，该模型具有以下形式：

$$p({\mathbf{z}}_t\mid {\mathbf{z}}_{t-1},\mathbf{\theta })=\mathcal{N}({\mathbf{z}}_t\mid \mathbf{F}{\mathbf{z}}_{t-1},\mathbf{Q})\text{\hspace{1em}(2)}$$

$$p(y_t\mid {\mathbf{z}}_t,\mathbf{\theta })=\mathcal{N}(y_t\mid \mathbf{H}{\mathbf{z}}_t+{\mathbf{\beta }}^{\top }{\mathbf{u}}_t,{\sigma }_y^2)\text{\hspace{1em}(3)}$$

其中，$\mathbf{F}$ 和 $\mathbf{Q}$ 是分块结构矩阵，每个分量对应一个分块。然后，向量 $\mathbf{H}$ 将每个分量的所着相关片段相加，以生成总体均值。注意，矩阵 $\mathbf{F}$ 和 $\mathbf{H}$ 是固定稀疏矩阵，两个矩阵均可以从模型相应分量的形式中推导出来，如下所述。因此，唯一的模型参数是方差项 $\mathbf{Q}$ 和 ${\sigma }_y^2$，以及可选的回归系数 $\mathbf{\beta }$。我们可以将其推广到向量值观测的情况，如下所示：

$$p({\mathbf{z}}_t\mid {\mathbf{z}}_{t-1},\mathbf{\theta })=\mathcal{N}({\mathbf{z}}_t\mid \mathbf{F}{\mathbf{z}}_{t-1},\mathbf{Q})\text{\hspace{1em}(4)}$$

$$p({\mathbf{y}}_t\mid {\mathbf{z}}_t,\mathbf{\theta })=\mathcal{N}({\mathbf{y}}_t\mid \mathbf{H}{\mathbf{z}}_t+\mathbf{D}{\mathbf{u}}_t,\mathbf{R})\text{\hspace{1em}(5)}$$

源与流

结构化构建块

我们将讨论常见的结构时间序列模型的构建块。

局部水平模型

最简单的潜在动力过程被称为局部水平模型（local level model）。局部水平模型假设观测值 $y_t\in \mathbb{R}$ 是由具有（潜在）均值 ${\mu }_t$ 的高斯分布产生的，该高斯分布根据随机游走随时间进行如下的演变：

$$y_t={\mu }_t+{ϵ}_{y,t},{ϵ}_{y,t}\sim \mathcal{N}(0,{\sigma }_y^2)\text{\hspace{1em}(6)}$$

$${\mu }_t={\mu }_{t-1}+{ϵ}_{\mu ,t},{ϵ}_{\mu ,t}\sim \mathcal{N}(0,{\sigma }_{\mu }^2)\text{\hspace{1em}(7)}$$

我们还假设 ${\mu }_1\sim \mathcal{N}(0,{\sigma }_{\mu }^2)$。因此，任何未来步骤的潜在均值都服从分布 ${\mu }_t\sim \mathcal{N}(0,t{\sigma }_{\mu }^2)$，这意味着方差随时间增长。我们也可以使用自回归（AR）过程，即 ${\mu }_t=\rho {\mu }_{t-1}+{ϵ}_{\mu ,t}$，其中 $|\rho |<1$，从而具有平稳分布 ${\mu }_{\infty }\sim \mathcal{N}\left(0,\frac{{\sigma }_{\mu }^2}{1-{\rho }^2}\right)$，因此不确定性增长到有限的渐近值，而不是无限制地增长。

局部线性模型

许多时间序列呈现出向上或向下的线性趋势（至少在局部是这样的趋势）。我们可以通过让水平值 ${\mu }_t$ 在每一步中变化 ${\delta }_{t-1}$（表示直线在 $\Delta t=1$ 间隔内的斜率）来对此进行建模：

$${\mu }_t={\mu }_{t-1}+{\delta }_{t-1}+{ϵ}_{\mu ,t}\text{\hspace{1em}(8)}$$

斜率本身也遵循随机游走：

$${\delta }_t={\delta }_{t-1}+{ϵ}_{\delta ,t}\text{\hspace{1em}(9)}$$

并且 ${ϵ}_{\delta ,t}\sim \mathcal{N}(0,{\sigma }_{\delta }^2)$。这被称为局部线性趋势（local linear trend）模型。

我们可以通过定义以下动力学模型将这两个过程结合起来：

$$\underset{{\mathbf{z}}_t}{\underbrace{\left(\begin{array}{c}{\mu }_t\\ {\delta }_t\end{array}\right)}}=\underset{\mathbf{F}}{\underbrace{\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)}}\underset{{\mathbf{z}}_{t-1}}{\underbrace{\left(\begin{array}{c}{\mu }_{t-1}\\ {\delta }_{t-1}\end{array}\right)}}+\underset{{\mathbf{ϵ}}_t}{\underbrace{\left(\begin{array}{c}{ϵ}_{\mu ,t}\\ {ϵ}_{\delta ,t}\end{array}\right)}}\text{\hspace{1em}(10)}$$

对于发射模型，我们有：

$$y_t=\underset{\mathbf{H}}{\underbrace{\left(\begin{array}{cc}1& 0\end{array}\right)}}\underset{{\mathbf{z}}_t}{\underbrace{\left(\begin{array}{c}{\mu }_t\\ {\delta }_t\end{array}\right)}}+{ϵ}_{y,t}\text{\hspace{1em}(11)}$$

我们也可以使用斜率的自回归模型。即

$${\delta }_t=D+\rho ({\delta }_{t-1}-D)+{ϵ}_{\delta ,t}\text{\hspace{1em}(12)}$$

其中，$D$ 是 $\delta$ 将恢复到的长期斜率。这被称为半局部线性趋势（semilocal linear trend）模型，这种模型比较适用于长期预测。

在统计学界，将动力学记作 ${\mathbf{\alpha }}_t={\mathbf{T}}_t{\mathbf{\alpha }}_{t-1}+{\mathbf{c}}_t{\mathbf{R}}_t{\mathbf{\eta }}_t$，将观测结果记作 ${\mathbf{y}}_t={\mathbf{Z}}_t{\mathbf{\alpha }}_t+{\mathbf{\beta }}^{\top }{\mathbf{x}}_t+{\mathbf{H}}_t{\mathbf{ϵ}}_t$，其中 ${\mathbf{\eta }}_t\sim \mathcal{N}(0,\mathbf{I})$ 和 ${\mathbf{ϵ}}_t\sim \mathcal{N}(0,\mathbf{I})$。

添加协变量

我们可以很容易地将协变量 ${\mathbf{u}}_t$ 包含到模型中，以提高预测准确率。如果使用线性模型，则我们有

$$y_t={\mu }_t+{\mathbf{\beta }}^{\top }{\mathbf{u}}_t+{ϵ}_{y,t}\text{\hspace{1em}(13)}$$

注意，在预测未来时，我们需要某种方法来预测未来的输入值 ${\mathbf{u}}_{t+h}$。一个简单的方法就是假设未来的输入与现在相同，即 ${\hat{\mathbf{u}}}_{t+h}={\mathbf{u}}_t$。

对季节性进行建模

许多时间序列也表现出季节性（seasonality），即周期性波动。这可以通过创建一个由一系列周期项 $s_t$ 组成的潜在过程来建模。为了对周期性进行建模，我们确保在 $S$ 个步骤的完整周期内，这些项的总和（平均值）为零：

$$s_t=-\sum _{k=1}^{S-1}{s}_{t-k}+{ϵ}_{s,t},{ϵ}_{s,t}\sim \mathcal{N}(0,{\sigma }_s^2)\text{\hspace{1em}(14)}$$

例如，对于 $S=4$，我们有 $s_t=-(s_{t-1}+s_{t-2}+s_{t-3})+{ϵ}_{s,t}$。我们可以通过将最后的 $S-1$ 个季节叠加到状态向量中，将其转换为一阶模型。

将各种潜在过程结合起来

我们可以将各种潜在过程（局部水平、线性趋势和季节周期）组合成单个线性高斯状态空间模型，因为稀疏图结构可以由稀疏矩阵编码。更准确地说，转移模型变成：

$$\underset{{\mathbf{z}}_t}{\underbrace{\left(\begin{array}{c}{s}_t\\ s_{t-1}\\ s_{t-2}\\ {\mu }_t\\ {\delta }_t\end{array}\right)}}=\underset{\mathbf{F}}{\underbrace{\left(\begin{array}{ccccc}-1& -1& -1& 0& 0\\ 1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 1\\ 0& 0& 0& 0& 1\end{array}\right)}}\underset{{\mathbf{z}}_{t-1}}{\underbrace{\left(\begin{array}{c}{s}_{t-1}\\ s_{t-2}\\ s_{t-3}\\ {\mu }_{t-1}\\ {\delta }_{t-1}\end{array}\right)}}+\mathcal{N}\left(0,\mathrm{diag}([{\sigma }_s^2,0,0,{\sigma }_{\mu }^2,{\sigma }_{\delta }^2])\right)\text{\hspace{1em}(15)}$$

定义了模型之后，我们可以使用卡尔曼滤波器来计算 $p({\mathbf{z}}_t\mid {\mathbf{y}}_{1:t})$，然后通过在潜在空间中向前滚动来及时进行向前预测，然后预测输出：

$$p({\mathbf{y}}_{t+1}\mid {\mathbf{y}}_{1:t})=\int p({\mathbf{y}}_{t+1}\mid {\mathbf{z}}_{t+1})p({\mathbf{z}}_{t+1}\mid {\mathbf{z}}_t)p({\mathbf{z}}_t\mid {\mathbf{y}}_{1:t})d{\mathbf{z}}_t\text{\hspace{1em}(16)}$$

这可以以闭合形式进行计算。

模型拟合

一旦指定了模型的形式，我们就需要学习模型参数 $\mathbf{\theta }=(\mathbf{D},\mathbf{R},\mathbf{Q})$，因为 $\mathbf{F}$ 和 $\mathbf{H}$ 固定在结构块指定的值上，并且 $\mathbf{B}=0$。常见的方法基于最大似然估计和贝叶斯推理。后一种方法被称为贝叶斯结构化时间序列（Bayesian Structural Time Series, BSTS）建模，通常使用以下共轭先验：

$$p(\mathbf{\theta })=\mathrm{MNIW}(\mathbf{R},\mathbf{D})\mathrm{IW}(\mathbf{Q})\text{\hspace{1em}(17)}$$

或者，如果存在大量协变量，我们可以对回归系数 $\mathbf{D}$ 使用稀疏性促进先验分布（例如，“尖峰”和“平板”分布）。

预测

一旦在历史数据集上估计了参数，我们就可以对新的时间序列进行推理，使用卡尔曼滤波器计算 $p({\mathbf{z}}_t\mid {\mathbf{y}}_{1:t},{\mathbf{u}}_{1:t},\mathbf{\theta })$。在给定当前后验的情况下，我们可以及时"向前滚动"，通过计算 $p({\mathbf{y}}_{t+h}\mid {\mathbf{y}}_{1:t},{\mathbf{u}}_{1:t+h},\mathbf{\theta })$ 来预测未来 $h$ 步之后的观测值。如果参数不确定，我们可以从后验 $p(\mathbf{\theta }\mid {\mathbf{y}}_{1:t},{\mathbf{u}}_{1:t})$ 中进行采样，然后对预测进行蒙特卡罗平均。

时间序列干预的因果影响

我们讨论了在给定一些准实验(quasi-experimental)时间序列数据的情况下，如何估计干预的因果效应。（术语"准实验性"是指在干预下收集数据，但不使用随机分配。）例如，假设 $y_t$ 是某家公司在时间 $t$ 的网页点击率（Click Through Rate, CTR）。该公司在时间 $n$ 发起广告活动，并观察干预前的结果 ${\mathbf{y}}_{1:n}$ 和干预后的结果 ${\mathbf{y}}_{n+1:m}$。（这是一个中断时间序列(interrupted time series)的例子，因为"自然"过程在某个时刻受到了干扰。）

一个自然的问题是：如果公司没有开展广告活动，网页点击率会是什么？这是一个反事实(counterfactual)的问题。如果我们可以预测这个反事实时间序列 ${\mathbf{y}}_{n+1:m}$，那么我们就可以将实际值 $y_t$ 与预测值 $y_t$ 进行比较，并使用比较结果来估计干预的因果影响(causal impact)。

为了预测反事实结果，我们使用结构时间序列(STS)模型。结构时间序列模型是线性高斯状态空间模型，其中箭头根据时间箭头(arrow of time)具有自然的因果解释；因此结构时间序列是一种结构方程模型，同时也是一种结构因果模型。结构因果模型的使用使得我们能够在给定观测数据的情况下推理噪声变量的潜在状态；然后，我们可以将时间"回滚"到干预点，在那里，我们使用孪生网络方法进行反事实推理，通过在新版本的模型中"在时间上向前滚动"，来探索一条与实际路径不同的"岔路"。

计算反事实预测

为了更详细地解释该方法，请考虑孪生网络。干预发生在时间 $n=2$ 之后，总共有 $m=4$ 个观察结果。我们观察干预前的两个数据点 ${\mathbf{y}}_{1:2}$ 和干预后的两个点 ${\mathbf{y}}_{3:4}$。我们假设观测是由随时间演变的潜在状态 ${\mathbf{z}}_{1:4}$ 产生的。这些状态受外生噪声项的影响，这些噪声项可以表示任何一组未建模的因素，例如经济状态。此外，我们还有外生协变量 ${\mathbf{x}}_{1:m}$。

考虑用于估计在时间步长 $n=2$ 之后发生的干预的因果影响的孪生网络状态空间模型。我们有 $m=4$ 个实际观测值，表示为 $y_{1:4}$。由于我们假设 ${\mathbf{z}}_{3:T}$ 来自不同的分布（即干预后分布），因此我们将传入弧切割为 ${\mathbf{z}}_3$。然而，在图底部所示的反事实世界中（带有波浪符号），我们假设分布与过去相同，因此信息沿着链不间断地流动

为了预测如果我们没有进行干预（一个使用 $a=0$ 表示的事件）会发生什么，我们复制干预后发生的模型部分，并使用模型进行预测。目标是计算反事实分布 $p({\mathbf{y}}_{n+1:m}\mid {\mathbf{y}}_{1:n},{\mathbf{x}}_{1:m})$，其中 $y_t$ 表示当动作为 $a=0$ 时的反事实结果。我们可以按下式计算这种反事实分布：

$$p({\mathbf{y}}_{n+1:m}\mid {\mathbf{y}}_{1:n},{\mathbf{x}}_{1:m})=\int p({\mathbf{y}}_{n+1:m}\mid {\mathbf{z}}_{n+1:m},{\mathbf{x}}_{n+1:m},\mathbf{\theta })p({\mathbf{z}}_{n+1:m}\mid {\mathbf{z}}_n,\mathbf{\theta })\times \text{\hspace{1em}(18)}$$

$$p({\mathbf{z}}_n,\mathbf{\theta }\mid {\mathbf{x}}_{1:n},{\mathbf{y}}_{1:n})d\mathbf{\theta }d{\mathbf{z}}_{n}d{\tilde{\mathbf{z}}}_{n+1:m}\text{\hspace{1em}(19)}$$

其中，

$$p({\mathbf{z}}_n,\mathbf{\theta }\mid {\mathbf{x}}_{1:n},{\mathbf{y}}_{1:n})=p({\mathbf{z}}_n\mid {\mathbf{x}}_{1:n},{\mathbf{y}}_{1:n},\mathbf{\theta })p(\mathbf{\theta }\mid {\mathbf{x}}_{1:n},{\mathbf{y}}_{1:n})\text{\hspace{1em}(20)}$$

对于线性高斯状态空间模型，$p({\mathbf{z}}_n\mid {\mathbf{x}}_{1:n},{\mathbf{y}}_{1:n},\mathbf{\theta })$ 可以使用卡尔曼滤波计算，并且 $p(\mathbf{\theta }\mid {\mathbf{y}}_{1:n},{\mathbf{x}}_{1:n})$ 可以使用马尔可夫链蒙特卡罗方法或变分推理计算。

我们可以使用来自上述后验预测分布的样本来计算每个时间步长的处理效果(treatment effect)分布的蒙特卡罗近似值 ${\tau }_t^i=y_t-{\tilde{y}}_t^i$，其中索引 $i$ 指向后验样本。我们还可以使用 ${\sigma }_t^i={\sum }_{s=n+1}^t{\tau }_s^i$ 来近似累积因果影响的分布。（这些量中存在不确定性，既有关于控制模型真实参数的认识不确定性，也有由于系统和观测噪声引起的随机不确定性。

方法的有效性所需的假设

方法的有效性基于以下三个假设：

• 可预测性：我们假设在现有数据的情况下，我们的模型可以充分预测结果。我们可以通过使用回溯(backcasting)来检查可预测性，在回溯中，我们对部分历史数据进行预测。
• 无影响性：我们假设干预不会改变未来的协变量 ${\mathbf{x}}_{n+1:m}$。我们可以通过将每个协变量作为结果变量来运行该方法，从而验证无影响性。
• 稳定性：我们假设，如果没有进行干预，处理前阶段的结果模型将在处理后阶段继续适用。我们可以通过观察如果处理在时间上提前，是否可以预测效果来检查稳定性。

与综合控制的比较

使用其他"捐赠"（donor）时间序列 ${\mathbf{\beta }}^{\top }{\mathbf{x}}_t$ 的线性组合本质上类似于合成控制（synthetic control）的概念。然而，我们并不局限于捐赠时间序列的凸组合。事实上，当我们有许多协变量时，我们可以使用"尖峰"和"平板"先验或马蹄先验来选择相关的协变量。此外，即使我们只观察结果序列 $y_t$，而没有任何其他平行时间序列，结构时间序列方法也可以应用。

Prophet

Prophet（先知）是 Facebook 上一个流行的时间序列预测库。该时间序列预测库拟合以下形式的广义加性模型：

$$y(t)=g(t)+s(t)+h(t)+{\mathbf{w}}^{\top }\mathbf{x}(t)+{ϵ}_t\text{\hspace{1em}(21)}$$

其中，$g(t)$ 是趋势函数，$s(t)$ 为季节性波动（使用应用于正弦基集的线性回归建模），$h(t)$ 为由稀疏"假日效应"组成的可选集合，$\mathbf{x}(t)$ 是由（可能滞后的）协变量组成的可选集，$\mathbf{w}$ 为回归系数，$ϵ(t)$ 则为残差噪声项（假设为独立同分布的高斯分布）。

Prophet 是一个回归模型，而不是自回归模型，因为该模型在给定时间戳 $t$ 和协变量 ${\mathbf{x}}_{1:T}$ 的情况下，预测时间序列 ${\mathbf{y}}_{1:T}$，但不以过去对 $y$ 的观测为条件。为了对时间的依赖性进行建模，假设趋势函数是具有 $S$ 个变化点的分段线性趋势，这些变化点在时间上均匀分布。也就是说，该模型具有以下形式：

$$g(t)=(k+\mathbf{a}(t{)}^{\top }\mathbf{\delta })t+(m+\mathbf{a}(t{)}^{\top }\mathbf{\gamma })\text{\hspace{1em}(22)}$$

其中，$k$ 是增长率，$m$ 是偏移量，$a_j(t)=\mathbb{I}(t\ge s_j)$，其中 $s_j$ 是第 $j$ 个变化点的时间，${\delta }_j\sim \mathrm{Laplace}(\tau )$ 是变化的幅度，${\gamma }_j=-s_j{\delta }_j$ 使函数连续。$\mathbf{\delta }$ 上的拉普拉斯先验确保了最大后验参数估计是稀疏的，因此跨变化点边界的差异通常为 0。