时间序列分析:三重指数平滑

  双重指数平滑引入了趋势,但还不能处理序列中存在的季节性。为了体现序列中的趋势,对 $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)

图 1

可见该数据具有季节性,周期振幅(方差)随时间增长,具有乘性周期。

  对数据进行 $\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.6527088842.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 1961Jan 1962 的预测值与前面手工计算的结果一致。

  绘制预测结果如图 3 所示。

plot(model.forecast)

图 3

注意这里的预测值都是经过 $\log$ 变换的,需要进行反变换得到真正的预测值。