用R语言做时间序列分析.docx
gCM用R语言做时间序列分析时间序列（time series ）是一系列有序的数据。通常是等时间间隔的采样数据。如果不是等间隔，则一般会标注每个数据点的时间刻度。下面以time series普遍使用的数据airline passenger为例。这是八一年的每月乘客数量，单位是千人次。ooco200220042006200820102012Time如果想尝试其他的数据集，可以访问这里:可以很明显的看出，airli ne passe nger的数据是很有规律的。time series data mining主要包括decompose（分析数据的各个成分，例如趋势，周期性）prediction（预测未来的值），classificatio n（对有序数据序列的feature提取与分类），clusteri ng（相似数列聚类）等。这篇文章主要讨论 prediction （ forecast，预测）问题。即已知历史的数据，如何准确预测未来的数据。先从简单的方法说起。给定一个时间序列，要预测下一个的值是多少，最简单的思路是什么呢？ exponential smoothing（1） mean （平均值）：未来值是历史值的平均。（指数衰减）：当去平均值得时候，每个历史点的权值可以不一样。最自然的就是越近的点赋予越大的权重。=aX± + C/X2 + a3Xs + .或者，更方便的写法，用变量头上加个尖角表示估计值X"= aX + (1 - )及（3） snaive ：假设已知数据的周期，那么就用前一个周期对应的时刻作为下一个周期对应时刻的预测值（4） drift :飘移，即用最后一个点的值加上数据的平均趋势Xt+h 11 =禺+占2心-斗-JL= Xt+占(良-如Tt = 介绍完最简单的算法，下面开始介绍两个time series 里面最火的两个强大的算法:Holt-Winters 和ARIMA。上面简答的算法都是这两个算法的某种特例。（5） Holt-Winters ： 三阶指数平滑Holt-Winters的思想是把数据分解成三个成分：平均水平（ level），趋势（trend ），周期性（seasonality ）。R里面一个简单的函数stl就可以把原始数据进行分解:022009EPpu上WElu一阶Holt Winters假设数据是stationary的（静态分布），即是普通的指数平滑。二阶算法假设数据有一个趋势，这个趋势可以是加性的（addtive,线性趋势），也可以是乘性的（multiplicative,非线性趋势），只是公式里面一个小小的不同而已。三阶算法在二阶的假设基础上，多了一个周期性的成分。同样这个周期性成分可以是additive和multiplicative的。举个例子，如果每个二月的人数都比往年增加1000人，这就是additive ;如果每个二月的人数都比往年增加120%，那么就是multiplicative。Holt-Winters Exponential smoothing数)tq=CtZf + ( Q)S" At) Double exponential smoothing51 = Ti龄= art + (1 -+ bi)52 = J 1 一工 0 /3(5f 5(_1) + (1 Ftm - St + mbtTriple exponential smoothing死=7。唯上 ft + (1 n)(的i + 如-)L-#(th. - 5(J) + (1 -如-iQ 7 h (1-7) Q-l斤 + m =临 + 77也)G1+l+gi) TTind £ JR里面有Holt-Winters的实现，现在就可以用它来试试效果了。我用前十年的数据去预测最后一年的数据。性能衡量采用的是RMSE。当然也可以采用别的metrics :r>t=LMAE 二广 1 工 Yt - aftlt=iMSE=人(y"M RMSE 二、t=lrpMARE = 100n-1 工 IM ftl/lytlt二i预测结果如下：Forecasts from Holt-Winters multiplieath/e method结果还是很不错的。(6) ARIMA : AutoRegressive Integrated Moving AverageARIMA是两个算法的结合：AR和MA。其公式如下：血二 C +AiQ-i-是白噪声，均值为0, C是常数。ARIMA的前半部分就是Autoregressive ：,后半部分是movingaverage :_° AR实际上就是一个无限脉冲响应滤波器(infinite impulse resopnse ) , MA 是 一个有限脉冲响应(finite impulse resopnse )，输入是白噪声°ARIMA 里面的 I 指 Integrated(差分)。ARIMA ( p,d,q )就表示p阶AR, d次差分，计特性（mean，varianeeq阶MA。为什么要进行差分呢？ARIMA的前提是数据是stationary 的，也就是说统correlation等)不会随着时间窗口的不同而变化。用数学表 示就是联合分布相同：+T?，1八泳)一网(赴is*i龙起开当然很多时候并不符合这个要求，例如这里的airli ne passe nger 数据。有很多方式对原始数据进行变换可以使之statio nary ：(1)差分，即Integrated。例如一阶差分是把原数列每一项减去前一项的值。二阶差分是一阶差分基础上再来一次差分。这是最推荐的做法(2 )先用某种函数大致拟合原始数据，再用ARIMA处理剩余量。例如，先用一条直线拟合airline passenger 的趋势，于是原始数据就变成了每个数据点离这条直线的偏移。再用ARIMA去拟合这些偏移量。(3)对原始数据取log或者开根号。这对varianee不是常数的很有效。如何看数据是不是stationary呢？这里就要用到两个很常用的量了 ：ACF (autoeorrelation function )禾口 PACF (patial auto eorrelation function)。又寸于non-stationary 的数据，ACF图不会趋向于0，或者趋向0的速度很慢。下面是三张ACF图，分别对应原始数据，一阶差分原始数据，去除周期性的一阶差分数据：acfLagu.<diff act<00051 015cpodiff act andemave seasonalityCM确保stationary之后，下面就要确定p和q的值了。定这两个值还是要看ACF和PACF ：ModelAR modek用PACF廛|去确认p的倩MF开乞犬才旨数寰减至0正员交替，衰减至01个或者多个倩大于6其余大致为0在某个点之后开始指数衰减至0所有点都罡0固定点上有大于。的值不衰减至0MA model, q等于知:F中大于0的个数AR和MA模型的结合原始数据完全随机原始数据有周期性-原始数据不罡占tmtiongryl确定好P和q之后，就可以调用R里面的arime函数了。值得一提的是，R里面有两个 很强大的函数：ets和auto.arima。用户什么都不需要做，这两个函数会自动挑选一个 最恰当的算法去分析数据。在R中各个算法的效果如下:meanf(train, h=12226.2657naive(tfaint h 二 12)snaiveftrain, h-12)102.976550.70832rwfftrairt, h-12, drift=T)92.66636sesftrainA ti=12Jnitial='simplef alpha=0.2)9.77035holtftrain, h=12/damped=F/ initiakAsiimpfen, beta=0.65)76.86677hw(train3 h=123seasonal-multiplicative1)16.36156ets(train)stlf(train|2439025222.07215auto.arimaftrain)23.538735arimpftrain, order = c(0? 15 3), seasonal=list(order=c(0JT3)5period=12)1835567代码如下:passenger = read.csv(pv-unlist(passenger)pt<-ts(p,frequency=Plot(pt)train<-window(pt,start='passenger.csv* ,header=F,sep= ")12 ,start= 2001 )2001 ,end= 2011 +11/12) test <-window(pt,start=2012)library(forecast)pred_meanf<-meanf(train,h=rmse( test ,pred_meanf$mean) #12)226 2657pred_naive<-naive(train,h=rmse(pred_naive$mean,12)test )#102.9765pred_snaive<-snaive(train,h=12)rmse(pred_snaive$mean,test )# 50. 70832pred_rwf<-rwf(train,h=12, drift=T)rmse(pred_rwf$mean, test )# 92.66636pred_ses <- ses(train,h=12 ,initial= 'simple' ,alpha= 0. 2)rmse(pred_ses$mean, test) # 89.77035
