R语言与数据分析之七：时间序列简单指数平滑

上篇我们对时间序列数列有了整体的认识并将时间序列进行了分解，今天和小伙伴们分享常用预测算法中相对最简单的：简单指数平滑法。简单指数平滑适用于可用相加模型描述，并且处于恒定水平和没有季节变动的时间序列地短期预测。

简单指数平滑法提供了一种方法估计当前时间点上的水平。为了更加准确的估计当前时间的水平，我们使用alpha参数来控制平滑，alpha的取值在0-1之间。当alpha越接近0，临近预测的观测值在预测中的权重就越小。

我们采用伦敦1813年到1912年全部的每年每英尺降雨量来做分析对象，首先读入相关数据和绘制出序列图：

rain <- scan("http://robjhyndman.com/tsdldata/hurst/precip1.dat",skip=1)
rainseries <- ts(rain,start=c(1813))
plot.ts(rainseries)

从该图中可以看出整个曲线处于大致不变的水平，且随机变动在整个时间序列范围内也可以认为是大致不变的，所以该序列可以大致被描述为一个相加模型，因此我们可以使用简单指数平滑法进行预测。我们采用R中的HoltWinters()函数，为了能够使用HoltWinters中的指数平滑，我们需要进行参数设置：beta=FALSE和gamma=FALSE,预测结果如下图：

HoltWinters()告诉我们alpha参数的估计值约为0.024，非常接近0，说明该序列比较平稳，默认情况下HoltWInters仅会给出原始时间序列所覆盖时期内的预测，预测值存在房子名为为fitted的变量中，我们可以通过rainseriesforecasts$fitted来获取这些值。

此外我们可以画出预测值和实际值，看看预测效果：

plot(rainseriesforecasts)

从之前的alpha和上图，可见我们预测的过于平滑，R提供了样本预测误差平方和（SSE）来衡量预测效果。可以通过rainseriesforecasts$SSE来获取该值。另外我们还可以指定时间序列水平的初始值，常见的做法为采用时间序列第一个值,如这里采用1813年的23.56，具体实现如下：

HoltWinters(rainseries, beta=FALSE, gamma=FALSE, l.start=23.56)

为了预测未来时期的值，我们需要借助forecast包，假设我们需要对未来8年的下雨量进行预测，具体实现如下:

library("forecast")
rainseriesforecasts2 <-forecast.HoltWinters(rainseriesforecasts, h=8)

同时forecast.HoltWinters()函数给出了80%和95%的预测区间，为了更方便地查看预测结果，我们画出了预测效果图，具体实现和结果如下：

plot.forecast(rainseriesforecasts2)

其中蓝线是预测1913-1920间的降雨量，深灰色阴影区域为80%的预测区间，浅灰色阴影区域为95%的预测区间。

同样forecast也提供了”预测误差”的统计指标(residuals),借助”预测误差“我们可以评估预测是否有改进的可能性：如果预测误差是相关的，则很可能是简单指数平滑预测可以被另外一种预测技术优化。为了验证预测误差的相关性（act()），我们获取样本误差中1-20阶（期数）的相关图，具体R的实现和结果如下：

acf(rainseriesforecasts2$residuals, lag.max=20)#可以使用lag.max 指定我们想要看到的最大阶数

观看上图可以发现自相关系数在第3期的时候达到了置信界限。为了验证在滞后1-20阶时非0自相关属性是否显著，可以借助R语言的（Box.test()）Ljung-Box检验。具体实现如下：

统计量为17.4，并且P值是0.6（置信度只有40%）这样的值不足以拒绝“预测误差在1-20阶是非零自相关

证明预测误差在1-20阶是非零自相关的。

Ps: p值大，说明为纯随机序列。p值小，非纯随机序列，置信水平(1-p)

除了上述的验证方式，我们还可以验证预测误差是正太分布，并且均值为0，方差不变，为了做到这点，我们可以绘制出一个样本内预测误差图：

plot.ts(rainseriesforecasts2$residuals)

观察上图，我们可以归纳出整个时间区间内预测误差属于大致不变的正态分布，且均值接近0，为了更具体的展现，我们需要借助少量的代码，首先构建函数plotForecastErrors：

plotForecastErrors <- function(forecasterrors) 
{ 
	# make a red histogram of the forecast errors: 
	mysd <- sd(forecasterrors) 
	hist(forecasterrors, col="red", freq=FALSE) 
	# freq=FALSE ensures the area under the histogram = 1
	# generate normally distributed data with mean 0 and standard deviation mysd
	mynorm <- rnorm(10000, mean=0, sd=mysd) 
	myhist <- hist(mynorm, plot=FALSE) 
	# plot the normal curve as a blue line on top of the histogram of forecast errors:
	points(myhist$mids, myhist$density, type="l", col="blue", lwd=2)	
}

在控制台调用该函数：

source("plotForecastErrors.R")
plotForecastErrors(rainseriesforecasts2$residuals)

结果如图：

观察上图可见，预测的误差呈现均值为0 ，方差恒定的正态分布，可见预测算法无法改进了。