8 指数平滑

8.1 简单指数平滑

指数平滑最早是来自一种简单的预测方法： 用历史数据的线性组合预测下一时间点的值， 线性组合系数随距离变远而按负指数（几何级数）衰减：

\[ \hat x_h(1) \approx w x_h + w^2 x_{h-1} + \dots = \sum_{j=1}^\infty w^j x_{h+1-j} \] 其中\(0<w<1\)， \(w\)越小， 距离远的历史观测对预测的贡献越小。

因为是加权平均， 所以所有加权的和应该等于零， 注意到 \[ \sum_{j=1}^\infty w^j = \frac{w}{1-w} \] 所以第\(j\)个权重应为 \[ \frac{w^j}{\frac{w}{1-w}} = (1-w) w^{j-1}, \ j=1,2,\dots \] 于是 \[\begin{equation} \hat x_h(1) = (1-w)(x_h + w x_{h-1} + w^2 x_{h-2} + \dots) = (1-w) \sum_{j=0}^\infty w^j x_{h-j} \tag{8.1} \end{equation}\] 这种预测方法叫做指数平滑方法(exponential smoothing method)， 在早期应用中权重\(w\)是凭经验选取的。

经研究， ARIMA(0,1,1)的\(\hat x_h(1)\)的公式恰好具有如上形式。 对如下的ARIMA(0,1,1)模型： \[ (1-B) X_t = (1 - \theta B) \varepsilon_t \] 两边除以\(1 - \theta B\)，由于 \[ \frac{1}{1 - \theta z} = \sum_{j=0}^\infty \theta^j z^j \] 所以模型可以化为 \[ \begin{aligned} & (1 - B) \sum_{j=0}^\infty \theta^j X_{t-j} = \varepsilon_t \\ & \sum_{j=0}^\infty \theta^j X_{t-j} - \sum_{j=0}^\infty \theta^j X_{t-j-1} = \varepsilon_t \\ & X_t + \sum_{j=1}^\infty \theta^j X_{t-j} - \sum_{j=1}^\infty \theta^{j-1} X_{t-j} = \varepsilon_t \\ & X_t - (1-\theta) \sum_{j=1}^\infty \theta^{j-1} X_{t-j} = \varepsilon_t \\ & X_t - (1-\theta) \sum_{j=0}^\infty \theta^j X_{t-1-j} = \varepsilon_t \end{aligned} \] 以\(t=h+1\)代入得 \[ X_{h+1} = (1-\theta) \sum_{j=0}^\infty \theta^j X_{h-j} + \varepsilon_{h+1} \] 于是 \[ \hat x_h(1) = E(X_{h+1} | \mathscr F_h) = (1-\theta) \sum_{j=0}^\infty \theta^j X_{h-j} \] 这就是\(w=\theta\)的指数平滑预测公式。

因为指数平滑预测与ARIMA(0,1,1)预测的等价性， 可以用ARIMA建模方法估计权重\(w\)。 注意arima()函数估计的MA系数是\(1 + \theta_1 z + \dots + \theta_q z^q\)形式的系数。 也可以用ARIMA模型的充分性检验来判断指数平滑预测是否有意义。

例8.1 考虑CBOE的波动率指数（VIX）2004-01-02到2011-11-21的日收盘价的对数序列， 用指数平滑方法作一步预测。 共1988个观测。

读入数据，从中计算日收盘价对数值序列：

da <- read_table("d-vix0411.txt", col_types=cols(.default = col_double()))
xts.vix <- xts(da[,-(1:3)], make_date(da$year, da$mon, da$day))
str(xts.vix)

## An 'xts' object on 2004-01-02/2011-11-21 containing:
##   Data: num [1:1988, 1:4] 18 18.4 17.7 16.7 15.4 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:4] "Open" "High" "Low" "Close"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
##  NULL

vix <- log(xts.vix[,"Close"])
delta.vix <- diff(vix)[-1]

作日对数收盘价的时间序列图：

plot(vix, type="l", main="VIX Daily Log Close",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="auto")

图8.1: VIX日对数收盘价

明显不平稳且没有固定趋势。

作日对数收益率的时间序列图：

plot(delta.vix, type="l", main="VIX Daily Log Return",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="auto")

图8.2: VIX日对数收益率

日对数收益率的ACF:

forecast::Acf(c(coredata(delta.vix)), main="")

图8.3: VIX日对数收益率的ACF

虽然在滞后10的位置ACF也超界， 但是在低阶就只有滞后1明显超界， 所以对数收盘价用ARIMA(0,1,1)还是合理的。 ARIMA建模：

resm <- arima(c(coredata(vix)), order=c(0,1,1)); resm

