时间序列分析:三重指数平滑
双重指数平滑引入了趋势,但还不能处理序列中存在的季节性。为了体现序列中的趋势,对 $x_{n+1}$ 的预测 $\hat x_{n+1}$ 应由局部水平 $L_n$、趋势 $T_n$ 以及季节项 $S_{n+1-m}$ 组成,其中 $m$ 为季节周期。
\begin{equation}
\hat x_{n+1} = L_n + T_n + S_{n+1-m} \tag{1}
\end{equation}
其中局部水平 $L_n$ 需要移除季节项的干扰,是对 $x_n – S_{n-m}$ 的指数平滑,即
\begin{equation}
L_n = \alpha (x_n – S_{n-m}) + (1 – \alpha)(\hat x_n – S_{n-m}) \tag{2}
\end{equation}
其中
\begin{equation}
\hat x_{n} = L_{n-1} + T_{n-1} + S_{n-m}
\end{equation}
于是
\begin{align}
L_n = \alpha (x_n – S_{n-m}) + (1 – \alpha)(L_{n-1} + T_{n-1}) \tag{3}
\end{align}
趋势项 $T_n$ 的定义与双重指数平均相同,为
\begin{equation}
T_n = \beta \cdot (L_n – L_{n-1}) + (1 – \beta) \cdot T_{n-1} \tag{4}
\end{equation}
季节项 $S_n$ 是对季节的平滑,需要移除局部水平的干扰,即对 $x_n – L_n$ 进行指数平滑
\begin{equation}
S_n = \gamma (x_n – L_n) + (1 – \gamma )S_{n-m} \tag{5}
\end{equation}
式 $(1)$、$(3)$、$(4)$、$(5)$ 即构成了三重滑动平均模型。当使用前 $n$ 个数据对 $n+h$ 时刻进行预测时,有
\begin{equation}
\hat x_{n+h} = L_n + h \cdot T_n + S_{n+h-m} \tag{6}
\end{equation}
前面引入的季节性都是加性的,如果存在乘性的季节性,则式 $(1)$ 变为
\begin{equation}
\hat x_{n+1} = (L_n + T_n) \cdot S_{n+1-m} \tag{7}
\end{equation}
式 $(6) $ 变为
\begin{equation}
\hat x_{n+h} = (L_n + h \cdot T_n) \cdot S_{n+h-m} \tag{8}
\end{equation}
2. 一个例子
载入国际航空旅客人数的月度数据(单位为千人),绘制图像如图 1。
plot(AirPassengers)
可见该数据具有季节性,周期振幅(方差)随时间增长,具有乘性周期。
对数据进行 $\log$ 变换后,绘制图像如图 2。
data.ts.log <- log10(data.ts) plot(data.ts.log)
可见此时的季节性就变成了加性的了。
使用 HoltWinters
函数进行拟合
model <- HoltWinters(data.ts.log) model
得到的结果为
Holt-Winters exponential smoothing with trend and additive seasonal component. Call: HoltWinters(x = data.ts.log) Smoothing parameters: alpha: 0.326612 beta : 0.005744246 gamma: 0.8207255 Coefficients: [,1] a 2.680598830 b 0.003900787 s1 -0.031790733 s2 -0.061224237 s3 -0.015941495 s4 0.006307818 s5 0.014138008 s6 0.067260071 s7 0.127820295 s8 0.119893006 s9 0.038321663 s10 -0.014181699 s11 -0.085995400 s12 -0.044672707
手动计算 $h = 1$(1961 年 1 月)和 $h = 13$(1962 年 1 月)时的预测值
2.680598830 + 1 * 0.003900787 - 0.031790733 2.680598830 + 13 * 0.003900787 - 0.031790733
结果分别为 2.652708884
和 2.699518328
。
使用模型进行预测
library(forecast) model.forecast <- forecast(model) model.forecast
结果为
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 Jan 1961 2.652709 2.630898 2.674520 2.619351 2.686066 Feb 1961 2.627176 2.604218 2.650134 2.592065 2.662287 Mar 1961 2.676360 2.652297 2.700422 2.639560 2.713160 Apr 1961 2.702510 2.677380 2.727640 2.664077 2.740942 May 1961 2.714241 2.688076 2.740406 2.674225 2.754257 Jun 1961 2.771264 2.744092 2.798436 2.729708 2.812820 Jul 1961 2.835725 2.807571 2.863878 2.792667 2.878782 Aug 1961 2.831698 2.802586 2.860811 2.787174 2.876222 Sep 1961 2.754028 2.723977 2.784079 2.708069 2.799987 Oct 1961 2.705425 2.674454 2.736396 2.658059 2.752791 Nov 1961 2.637512 2.605638 2.669386 2.588765 2.686259 Dec 1961 2.682736 2.649974 2.715497 2.632631 2.732840 Jan 1962 2.699518 2.661306 2.737731 2.641078 2.757959 Feb 1962 2.673986 2.635014 2.712957 2.614383 2.733588 Mar 1962 2.723169 2.683445 2.762894 2.662416 2.783923 Apr 1962 2.749319 2.708848 2.789790 2.687424 2.811214 May 1962 2.761050 2.719838 2.802262 2.698022 2.824078 Jun 1962 2.818073 2.776126 2.860020 2.753921 2.882226 Jul 1962 2.882534 2.839857 2.925211 2.817265 2.947803 Aug 1962 2.878508 2.835105 2.921910 2.812129 2.944886 Sep 1962 2.800837 2.756714 2.844960 2.733356 2.868318 Oct 1962 2.752234 2.707395 2.797074 2.683658 2.820811 Nov 1962 2.684322 2.638770 2.729873 2.614656 2.753987 Dec 1962 2.729545 2.683285 2.775805 2.658796 2.800294
可见 Jan 1961
和 Jan 1962
的预测值与前面手工计算的结果一致。
绘制预测结果如图 3 所示。
plot(model.forecast)
注意这里的预测值都是经过 $\log$ 变换的,需要进行反变换得到真正的预测值。