时间序列 R 07 时间序列分解 Time series decomposition

一个时间序列可以分解为多个模型的组合

1.1 时间序列的组成

1.1.1 时间序列组成模式

三种时间序列模式（不计剩余残差部分）
1. 趋势Tend ：比如线性趋势，先增加后降低的整体趋势
2. 季节性Seasonal ：以时间为固定周期，呈现循环的特性
3. 周期性Cyclic：在以不固定周期不断震荡，通常周期性至少持续2年
下图就是讲时间序列分解之后的结果，应该比较容易理解上面的定义

下图是周期性的表现之一：

每个周期的震荡持续了6-10年，整体没有什么明显的趋势，第一幅图中trend包含了trend与cyclic

1.1.2 时间序列分解

可以将时间序列按下式分解为上面的几种模式，分解方式有相加方式或者相乘的方式如下：

$$y_{t} = S_{t} + T_{t} + E_{t},$$
或者
$$y_{t} = S_{t} \times T_{t} \times E_{t}.$$
其中 $S_{t}$为季节性模式部分， $T_{t}$为包含趋势与周期性的模式部分， $E_{t}$为剩余残差部分
其各自的适用范围：
- 季节性和trend-cycle变动的量级不随时间变化，用加法模式
- 两者的波动随时间成比例，适用于乘法模式
- 乘法模式可以通过log变为加法模式： $$y_{t} = S_{t} \times T_{t} \times E_{t} \quad \Longleftrightarrow\quad \log y_{t} = \log S_{t} + \log T_{t} + \log E_{t}$$
还是看这幅图

注意右边有一个灰色的比例条，他表示同样尺寸的线段在不同图上的比例，上图可见残差的长度较小。
下面介绍几种图;
1. seasonal sub-series plots

将每个seasonal里面的子周期呈现
2. Seasonally adjusted data
不考虑seasonal的数据

1.2 Moving averages 移动平均

移动平均用于估计时间序列的趋势和周期

1.2.1 Moving average smoothing

简单来说就是用数值前后k个时间点的值的平均值来代替当时值

$$\hat{T}_{t} = \frac{1}{m} \sum_{j=-k}^k y_{t+j}$$

m=2*k+1
称为m -MA
要注意的是，平滑之后，会使前面k个值和后面的k个值变为null.下图红色是平滑之后的曲线

1.2.2 Moving averages of moving averages

简单来说就是在已经平滑的基础上再一次平滑，这样可以使不同远近的数据有不同的权重，例如下例先使用4-ma平滑再使用2-ma平滑，称为2x4-MA：

公式为：

$$\begin{align*} \hat{T}_{t} &= \frac{1}{2}\Big[\frac{1}{4} (y_{t-2}+y_{t-1}+y_{t}+y_{t+1}) + \frac{1}{4} (y_{t-1}+y_{t}+y_{t+1}+y_{t+2})\Big] \\ &= \frac{1}{8}y_{t-2}+\frac14y_{t-1}+\frac14y_{t}+\frac14y_{t+1}+\frac18y_{t+2}. \end{align*}$$
可见最后的结果是取相近5个值的加权平均数，
其特征为第一个和最后一个值的权重为1/（2m）,其他值为1/m
另外还需注意，也可使用其他组合但是奇数需要和奇数组合，偶数与偶数组合

1.2.3seasonal data处理

• 如果季节性周期是偶数m，使用2*m-MA来估计trend-cyclic
• 如果季节性周期是奇数m，使用m-MA来估计trend-cyclic
• 例如使用2*12-MA来平滑以月为季节性周期的数据，使用7-ma来平滑以天为周期的数据，其他的选择可能会使trend-cyclic受到seasonal的影响。

1.2.4 权重移动平滑

故名，将原来的moving average 的每一相前面加一个权重，2*4-MA也相当于加了权重，文中提到了Spencer和Henderson权重，常见权重系数如下图：

S = Spencer’s weighted mov­ing average
H = Henderson’s weighted mov­ing average

1.3 经典分解

1.3.1 分解方法

经典分解分为加法分解和乘法分解
假设所用数据为季节性周期为m的时间序列,m也成为季节性指数
其步骤如下：
1. 如果季节性周期是偶数m，使用2*m-MA来估计trend-cyclic $\hat{T}_t$如果季节性周期是奇数m，使用m-MA来估计trend-cyclic $\hat{T}_t$
2. 计算去趋势序列$y-\hat{T}_t$
3. 计算季节趋势中每一个元素的平均值，如以一年中12月为周期，那么计算所有的1月的平均值作为1月的值，计算12个月的作为周期的$\hat{S}_t$
4. 计算$\hat{E}_t = y_t -\hat{T}_t - \hat{S}_t$
乘法未解与上面的模式类似
只不过是将减法变为除法，不再赘述

1.3.2 经典分解的缺点

