用R语言做时间序列分析.docx
《用R语言做时间序列分析.docx》由会员分享，可在线阅读，更多相关《用R语言做时间序列分析.docx（12页珍藏版）》请在冰点文库上搜索。
用R语言做时间序列分析
用R语言做时间序列分析
时间序列（timeseries）是一系列有序的数据。
通常是等时间间隔的采样数据。
如果不是等间隔，则一般会标注每个数据点的时间刻度。
下面以timeseries普遍使用的数据airlinepassenger为例。
这是^一年的每月乘客
数量，单位是千人次。
Time
如果想尝试其他的数据集，可以访问这
里：
可以很明显的看出，airlinepassenger的数据是很有规律的。
timeseriesdatamining主要包括decompose（分析数据的各个成分，例如趋势，周
期性），prediction（预测未来的值），classification（对有序数据序列的feature提取
与分类），clustering（相似数列聚类）等。
这篇文章主要讨论prediction（forecast，预测）问题。
即已知历史的数据，如何准确预
测未来的数据。
先从简单的方法说起。
给定一个时间序列，要预测下一个的值是多少，最简单的思路是什么
呢？
（1）mean（平均值）：
未来值是历史值的平均。
（2）exponentialsmoothing（指数衰减）：
当去平均值得时候，每个历史点的权值可
以不一样。
最自然的就是越近的点赋予越大的权重。
=aX±+ct^X2+a3X3+…
或者，更方便的写法，用变量头上加个尖角表示估计值
Xt+1=aXt+（1-a）Xt
（3）snaive:
假设已知数据的周期，那么就用前一个周期对应的时刻作为下一个周期对
应时刻的预测值
（4）drift:
飘移，即用最后一个点的值加上数据的平均趋势
t
Xt+h|t=禺+占2心-斗-丄=xt+占（罠-如
Tt=•
介绍完最简单的算法，下面开始介绍两个timeseries里面最火的两个强大的算法：
Holt-Winters和ARIMA。
上面简答的算法都是这两个算法的某种特例。
（5）Holt-Winters:
三阶指数平滑
Holt-Winters的思想是把数据分解成三个成分：
平均水平（level），趋势（trend），周
期性（seasonality）。
R里面一个简单的函数stl就可以把原始数据进行分解：
200220042006200820102012
time
一阶Holt—Winters假设数据是stationary的（静态分布），即是普通的指数平滑。
二阶
算法假设数据有一个趋势，这个趋势可以是加性的（additive,线性趋势），也可以是乘性的
（multiplicative,非线性趋势），只是公式里面一个小小的不同而已。
三阶算法在二阶的假
设基础上，多了一个周期性的成分。
同样这个周期性成分可以是additive和multiplicative
的。
举个例子，如果每个二月的人数都比往年增加1000人，这就是additive;如果每个
二月的人数都比往年增加120%，那么就是multiplicative。
Holt-Winters
•Exponentialsmoothing
数）—Tq
=CtZf+（]—Q）S[,扌At）
•Doubleexponentialsmoothing
si=Ti龄=art+（1-+b—i）
5=J1—工0—/3（5f—-5（_1）+（1—
Ft^m=St+mbt
•Tripleexponentialsmoothing
死=Tq
啣=ft——+（1—n）（的-i+如-[）
L
=#（軌一5（_i）+（1-①如-i
Q=7—+（1-7）Q-l
斤+m=临+77也）G-1+l+g-i）TTind£J
R里面有Holt-Winters的实现，现在就可以用它来试试效果了。
我用前十年的数据去预测最后一年的数据。
性能衡量采用的是RMSE。
当然也可以采用别的metrics:
MAE=广1工\Yt-
t=i
■ftl
n
MSE=^（yf-
t=l
■M2
RMSE=、
r>
t=L
rp
MAPE=100n_1工IM一
t=i
ftl/lytl
预测结果如下:
ForecastsfromHolt-Winters・multiplieath/emethod
结果还是很不错的。
（6）ARIMA:
AutoRegressiveIntegratedMovingAverage
ARIMA是两个算法的结合：
AR和MA。
其公式如下：
•血=c+^iQ-i
•是白噪声，均值为0,C是常数。
ARIMA的前半部分就是Autoregressive:
'_'—■'，后半部分是moving
y_j,t,_Vs?
n.
average:
_。
AR实际上就是一个无限脉冲响应滤波器
（infiniteimpulseresopnse）,MA是一-个有限脉冲响应（finiteimpulseresopnse），
输入是白噪声。
ARIMA里面的I指Integrated（差分）。
ARIMA（p,d,q）就表示p阶AR,d次差分，
q阶MA。
为什么要进行差分呢？
ARIMA的前提是数据是stationary的，也就是说统
计特性（mean，varianee，correlation等）不会随着时间窗口的不同而变化。
用数学表示就是联合分布相同：
+t?
■*-1^tk）—Fx（赴is*i龙起开
当然很多时候并不符合这个要求，例如这里的airlinepassenger数据。
有很多方式对原
始数据进行变换可以使之stationary:
（1）差分，即Integrated。
例如一阶差分是把原数列每一项减去前一项的值。
二阶差
分是一阶差分基础上再来一次差分。
这是最推荐的做法
（2）先用某种函数大致拟合原始数据，再用ARIMA处理剩余量。
例如，先用一条直线拟
合airlinepassenger的趋势，于是原始数据就变成了每个数据点离这条直线的偏移。
再
用ARIMA去拟合这些偏移量。
（3）对原始数据取log或者开根号。
这对varianee不是常数的很有效。
如何看数据是不是stationary呢？
这里就要用到两个很常用的量了：
ACF（auto
eorrelationfunction）禾口PACF（patialautoeorrelationfunction）。
对于
non-stationary的数据，ACF图不会趋向于0,或者趋向0的速度很慢。
下面是三张ACF
图，分别对应原始数据，一阶差分原始数据，去除周期性的一阶差分数据：
Lag
diffact
00
05
10
15
diffactand『emaveseasonality
Lag
确保stationary之后，下面就要确定p和q的值了。
定这两个值还是要看ACF和PACF:
MF形状
扌旨数寰减至0
正员交替，衰减至0
1个或者多个倩大于6其余大致为0
在某个点之后开始指数衰减至0
所有点都罡0
固定点上有大于。
的值
不衰减至0
Model
ARmodek用PACF團去确认p的倩
MAmodel,q等于知:
F中大于0的个数
AR和MA模型的结合
原始数据完全随机
原始数据有周期性—
■■原始数据不罡占tmtiongryl
确定好p和q之后，就可以调用R里面的arime函数了。
值得一提的是，R里面有两个很强大的函数：
ets和auto.arima。
用户什么都不需要做，这两个函数会自动挑选一个最恰当的算法去分析数据。
在R中各个算法的效果如下:
；prrn
meanf（train,h=12}
226.2657
naive（tfainfh二12）
102.9765
snaiveftrain,h-12）
50.70832
rwfftrairt,h-12,drift=T）
92.66636
sesftrain^ti=12,initial=lsimple,falpha=0.2）
89.77035
holtftrain,h=12/damped=F/initiak^siimpfe",beta=0.65）
76.86677
hw（train,h=12,seasonal-multiplicative1）
16.36156
ets（train）
24390252
stlf（train|
22.07215
auto.arimaftrain）
23.538735
arimpftrain,order=c（0?
1,3）,seasonal=list（order=c（0,lT3）,period=12））
1835567
代码如下:
passenger=read.csv（'passenger.csv',header=F,sep=''）
pv-unlist（passenger）
pt<-ts（p,frequency=12,start=2001）
plot（pt）
train<-window（pt,start=2001,end=2011+11/12）test<-window（pt,start=2012）
library（forecast）
pred_meanf<-meanf（train,h=12）
rmse（test,pred_meanf$mean）#226.2657
pred_naive<-naive（train,h=12）
rmse（pred_naive$mean,test）#102.9765
pred_snaive<-snaive（train,h=12）
rmse（pred_snaive$mean,test）#50.70832
pred_rwf<-rwf（train,h=12,drift=T）
rmse（pred_rwf$mean,test）#92.66636
pred_ses<-ses（train,h=12,initial='simple',alpha=0.2）
rmse（pred_ses$mean,test）#89.77035
rmse（pred_holt$mean.
test）#76.86677withoutbeta=
0.65itwouldbe
84.41239
pred_hw<-hw（train,h=
12,seasonal='multiplicative'
rmse（pred_hw$mean,
test）#16.36156
fitv-ets（train）
pred_stlf<-stlf（train）
rmse（pred_stlf$mean,test）#22.07215
plot（stl（train,s.window=
byLoessfitv-auto.arima（train）
accuracy（forecast（fit,h=12）,test）#23.538735
ma=arima（train,order=c（0,1,3）,seasonal=list（order=c（0,1,3）,
period=12））
pv-predict（ma,12）
accuracy（p$pred,test）#18.55567BT=Box.test（ma$residuals,lag=30,type=
"Ljung-Box",fitdf=2）
升级会员 