时间序列分解算法和d ecompose 函数实现
李思亮 55531469@qq.com
目录
时间序列分解算法和decompose 函数实现 ............................................................................................................. 1 1 数据读入并生成时间序列 ................................................................................................................................... 2 2 数据可视化 ........................................................................................................................................................... 4 3 时间序列分解 ....................................................................................................................................................... 7
在时间序列分析的过程中，往往需要对时间序列作出初步分析，本文主要采用R 语言作为分析平台，从数据的读入，可视化图，分解（decompose ）为趋势项，季节项，随机波动等角度对数据开展分析的几个案例。最后对分解算法作出初步描述并探讨其预测预报中的潜在应用。本文的数据和部分内容主要采用中的内容，感兴趣的读者可以参考。
1 数据读入并生成时间序列
对于数据分析来讲，数据读入是一个比较关键的步骤。常用的数据读入函数有scan ，read.table 等。下面列举了几种常见的数据。
首先是http://robjhyndman.com/tsdldata/misc/kings.dat，中包含了英国国王的寿命从William 开始，数据来源（Hipel and Mcleod, 1994）。
> kings kings
[1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56
上述例子中，读入了连续42个公国国王的寿命并将其赋给变量‘kings ’
如果我们希望对读入数据开展分析，下一步就是将其转化为时间序列对象（时间序列类），R 提供了很
多函数用于分析时间序列类数据。可以使用ts 函数将变量转化为时间序列类。
> kingsts
> kingsts Time Series: Start = 1 End = 42
Frequency = 1
[1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56
有时候我们会按照一定的时间周期来收集数据，这个周期可能是季度，月，日，小时，分。在大数据时代，有些情况下的数据是按照秒来采集收集。这种情况下，我们需要对数据的周期或频率进行设置。这里采用ts 函数中的frequency 参数可以实现这种功能。比方说，若按1年为一个周期，我们的月度时间
对于上述数据操作的好处是将数据转化为特定的“时间序列类”便于我们使用R 中的函数分析数据。
序列数据应为frequency=12，若为季度时间序列数据，则可设置frequency=4。
另外，还可以利用start 参数，设置时间序列的起点，比如若我们一个周期共4个观测，而第一个数据对应为1986年的起的第2个观测，则可使用start=c(1986,2)。
有一个纽约市月出生数量的数据集，从1946年1月至1959年12月。数据可通过
http://robjhyndman.com/tsdldata/data/nybirths.dat 获取。我们将其读入至R 中。 > births
> birthsTS
> birthsTS
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227 21.672 21.870
1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110 21.759 22.073
1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142 21.059 21.573
1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907 21.519 22.025
1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252 22.084 22.991
1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110 22.964 23.981
1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199 23.162 24.707
1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462 25.246 25.180
1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379 24.712 25.688
1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784 25.693 26.881
1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136 26.291 26.987
1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484 26.634 27.735
1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945 25.912 26.619
1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012 26.992 27.897
类似的，文件http://robjhyndman.com/tsdldata/data/fancy.dat 中包含了澳大利亚昆士兰州的一个海滩小镇纪念品商店的月销量数据，时间从1987/01~1993/12(数据来源Wheelwright and Hyndman, 1998)。我们用下述代码将其读入至R 中。
> souvenir
> souvenirTS
> souvenirTS
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74 4349.61 3566.34 5021.82 6423.48 7600.60 19756.21
1988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12 4752.15 5496.43 5835.10 12600.08 28541.72
1989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62 8176.62 8573.17 9690.50 15151.84 34061.01
1990 5921.10 5814.58 12421.25 6369.77 7609.12 7224.75 8121.22 7979.25 8093.06 8476.70 17914.66 30114.41
1991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55 12552.22 11637.39 13606.89 21822.11 45060.69
1992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78 19888.61 23933.38 25391.35 36024.80 80721.71
1993 10243.24 11266.88 21826.84 17357.33 15997.79 18601.53 26155.15 28586.52 30505.41 30821.33 46634.38 104660.67
2 数据可视化
上述章节中我们读入并生成了3个时间序列数据分别为，kingsts ，birthsTS 和souvenirTS 。数据分析的后续过程为可视化过程，让有关人员可以从直观角度来认识数据的变化规律。plot 函数提供了各种类的绘图功能
> plot(kingsts)
上图中的时间序列看来的话，没有比较明显的周期性规律，直观判断上述时间序列类似于一个随机波动序列（平稳时间序列），后续我们会探讨如何判断序列为平稳时间序列。
我们对纽约市月出生数量画图，结果为
> plot(birthsTS)
从上述结果看来，每个月的出生数量似乎含有一定的季节变化（周期性），在夏季较高，在冬季较低。另外每个周期内的波动幅度（振幅）似乎是固定的，并不随时间序列的水平（趋势）变化而变化。随机波动随时间变化也不明显。
我们同样对澳大利亚昆士兰州海滩旅游小镇纪念品商店的月销量数据画图。
从上图中可以看出误差波动和振幅都随时间水平增加而增加。因此，我们需要对原始时间序列作出对
数变换（取log ），以便使用additive 模型对数据开展分析预测。
直观上看来，上述的对数化时间序列的季节变化和随机变化随时间稳定在一个范围内，不随时间序列变化而变化。因此，时间序列对数化可以用于additive 模型。
3 时间序列分解
时间序列分解主要针对非平稳时间序列，非平稳时间序列包含4种主要趋势：长期趋势（Trend ），循环趋势，季节趋势（Seasonal ）和随机趋势（Random ）。
可以利用additive 模型分解长期趋势项和季节项。使用R 语言的decompose 函数。该函数可以分解出趋势项，季节项和随机项并利用additive 模型描述时间序列。
时间序列分解的主要参考的函数有decompose 函数和stl 函数。在这里我们先介绍decompose 函数的主要思路。
相加分解模型（additive ）的主要形式为：
x t =Tt +St +et
另外还有相乘模型（multiplicative ），主要形式为：
x t =Tt *St *et
其中，在t 时刻，x t 为该时刻的观测值，T t 为该时刻的趋势值，S t 为该时刻的季节（周期）效应，e t
为该时刻的无规则部分，也就是剩余项。一般而言，理想的分解模型结果中的剩余项应该为一个均值为0
的随机变量。
我们先利用R 的decompose 函数，看看有关时间序列分解的结果，再从结果入手分析其具体数据分析处理流程。分解我们之前的birthsTS 数据，可视化后的结果为：
> birthsTSDecompose plot(birthsTSDecompose)
Decomposition of additive time series
30
o b s e r v e d t r e n d s e a s o n a l r a n d o m
-1. 5
1946
-0. 50. 51.
5-2. 0-1. 00. 01. [***********]42628
[***********]5619581960
Time
decompose 函数提供了additive 和multiplicative 2种方法来估算各分解项。本文主要分析additive 模型的分解步骤。
现在的已知条件：时间序列，包含了时间序列的频率（frequency ），起始时间（start ）。
Step1：分解出趋势项。采用中心化移动平均值来计算趋势项，
x T t =
t -(
f -1
) 2
+x
t -(
f -1
) +12
+... +x
t +(
f -1
) 2
f
0. 5x
+x
+... +x f
, t ∈(
f +1f -1
, l -) 22
当f 为奇数时采用上述计算方法
f
t -()
2
f t -() +1
2
f t +() -1
2
+0. 5x
T t =
f t +()
2
, t ∈(
f f +1, l -) 22
当f 为偶数时采用上述计算方法，其中Tt 为趋势项，f 为时间序列频率，l 为时间序列长度。结果为长度为l 的时间序列，为便于后续的向量计算，当t 超出上述下标的定义域时，其值为NA ，如T 1。
R 中提供filter 函数计算2个序列的线性滤波（卷积），上述时间序列xt 的权重恰好为一个对称序列，因此可以通过线性滤波函数实现中心化移动平均，见如下代码： > f
> Ts identical(Ts,birthsTSDecompose$trend) [1] TRUE
Step2：计算季节项。
1、 采用将原始时间序列减去趋势项
S t =x t -T t
2、 将各个周期内相同频率下的值平均化，得到季节项figure
figure t =∑S t +i*f/f,t∈（1，f），n
i =0
n
=l 对f 取整，即max (n,nf≤l)
3、 将figure 中心化，得到中心化的季节项Figure ，代码可表述为
figure=figure-mean(figure)
最终得到的的长度为f 的季节项。
假设时间序列的频率为f ，长度为l, 设x t -T t 为S t ，则上述算法可以表示为: 用表格可以表示为
注:-为NA 数据，其中St 序列本身或已包含了NA 数据。
Step3：计算周期性季节项Seasonal ，周期性季节项的计算则为将figure 重复至长度为l 的序列。
Seasonal t =figure t %%f , 其中%%为取余运算
Step4：计算随机项Random ，随机项的计算方法为
Random =X t -T t -Seasonal t
4 总结和展望
根据上述的分析可以看出R 内置的decompose 函数的核心为趋势项（step1）和季节项（step2）。关于趋势项的计算采用了中心化移动平均的方法。这种方法计算Ts 时采用的x 权重系数都是相同的，可以考虑采用自适应滤波的方法将各个权重系数做实时调整。主要流程如下：
此外，若不采用自适应滤波法调整计算趋势项的x 各项权重。可以考虑对误差开展时间序列分析，进一步分解误差。这个的讨论已经超出了本文的主题，后续将另作研究。
关于基本模型的设计方面，decompose 采用了加法模型和乘法模型。另外可以考虑采用混合模型，如
x t =Tt *St +et
基于上述模型建模，不在本篇讨论。
最后我们可以将分解的时间序列趋势项，季节项或者误差项作出预测。如Holt Winters 算法和stl 算法均体现了时间序列分解在预测方面的应用。