## 
## Call:
## arima(x = c(coredata(vix)), order = c(0, 1, 1))
## 
## Coefficients:
##           ma1
##       -0.1500
## s.e.   0.0244
## 
## sigma^2 estimated as 0.004561:  log likelihood = 2535.71,  aic = -5067.42

建立的模型为 \[ X_t = X_{t-1} + \varepsilon_t - 0.1500 \varepsilon_{t-1} \] 其中\(\hat\sigma^2 = 0.004561\)。 \(\theta = 0.1500\)。

检验残差是否白噪声：

Box.test(resm$residuals, lag=10, fitdf=1)

## 
##  Box-Pierce test
## 
## data:  resm$residuals
## X-squared = 47, df = 9, p-value = 3.925e-07

白噪声检验很显著， 说明模型有所不足。 从差分序列的\(\hat\rho_{10}\)很大可以预见这个问题。

另行试验ARIMA(0,1,10):

resm2 <- arima(c(coredata(vix)), order=c(0,1,10)); resm2

## 
## Call:
## arima(x = c(coredata(vix)), order = c(0, 1, 10))
## 
## Coefficients:
##           ma1      ma2      ma3      ma4      ma5      ma6      ma7      ma8
##       -0.1493  -0.0716  -0.0484  -0.0258  -0.0303  -0.0401  -0.0131  -0.0181
## s.e.   0.0223   0.0226   0.0228   0.0225   0.0226   0.0236   0.0231   0.0219
##          ma9    ma10
##       0.0092  0.0966
## s.e.  0.0245  0.0230
## 
## sigma^2 estimated as 0.004453:  log likelihood = 2559.57,  aic = -5097.13

Box.test(resm2$residuals, lag=20, fitdf=10)

## 
##  Box-Pierce test
## 
## data:  resm2$residuals
## X-squared = 9.5782, df = 10, p-value = 0.4782

这个模型的白噪声检验通过，而且AIC也更好。

注意我们用了8年的数据， 教材(Tsay 2013)用的数据是2008-05-01到2010-04-19为两年数据。 两年数据可能ARIMA(0,1,1)就足够。

○○○○○

8.2 更一般的平滑方法

简单指数平滑， 相当于只用一个变化的水平值\(\ell_t\)拟合数据\(y_t\)， \[\begin{equation} \ell_t = \alpha y_t + (1 - \alpha) \ell_{t-1}, \tag{8.2} \end{equation}\] 预测为 \[\begin{equation} \hat y_{t+h|t} = \ell_t . \tag{8.3} \end{equation}\] 易见 \[\begin{equation} \ell_t \approx \sum_{j=0}^\infty \alpha (1-\alpha)^j y_{t-j} . \tag{8.4} \end{equation}\] 相当于(8.1)中\(w = 1-\alpha\)。

可以增加一个线性的趋势增长项\(b_t\)，变成： \[\begin{equation} \begin{aligned} \ell_t =& \alpha y_t + (1 - \alpha) (\ell_{t-1} + b_{t-1}), \\ b_t =& \beta (\ell_t - \ell_{t-1}) + (1 - \beta) b_{t-1} , \\ \hat y_{t+h|t} =& \ell_t + h b_t. \end{aligned} \tag{8.5} \end{equation}\]

可以增加季节项\(s_t\)，变成： \[\begin{equation} \begin{aligned} \ell_t =& \alpha (y_t - s_{t-m}) + (1 - \alpha) (\ell_{t-1} + b_{t-1}), \\ b_t =& \beta (\ell_t - \ell_{t-1}) + (1 - \beta) b_{t-1} , \\ s_t =& \gamma (y_t - \ell_{t-1} - b_{t-1}) + (1 - \gamma) s_{t-m}, \\ \hat y_{t+h|t} =& \ell_t + h b_t + s_{t - m + h_m^+} . \end{aligned} \tag{8.6} \end{equation}\] 其中\(m\)是季节频率，如月度数据为12， 季度数据为4。 \(t - m + h_m^+\)是\(t, t-1, \dots, t-m+1\)中与\(t-m+h\)同月（季度）的值。

各项可以是相乘的， 或者趋势水平项与季节性相乘。 对于带有指数增长模式的序列， 可以先作对数变换再建立加性模型。

8.3 时间序列成分分解

一个时间序列\(y_t\)， 常常将其分解为趋势部分\(T_t\)， 季节性变换\(S_t\)， 商业周期性变化\(C_t\)， 以及随机项\(R_t\)。 可以是加性的， 即 \[ Y_t = T_t + S_t + C_t + R_t . \] 也可以是乘性的， 如 \[ Y_t = T_t \times S_t \times C_t \times R_t , \] 或 \[ Y_t = T_t \times S_t \times C_t + R_t . \] 前面用过的航空乘客数的例子如果分解， 就应该使用乘性的模型。 乘性模型的数据常常可以通过进行对数变换， 转换为加性模型。

