时间序列分解论文STL: A Seasonal-Trend Decomposition Procedure Based on Loess

摘要

  STL是一种滤波方式可以把时间序列分解为趋势项、季节项和剩余项。STL包含一系列局部加权回归平滑器，计算速度比较快，可以应对非常大的的时间序列数据。

1 引言

  STL是一个可以把一个季节性时间序列分解为三项：趋势项、季节性和剩余项。
  下图是原始序列：

  下图是趋势项：数据中非平稳的长期变化的低频变量

  下图是季节项：

  下图是剩余项：

  分别用$Y_v,T_v,S_v$和$S_v$表示时间序列数据、趋势项、季节项和剩余项，所以有：$$Y_v=T_v+S_v+R_v$$
  以上四张图片就是文中的figure1，它是一个二氧化碳数据，横坐标以天为单位，共有4193个数据。

1.1 设计目标

  本文的第二部分介绍loess（局部移动平均）
  本文的第三部分介绍如何选择STL中的参数。有些参数必须基于一些必要的先验知识，而有些参数的选择需要基于时间序列的性质，第三部分还给出了一些选择参数的辅助方法。
  本文的第四部分介绍了STL方法的实现。
  本文的第五部分阶介绍了特征值和频率分析。
  本文的第六部分是全文的总结部分。

2 STL的定义

  这一部分介绍loess平滑器和STL的操作。目标是针对STL给出简单直接的解释。

2.1 Loess

  假定$x_i$和$y_i$分别是自变量和因变量，$g(x)$是局部加权回归曲线。实际上，loess可以作为函数去平滑任意数量的自变量，但是对于STL，只需要针对一个自变量的情况。
  $g(x)$的计算步骤如下：选定一个正整数$q$，暂时假定$q\le n$，最接近$x$的$q$个值被选中，并根据它们距离$x$的距离远近为其赋予一个邻接权重。 用${\lambda }_q(x)$表示选中的$q$个值中距离$x$最远的距离值。用$W$表示三次权重函数：$$for0\le u<1,W(u)=(1-u^3{)}^3$$,$$foru\ge 1,W(u)=0$$
所以对于任意的$x_i$得到邻接权重：$$v_i(x)=W(\frac{|x_i-x|}{{\lambda }_q(x)})$$
  所以距离$x$越近的$x_i$具有最大的邻接权重，距离$x$越远的$x_i$的权重越小，直到超过了最远的距离权重就变成了零。
   接下来，用一个$d$阶的多项式去拟合这些带有权重$v_i(x)$的数据。 这个局部拟合多项式在$x$处的值是$g(x)$。在这篇文章中，$d=1or2$，也就是局部线性或局部二次曲线。
  现在假设$q>n$。${\lambda }_n(x)$是从$x$到最远的$x_i$之间的距离。因为$q>n$，所以定义${\lambda }_q(x)$为：$${\lambda }_q(x)={\lambda }_n(x)\frac{q}{n}$$，然后我们用前面一样的方式得到邻接权重。
  随着$q$的增加，$g(x)$会变得平滑。当$q$趋向于无穷，权重$v_i(x)$趋近于1并且$g(x)$趋近于一个一般的$d$阶的最小二乘多项式。当$d=1$时，表示数据潜在的模式有非常小的曲率；当$d=2$时，表示数据存在非常显著的曲率，例如波峰、波谷。
  假定给定的每一个观测值$(x_i,y_i)$都有一个权重${\rho }_i$来表示这个观测值相比于其他观测值的可信度。举个例子，如果$y_i$的方差为${\sigma }^2{k}_i$且$k_i$已知，所以${\rho }_i$可能为$1/k_i$。我们可以用一种简单直接的方式把这些权重合并到局部加权回归中，其中用${\rho }_i{v}_i(x)$作为局部最小二乘的权重。正如我们看到的，这提供了一种我么可以向STL添加鲁棒性的机制。

2.2 整体设计：内循环和外循环

  STL包含有两个循环的过程：外循环中嵌套一个内循环。内循环每迭代一次，季节项和趋势项就更新一次；完整的内循环包括$n_{(i)}$次这样的迭代。外循环每迭代一次包括了内循环和鲁棒性权重的计算；这些权重在下一次内循环中将会被用来减少暂时的、离奇的点对趋势项和季节项的影响。最初的外循环迭代的时候所有的鲁棒性权重等于1，然后执行$n_{(o)}$次外循环。在第三部分，将会讨论$n_{(i)}$和$n_{(o)}$的选取方式。
  假定季节项中每一个周期或循环包含的观测值的数量为$n_{(p)}$。例如，如果时间序列是关于月份的每年为一个周期，所以$n_{(p)}=12$。我们称$n_{(p)}$的一个周期为周期子序列。

2.3 内循环

  内循环每迭代一次包括一次季节项平滑更新季节项，紧接着趋势项平滑更新趋势项。假定对于$v=1toN$，$S_v^{(k)}$和$T_v^{(k)}$分别是第$k$次迭代后的季节项和趋势项；这两项的定义是在所有的时间上$v=1toN$，即使某个时间点上的数据缺失。第$k+1$步，$S_v^{(k+1)}$和$T_v^{(k+1)}$通过以下的方式进行计算：

