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

摘要

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

1 引言

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

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

下图是季节项：

下图是剩余项：

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

1.1 设计目标

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

2 STL的定义

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

2.1 Loess

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

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

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

2.3 内循环

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

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

2.4 外循环

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

2.5 季节项的滞后平滑

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

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