季节项\(S_t\)是有固定变换频率的成分， 比如月度数据， 相同月份的\(Y_t\)一般有很大的相似之处。 \(C_t\)则是更长期的周期更长的变化， 而且一般没有固定周期。 实际建模时常常将\(T_t\)与\(C_t\)合并在一起， 仅使用\(T_t\)项。

支持时间序列分解模型的R函数：

• stats包的HoltWinters()函数提供这样的指数平滑方法。
• forecast包的ets()函数提供了自动选择合适的指数平滑方法并进行预报的功能。
• stats包的StructTS()使用最大似然估计方法估计这样的分解， 仅支持加性模型。
• stats包的stl()用局部多项式拟合方法进行分解。
• stats包的decompose()用简单的滑动平均方法估计分解成分。
• feasts包和fabletools包的model和STL()函数。

考虑著名的航空数据的建模。 AirPassengers时间序列为1949年到1960年某航空公司月乘客人数。

图形：

plot(AirPassengers)

数据包含趋势、季节项， 并且应为乘性模型， 或者取对数后建模。

8.3.1 简单的滤波方法分解

最简单的一种分解方法是使用滑动平均估计趋势项， 去掉趋势后同月份平均估计季节项， 然后估计随机项。 这种方法得到的季节项是不变的模式。

较早期的的stats::decompose()就是用滑动平均滤波的方法估计趋势。

al_dec <- decompose(AirPassengers, 
  type = "multiplicative")
plot(al_dec)

fabletools包的model()与feasts包的classical_decomposition()配合可以进行这种分解。

比如，进行对数变换用这种方法分解：

library(tsibble)
library(feasts)
library(fabletools)
as_tsibble(log(AirPassengers)) |>
  model(classical_decomposition(
    value, type = "additive") ) |>
  components() |>
  autoplot() +
  labs(title = "航空乘客数对数值分解")

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).

图中各部分依次是原始数据、趋势、季节项、随机项。 图中左侧的条形用来比较不同坐标系中y轴范围大小， 这些条形在坐标系中按y坐标计算的高度是相同的， 但y轴范围越大， 这些条形就越短。

因为用了对称的滑动平均估计趋势项， 所以开头和结尾的6项无法估计。 这种方法得到的季节项是完全重复的。 趋势估计有可能过于平滑， 对于某些局部的变化估计结果不够稳健。

不做对数变换， 用乘性的移动平均进行分解：

as_tsibble(AirPassengers) |>
  model(classical_decomposition(
    value, type = "multiplicative") ) |>
  components() |>
  autoplot() +
  labs(title = "航空乘客数乘性分解")

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).

乘性的分解在用滑动平均估计了趋势项后， 用除法去掉趋势， 然后再对同月的数据平均得到季节项， 用除法得到剩余的随机项。

8.3.2 HoltWinters函数分解

Holt-Winters方法是利用8.2这样的递推算法估计水平、季节、随机项， 进行分解。

用R的HoltWinters()建模分解，并作拟合图：

al_hw <- HoltWinters(
  AirPassengers,
  seasonal = "multiplicative")
plot(al_hw)

各个分解成分的滤波结果：

plot(al_hw$fitted, main="分解成分滤波结果")

各成分为拟合值、\(l_t\)、\(b_t\)、\(s_t\)的滤波估计。

8.3.3 用forecast::ets()分解

这个函数是将指数平滑分解方法表达为状态空间模型进行估计， 见28.1, 29.3。 forecast::ets()函数能够自动选择需要的分解类型。

(al_ets <- forecast::ets(AirPassengers))

## ETS(M,Ad,M) 
## 
## Call:
##  forecast::ets(y = AirPassengers) 
## 
##   Smoothing parameters:
##     alpha = 0.7096 
##     beta  = 0.0204 
##     gamma = 1e-04 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 120.9939 
##     b = 1.7705 
##     s = 0.8944 0.7993 0.9217 1.0592 1.2203 1.2318
##            1.1105 0.9786 0.9804 1.011 0.8869 0.9059
## 
##   sigma:  0.0392
## 
##      AIC     AICc      BIC 
## 1395.166 1400.638 1448.623

plot(al_ets)

这个模型中误差项和季节项是乘性的， 趋势项是加性的。

8.3.4 用StructTS()分解

这也是用状态空间模型表示趋势、季节项、随机项的分解， 作最大似然估计。 不支持乘性模型， 可以预先取对数。

(al_sts <- StructTS(log10(AirPassengers)))

## 
## Call:
## StructTS(x = log10(AirPassengers))
## 
## Variances:
##     level      slope       seas    epsilon  
## 0.0001456  0.0000000  0.0002635  0.0000000