• 步骤1：去趋势。计算去趋势序列$Y_v-T_v^{(k)}$。如果$Y_v$在这一点缺失，去趋势序列在这一点也缺失。
• 步骤2：周期子序列平滑。去趋势序列中的每一个周期子序列都是通过参数为$q=n_{(s)}$和$d=1$的局部加权回归进行平滑。所有时间点的平滑之都要进行计算，包括那些缺失的值和时间序列第一个点前面的一个值以及时间序列的最后一个点后面一个值。例如，假如有一个时间序列的时间是按月的，$n_{(p)}=12$，该序列是从1943年的1月到1985年的1月，并且1960年一月的数据缺失；然后平滑就是计算从1942年1月到1986年1月的所有平滑值。所有的周期子序列的平滑值的集合就是一个暂时的季节序列$C_v^{(k+1)}$，包括$N+2n_{(p)}$个值，从$v=-n_{(p)}+1$到$N+n_{(p)}$。figure1中的$n_{(s)}=35$，具体的$n_{(s)}$的选择方式将会在第三部分和第五部分讨论。
• 步骤3：平滑之后的周期子序列的低通滤波器。这个滤波器是应用于$C_v^{(k+1)}$的。这个滤波器包括一个长度为$n_{(p)}$移动平均，一个长度为3的移动平均，一个参数$d=1$和$q=n_{(l)}$的局部移动平均平滑。$n_{(l)}$的选择方式将会在第三部分和第五部分讨论。对于figure1来说，$n_{(l)}=365$。最后的输出$L_v^{(k+1)}$是定义在从$v=1$到N上。
• 步骤4：平滑之后的周期子序列去趋势。季节项中的第$(k+1)$次循环为$S_v^{(k+1)}=C_v^{(k+1)}-L_v^{(k+1)}$，对所有的$v=1$到$N$。这里的$L_v^{(k+1)}$被减去防止低频信息进入季节项。
• 步骤五：去季节项。计算去季节项序列$Y_v-S_v^{(k+1)}$。如果$Y_v$在特定的时间点有缺失值，则这个去季节项也缺失。
• 步骤六：趋势平滑。去季节项通过参数为$q=n_{(t)}$和$d=1$的局部加权回归进行平滑。从$v=1$到$N$在所有的时间点上都计算，即使这个点有缺失值。在第$k+1$次循环以后得到趋势项$T_v^{(k+1)}$即平滑项的集合。对于figure1来说，$n_{(t)}=573$，$n_{(t)}$的选择方法将会在第三部分和第五部分讨论。
    因此，内循环中的季节平滑项就是步骤2，3和4，趋势平滑项是步骤6。在执行之前，设趋势项的初值为0，即$T_v^{(0)}=0$。

2.4 外循环

  假设我们已经进行了内部循环的初始运行以获得估计值：$T_v$和$S_v$分别是趋势项和季节项，然后我们得到了季节项:$$R_v=Y_v-T_v-S_v$$我们会为每个观测值$Y_t$定义一个权重。这些鲁棒权重反映了这些残差值的有多极端。残差值非常大的奇点将会有一个非常小的或零的权重。让：$$h=6median(|R_v|)$$，然后就可以得到时间点$v$处的鲁棒性权重$${\rho }_v=B(|R_v|/h)$$其中B是二次平方权重函数：$$for0\le u<1,B(u)=(1-u^2{)}^2$$$$foru>1,B(u)=0$$
  在步骤2和步骤6的平滑操作中，时刻$v$的邻接权重会乘以一个鲁棒权重${\rho }_v$。这这种外循环上的鲁棒操作会执行总共$n_{(o)}$次。除了初始步骤以外，每次我们进入到内循环中，我们并没有设置$T_v^{(0)}=0$，而是使用前一次内循环中步骤6的趋势项。

2.5 季节项的滞后平滑

  考虑figure1中每天的$C{O}_2$序列。STL的步骤2，也就是构成季节项的最基础的一步有如下的形式：这365个去趋势序列被分别平滑，并且它们被放在一起构成了暂时的季节项。这也就意味着季节项中的每一个周期子序列将会被平滑在整个年的长的。例如，7月4日的季节项的值从一年到零一年变化得非常平滑。但是这种平滑并没有保证从一天到另一天的季节项变化也是平滑的。

  figure2的上面的图展示了每天$C{O}_2$序列的季节项，很明显局部非常粗糙。但是对于这个时间序列来说，我们想要得到一个项，这个项的前一天到后一天是平滑的。
  一个简单的平滑方式就是滞后平滑。平滑的值作为季节项的最终值。在进行此平滑操作时，我们希望确保使用局部二次拟合因为季节项中有很多波峰和波谷。最后，$q$很小或中等因为这些不平滑通常具有很小的方差。对于每天的$C{O}_2$数据，$q$取51。最终的结果就是figure2的下图，也就是figure1中展示的项。