作者：黄天元，复旦大学博士在读，热爱数据科学与开源工具（R），致力于利用数据科学迅速积累行业经验优势和科学知识发现，涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等，著有《R语言数据高效处理指南》（《R语言数据高效处理指南》(黄天元)【摘要 书评 试读】- 京东图书，《R语言数据高效处理指南》(黄天元)【简介_书评_在线阅读】 - 当当图书）。

简单指数平滑

指数平滑法是基于滑动平均法延伸而来的，滑动平均是把前N个序列的平均值作为N+1的预测，不断地往后计算，获得一条相对平滑的曲线。但是指数平滑法则认为，时间越靠近那么应该予以的权重就越大。赋予权重的大小一般以α表示，具体公式请参照7.1 Simple exponential smoothing | Forecasting: Principles and Practice，这里仅给出实现。

简单指数平滑的实现代码如下：



library(fpp2)
oildata <- window(oil, start=1996)
# Estimate parameters
fc <- ses(oildata, h=5)
# Accuracy of one-step-ahead training errors
round(accuracy(fc),2)
#>               ME  RMSE   MAE MPE MAPE MASE  ACF1
#> Training set 6.4 28.12 22.26 1.1 4.61 0.93 -0.03



其中，ses函数可以对后面的5个观测进行预测，但accuracy函数所给出的是在训练集中得到的准确度指标。可视化如下：



autoplot(fc) +
  autolayer(fitted(fc), series="Fitted") +
  ylab("Oil (millions of tonnes)") + xlab("Year")





r语言五日滑动平均法 r语言平滑处理_ci


拟合线和预测线都显示了出来，预测线还包含置信区间。这里对α进行了自动的估计，可以进行查看：


fc$model

Simple exponential smoothing 

Call:
 ses(y = oildata, h = 5) 

  Smoothing parameters:
    alpha = 0.8339 

  Initial states:
    l = 446.5868 

  sigma:  29.8282

     AIC     AICc      BIC 
178.1430 179.8573 180.8141


α越大，说明对新近观测赋予的权重越高。

Holt指数平滑（考虑趋势性）

简单指数平滑没有考虑序列的趋势性，因此引入了Holt指数平滑法：


air <- window(ausair, start=1990)
fc <- holt(air, h=5)
fc$model

Holt's method 

Call:
 holt(y = air, h = 5) 

  Smoothing parameters:
    alpha = 0.8302 
    beta  = 1e-04 

  Initial states:
    l = 15.5715 
    b = 2.1017 

  sigma:  2.3645

     AIC     AICc      BIC 
141.1291 143.9863 147.6083


可以看到这个方法比之前获得了一些新的参数，如beta和b。beta是趋势性的指数，它表示随着时间的变化，时间序列趋势所发生的变化。而b则是初始值。

该方法有一种变式，叫做Damped trend methods。它的假设是，一个序列在短期内是具有一定趋势的，但是在长期中却是围绕着一个均值在波动。这里引入一个新的参数ϕ，0<ϕ<1，一般在0.8到0.98之间，如果等于1，那么就与传统的Holt平滑法无异。实现如下：


fc <- holt(air, h=15)
fc2 <- holt(air, damped=TRUE, phi = 0.9, h=15)
autoplot(air) +
  autolayer(fc, series="Holt's method", PI=FALSE) +
  autolayer(fc2, series="Damped Holt's method", PI=FALSE) +
  ggtitle("Forecasts from Holt's method") + xlab("Year") +
  ylab("Air passengers in Australia (millions)") +
  guides(colour=guide_legend(title="Forecast"))


r语言五日滑动平均法 r语言平滑处理_r语言五日滑动平均法_02


可以看出来，红线可能更加合理，而蓝线很可能无止境地网上一直增长。观察一下模型参数：


fc2$model
Damped Holt's method 

Call:
 holt(y = air, h = 15, damped = TRUE, phi = 0.9) 

  Smoothing parameters:
    alpha = 0.763 
    beta  = 0.2526 
    phi   = 0.9 

  Initial states:
    l = 14.5766 
    b = 3.7517 

  sigma:  2.6165

     AIC     AICc      BIC 
145.3965 148.2536 151.8757


ϕ（读作phi）是用来削弱趋势性的参数，ϕ越大则对趋势削弱的效果越小，等于1的时候则完全没有削弱的效果。如果0<ϕ<1，那么随着预测后续变量越多，将会让预测趋近于ℓT+ϕbT/(1−ϕ)ℓT+ϕbT/(1−ϕ)。具体细节请看：7.2 Trend methods | Forecasting: Principles and Practice。

上面的预测中我们认为给定了phi，如果给出可以自动进行估计。

Holt-Winters’平滑法

考虑了趋势性之后，还可以继续考虑季节性因素，进而提出Holt-Winters’平滑法。其中，季节性可以分为可加和可乘两种，如果季节性随着时间变化而变化，那么应该使用可乘，否则使用可加。代码如下：


aust <- window(austourists,start=2005)
fit1 <- hw(aust,seasonal="additive")
fit2 <- hw(aust,seasonal="multiplicative")


其中fit1是可加性的，fit2是可乘性的。可以用accuracy函数来判断哪个效果更好：


rbind(accuracy(fit1),accuracy(fit2))

                      ME     RMSE      MAE           MPE     MAPE      MASE        ACF1
Training set 0.008115785 1.763305 1.374062 -0.2860248405 2.973922 0.4502579 -0.06272507
Training set 0.092062283 1.575631 1.254960 -0.0006505533 2.705390 0.4112302 -0.07955726


可以看到fit2的RMSE更小，因此可以认为可乘性的模型更佳。

Holt-Winters’平滑法一样可以考虑趋势衰减，是要设置damp参数为TRUE即可。事实上一种比较优秀的方法就是用Holt-Winters'季节可乘的趋势衰减平滑法，具体评论如下：“A method that often provides accurate and robust forecasts for seasonal data is the Holt-Winters method with a damped trend and multiplicative seasonality”（见7.3 Holt-Winters’ seasonal method | Forecasting: Principles and Practice）。实现如下：


hw(y, damped=TRUE, seasonal="multiplicative")  #y为要进行平滑的时间序列


自动化平滑模型选择与预测

介绍了很多方法，它们可能在不同的情况中适用。为了用一站式的服务，引入了ets函数。它不能进行预测，但是能够根据时间序列的状况，选择最佳的模型。简单操作：


aust <- window(austourists, start=2005)
fit <- ets(aust)
summary(fit)
#> ETS(M,A,M) 
#> 
#> Call:
#>  ets(y = aust) 
#> 
#>   Smoothing parameters:
#>     alpha = 0.1908 
#>     beta  = 0.0392 
#>     gamma = 2e-04 
#> 
#>   Initial states:
#>     l = 32.3679 
#>     b = 0.9281 
#>     s = 1.022 0.9628 0.7683 1.247
#> 
#>   sigma:  0.0383
#> 
#>   AIC  AICc   BIC 
#> 224.9 230.2 240.9 
#> 
#> Training set error measures:
#>                   ME  RMSE  MAE     MPE  MAPE   MASE
#> Training set 0.04837 1.671 1.25 -0.1846 2.693 0.4095
#>                ACF1
#> Training set 0.2006


相信上述的解释已经能够看出来什么意思，最上面模型选择最后是ETS(M,A,M) 表示（关于详细的模型展开式，参考7.5 Innovations state space models for exponential smoothing | Forecasting: Principles and Practice），残差用可乘、趋势用可加、季节用可乘的模型。如果要直接进行预测，可以使用forecast函数，它接收一个时间序列，自动使用ETS函数进行模型选择并预测。