plot(al_sts$fitted)

8.3.5 用STL方法

STL方法用loess， 即局部多项式平滑方法进行分解。 STL能够处理多种季节频率， 季节项可以随时间缓慢变化， 趋势项的光滑程度可以用参数调节， 可以进行稳健的分解。 但是，STL不提供工作日和节假日调整， 仅能进行加性的分解。

stats::stl()是较早期的STL实现， feasts包的STL()与fabletools包的model()配合是较近期提供的STL实现。

使用stats::stl()的例子：

al_stl <- stl(AirPassengers, s.window=7)
plot(al_stl)

使用feasts::STL()的例子：

as_tsibble(AirPassengers) |>
  model(stl = STL(value)) -> ap_stl

components(ap_stl)中包含了原始数据和分解结果，如：

head(components(ap_stl))

## # A dable: 6 x 7 [1M]
## # Key:     .model [1]
## # :        value = trend + season_year + remainder
##   .model    index value trend season_year remainder season_adjust
##   <chr>     <mth> <dbl> <dbl>       <dbl>     <dbl>         <dbl>
## 1 stl    1949 1月   112  123.      -17.0      5.83           129.
## 2 stl    1949 2月   118  124.      -17.2     11.5            135.
## 3 stl    1949 3月   132  124.        7.52     0.174          124.
## 4 stl    1949 4月   129  125.       -1.15     5.27           130.
## 5 stl    1949 5月   121  126.       -3.78    -0.822          125.
## 6 stl    1949 6月   135  126.       18.9    -10.2            116.

分解结果中包含趋势、季节项、随机项， 还包含季节调整结果。

可以用autoplot()对分解结果作图：

components(ap_stl) |>
  autoplot()

这是加性模型， 所以季节项的振幅随时间增加而变大。

原始数据与季节调整数据的叠加图形：

components(ap_stl) |>
  as_tsibble() |>
  autoplot(value, colour = "gray") +
  geom_line(aes(y=season_adjust), colour = "#0072B2") +
  labs(y = "乘客数 (千人)",
    title = "航空乘客季节调整")

因为使用了加性模型而非乘性模型， 所以季节调整的效果不佳。

feasts包通过比较趋势项、季节项、随机项的方差， 定义了趋势强度和季节性强度指标， 都取值于0和1之间。 如：

as_tsibble(AirPassengers) |>
  features(value, feat_stl) |>
  pivot_longer(everything(),
    names_to="Strength",
    values_to="Value")

## # A tibble: 9 × 2
##   Strength                  Value
##   <chr>                     <dbl>
## 1 trend_strength            0.991
## 2 seasonal_strength_year    0.941
## 3 seasonal_peak_year        7    
## 4 seasonal_trough_year     11    
## 5 spikiness                 3.03 
## 6 linearity              1325.   
## 7 curvature               131.   
## 8 stl_e_acf1                0.509
## 9 stl_e_acf10               0.930

8.3.6 用X11和SEATS季节调整方法分解

X11和SEATS都是美国、澳大利亚等国官方统计用于对月度、季度数据进行季节调整的算法。 这些算法经过多年研究改进， 已经比较成熟可靠。 下面的程序利用了R的seasonal包进行计算。

X11是一系列的滤波算法， 也能够适应工作日的变化、节日效应等。

library(seasonal)

## Warning: package 'seasonal' was built under R version 4.4.3

## 
## Attaching package: 'seasonal'

## The following object is masked from 'package:tibble':
## 
##     view

as_tsibble(AirPassengers) |>
  model(x11 = X_13ARIMA_SEATS(value ~ x11())) |>
  components() -> x11_ap
autoplot(x11_ap) + 
  labs(title = "航空乘客数X11分解")

对季节调整作图：

library(seasonal)
x11_ap |>
  ggplot(aes(x = index)) +
  geom_line(aes(y = value, colour = "Data")) +
  geom_line(aes(y = season_adjust,
                colour = "Seasonally Adjusted")) +
  geom_line(aes(y = trend, colour = "Trend")) +
  labs(y = "乘客数 (千人)", x = "",
       title = "航空乘客数")

X11得到的季节项是缓慢变化的， 分月作图：

x11_ap |>
  gg_subseries(seasonal)

SEATS是另外一套季节调整算法， 由西班牙统计机构开发。

as_tsibble(AirPassengers) |>
  model(seats = X_13ARIMA_SEATS(value ~ seats())) |>
  components() -> seats_ap
autoplot(seats_ap) +
  labs(title = "航空乘客数用SEATS分解")

B 参考文献

———. 2013. 金融数据分析导论：基于R语言. 机械工业出版社.