R语言2——时间序列分析

在统计研究中，常用按时间顺序排列的一组随机变量 来表示一个随机事件的时间序列，简记为 或 ，用 或 表示该随机序列的n个有序观察值。时间序列分析就是解释随机时序 的性质，它通过分析它的观察值序列 的性质，由观察值序列的性质来推断随机时序 的性质。

第一步，在R环境下将数据导入，生成时序对象。数据导入在上一节中已经有写过，在此不赘述。

tproduction<-ts(production,start=c(1962,1),frequency=12)
tproduction
plot(production)

ps：读取时间序列的部分数据采用代码：

tproduction.subset<-window(tproduction,start=c(1966,5),end=c(1972,6))

这样，时间序列已经生成。从图像上来看，序列总体呈上升趋势，序列属于非平稳序列。

第二步，对时序数据进行描述和可视化，通常采用对时序平滑化探究其总体趋势，并对其进行分解以观察时序中是否存在季节性因素。平滑处理是为了去除时间序列显著的随机或误差成分。存在季节性因素的时间序列数据可以通过季节性分解，分解为趋势因子、季节性因子和随机因子。其中趋势因子能捕捉到长期变化；季节性因子能捕捉到一年内的周期变化；随机因子能捕捉到那些不能被趋势或季节效应解释的变化。

时序数据中通常有很显著的随机或误差成分。为了辨明数据中的规律，撇开这些波动，最简单的办法就是进行简单移动平均。如果每个数据点都可用这一点和其前后两个点的平均值来表示，这就是居中移动平均。R中有几个函数可以做简单移动平均，包括TTR包中的SMA（）函数，zoo包中的rollmean（）函数，forecast包中的ma（）函数，还有R中自带的ma（）函数，都可以进行平滑处理。

library(forecast)
> plot(ma(tproduction,3))
> plot(ma(tproduction,7))
> plot(ma(tproduction,15))
> plot(ma(tproduction,25))

随着平均的观测值的个数的增大，曲线逐渐趋于平滑。25个观测值的平滑效果最好，还可以明显看出数据的上升趋势，在1972年到1974年间，牛奶销售表现出持平的效果，有微小的上升趋势。对于间隔大于1天的时序数据，数据可能会存在季节性因子，我们需要的不只是总体趋势。 季节性分析可以通过相加模型或相乘模型来分解数据。即三个因子相加或相乘。分解的常用方法是用LOESS光滑做季节性分解，可以通过R中的stl（）函数来实现：

stl(ts,s.window=,t.window=)

ts是将要分解的时序，参数s.window控制季节效应变化速度，t.window控制趋势项变化的速度，值越小变化越快。s.window=“periodic”可使得季节效应在各年间都一样。其中ts和s.window的参数是必须要设置的。stl函数只处理相加模型，将相乘模型可以进行对数转换变成相加模型。

注意：如果序列的波动于趋势成正比，那么相乘模型更适合这类情况。对于本实验的情况，序列的波动不与趋势成正比，因此采用相加模型即可。 选择好模型以后（相加模型直接进行下一步，相乘模型进行对数变换），开始分解时间序列。

fit<-stl(tproduction,s.window="period")
plot(fit)
fit$time.series

最后生成季节性分解图以及每个观测值各分解项的值（季节效应，趋势，随机波动）。通过对季节性分解图的观察发现，序列的趋势为单调增长，季节效应表明，牛奶在3月到7月的销量更高。由于每个图的y轴尺度不同，因此我们通过图中右侧的灰色长条来指示量级，即每个长条代表的量级一样。



第三步，通过对时序数据的描述和可视化，针对其性质对时序选择相应的方法进行预测。指数模型是用来预测时序未来值的最常用模型，他们的短期预测能力较好。单指数模型拟合的是只有常数水平项和时间点i处随即向的时间序列，它认为时间序列不存在趋势项和季节效应；双指数模型也叫Holt指数平滑，它拟合的是有水平项和趋势项的时序；三指数模型也叫Holt-Winters指数平滑，拟合的是有水平项、趋势项以及季节效应的时序。根据前面的数据描述的和可视化结果，该数据应采用方法三指数模型。R自带的HoltWinters()函数或者forecast包中的ets()函数可以拟合指数模型。

ets(ts,model="zzz")

其中，ts是要分析的时序，限定的模型的字母有三个，第一个为误差项（水平项），第二个为趋势项（斜率），第三个为季节项。可选的字母有A（相加模型），M（相乘模型），无（N），自动选择（Z）。

ets(ts,model="ANN") 单指数模型
ses(ts)
ets(ts,model="AAN") 双指数模型
holt(ts)
ets(ts,model="AAA") 三指数模型
hw(ts)

建立拟合模型的代码：

fit<-ets(tproduction,model="AAA")
fit

结果：

ETS(A,A,A) Call: ets(y = tproduction, model = "AAA") Smoothing parameters: alpha = 0.7222 beta = 1e-04 gamma = 2e-04 Initial states: l = 605.8538 b = 1.6751 s=-42.8457 -79.2868 -48.5788 -51.9412 -9.4986 32.5529 81.9608 110.5424 49.5335 34.5776 -58.3818 -18.6344 sigma: 7.1515 AIC AICc BIC 1555.844 1559.924 1608.951

拟合精度：

accuracy(fit)

精度结果：

ME RMSE MAE MPE MAPE MASE
Training set -0.03164542 7.151456 5.380847 -0.001196899 0.7350337 0.2527588
ACF1 Training set -0.004837951

未来值的预测：

pred<-forecast(fit,5)
pred

预测结果:

Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 1976 866.6590 857.4940 875.8239 852.6424 880.6755
Feb 1976 828.5864 817.2805 839.8923 811.2955 845.8773
Mar 1976 923.2189 910.1169 936.3210 903.1811 943.2568
Apr 1976 939.8496 925.1691 954.5302 917.3977 962.3016
May 1976 1002.5333 986.4278 1018.6388 977.9020 1027.1645

生成折线图：

plot(pred)

折线图结果：



蓝色线就是预测的曲线值

https://github.com/kabacoff/RiA2/blob/master/Ch15%20Time%20series.R

https://www.manning.com/books/r-in-action