1. 使用moving average时会使数据的前面和后面的数值缺失，变为null
2. 季节性的趋势可能会随时间变化，而经典中使用了平均值，忽略了这一变化
3. 其求平均值时，容易受到异常值的影响，如有一个月特别不正常，直接计算平均值容易得到不正确的结果

1.4 X-12-ARIMA

这是美国人口普查局的发明，克服了1.3.3中的缺点，但是只适用于季度或者月份为季节性周期的时间序列。
以乘法模式为例,分解步骤：
1. 以2*12-MA来计算$\hat{T}_t$
2. 计算集中率（centered rations）$y_t/\hat{T}_t$
3. 使用3*3-MA对$y_t/\hat{T}_t$的月份进行平滑以估计季节性$\hat{S}_t$
4. 用$y_t/\hat{T}_t$除以$\hat{S}_t$得到$\hat{E}_t$
5. 去掉不合理的$\hat{E}_t$，得到修正后 $\hat{E}_t$
6. 使用$\hat{E}_t$乘以$\hat{S}_t$得到修正后的centered rations
7. 重复步骤3得到修正后的$\hat{S}_t$
8. 用新的seasonal除源数据得到seasonally adjusted序列（$y_t/\hat{S}_t$）
9. 使用Henderson MA对上面的seasonal adjusted 序列平滑得到trend-cyclic $\hat{T}_t$，注意随机性越大使用的平滑次数越大，对于以月份为周期的序列可以使用9-, 13-, or 23-term Henderson
10. 重复步骤二得到新的$y_t/\hat{T}_t$centered rations
11. 使用新的数据，重复步骤3-6，其中使用3x5-MA代替3x3-MA
12. 使用3x5-MA重复步骤7
13. 重复步骤8得到seasonally adjusted序列（$y_t/\hat{S}_t$）
14. 上面得到的seasonally adjusted序列（$y_t/\hat{S}_t$）除以第9步得到的trend-cyclic $\hat{T}_t$，得到remainder component$\hat{E}_t$
15. 去除$\hat{E}_t$不合理的值
16. 完成
注意：是具体情况步骤11,12中的ma可以替换为3x3 3x5 3x9
R 中的相关包 https://cran.r-project.org/web/packages/x12/index.html

1.5 STL

stl=Seasonal and Trend decomposition using Loess
Loess是一种评估非线性关系的方法
1990年Cleveland发表的文章阐述了stl
http://cs.wellesley.edu/~cs315/Papers/stl%20statistical%20model.pdf

1.5.1 优点缺点

优点
1. 可以适用于任何季节周期的数据，（x12ARIMA只能是季度或者月份）
2. 季节性成分随时间的变化率可以被控制
3. trend-cycle 的平滑可以被控制
4. 更加鲁棒
缺点
1. 不能自动处理trading day or calendar variation
2. 只有加法模式
可以使用Box-Cox 来将乘法变为加法，将$\lambda$设为0为乘法，1为加法。
看一下例子


图6.10中将季节性变化设置为不变，得到的趋势更多变2010年之后成下降趋势
6.10的代码：

fit <- stl(elecequip, t.window=15, s.window="periodic", robust=TRUE)
plot(fit)

t.window, s.window是两个最终要的参数

1.6 分解预测

利用时间序列分解可以将时间序列分解为

$$y_t = \hat{S}_t + \hat{A}_t \\ \hat{A}_t = \hat{T}_t+\hat{E}_t$$
对于季节性的成分的预测直接使用了naive方法，具体见下面的链接
http://blog.csdn.net/bea_tree/article/details/51182002#t5，对于 $\hat{A}_t$的预测则有很多种方法，以后会慢慢叙述，这里还是使用了naive的方法。
栗子：
步骤1 分解，用naive来预测 $\hat{A}_t$

fit <- stl(elecequip, t.window=15, s.window="periodic", robust=TRUE)
eeadj <- seasadj(fit)
plot(naive(eeadj), xlab="New orders index",
  main="Naive forecasts of seasonally adjusted data")

,步骤2 用naive预测seasonal

fcast <- forecast(fit, method="naive")
plot(fcast, ylab="New orders index")

Further Reading
https://www.otexts.org/fpp/6/8

1.7 Further reading

Cleveland, R. B., W. S. Cleveland, J. E. McRae and I. J. Terpenning (1990). STL : A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics 6(1), 3–73.
Gomez, V. and A. Maravall (2001). “Seasonal adjustment and signal extraction in economic time series”. In: A course in time series analysis. Ed. by D. Pena, G.C. Tiao and R.S. Tsay. New York: John Wiley & Sons. Chap. 8, pp.202–246.
Ladiray, D. and B. Quenneville (2001). Seasonal adjustment with the X-11 method. Lecture notes in statistics. Springer-Verlag.
Miller, D. M. and D. Williams (2003). Shrinkage estimators of time series seasonal factors and their effect on forecasting accuracy. International Journal of Forecasting, 19(4), 669–684.
Theodosiou, M. (2011). Forecasting monthly and quarterly time series using STL decomposition. International Journal of Forecasting, 27(4), 1178–1195.