时间序列（time series ）是一系列有序的数据。通常是等时间间隔的采样数据。如果不是 等间隔，则一般会标注每个数据点的时间刻度。 下面以time series 数量，单位是千人次。
普遍使用的数据 airline passenger 为例。 这是^一年的每月乘客
Time
如果想尝试其他的数据集，可以访问这
里： https://datamarket.com/data/list/?q=provider:tsdl 可以很明显的看出，airli ne passe nger time series data mining 期性），prediction 与分类），clusteri ng
的数据是很有规律的。
（分析数据的各个成分，例如趋势，周
（对有序数据序列的 feature 提取
主要包括 decompose
（预测未来的值）， classificatio n （相似数列聚类）等。
这篇文章主要讨论 prediction （ forecast，预测）问题。 测未来的数据。
即已知历史的数据，如何准确预
先从简单的方法说起。 给定一个时间序列，要预测下一个的值是多少， 呢？
最简单的思路是什么
（1） mean （平均值）：未来值是历史值的平均。
（2） exponential smoothing （指数衰减）：当去平均值得时候，每个历史点的权值可
以不一样。最自然的就是越近的点赋予越大的权重。
= aX± + ct^X2 + a3X3 + …
或者，更方便的写法，用变量头上加个尖角表示估计值
Xt+1 = aXt + （1 - a）Xt （3） sn aive :假设已知数据的周期，那么就用前一个周期对应的时刻作为下一个周期对 应时刻的预测值
（4） drift :飘移，即用最后一个点的值加上数据的平均趋势
t
Xt+h|t =禺+占2心-斗-丄= xt +占（罠-如
Tt = •
介绍完最简单的算法，下面开始介绍两个 time series
里面最火的两个强大的算法：
Holt-Winters
和ARIMA 。上面简答的算法都是这两个算法的某种特例。
（5） Holt-Winters : 三阶指数平滑
Holt-Winters
的思想是把数据分解成三个成分：平均水平（
level ），趋势（trend ），周期性（seasonality ）。R里面一个简单的函数 stl就可以把原始数据进行分解：
2002 2004 2006 2008 2010 2012
time
一阶Holt — Winters 假设数据是stationary
的（静态分布），即是普通的指数平滑。二阶
（additive,线性趋势），也可以是乘性的
三阶算法在二阶的假
算法假设数据有一个趋势，这个趋势可以是加性的
（multiplicative, 非线性趋势），只是公式里面一个小小的不同而已。 设基础上，多了一个周期性的成分。同样这个周期性成分可以是 的。举个例子，如果每个二月的人数都比往年增加
additive 和multiplicative
1000人，这就是additive ;如果每个
二月的人数都比往年增加 120%，那么就是multiplicative 。
Holt-Winters
• Exponential smoothing
数)—TQ
= CtZf + ( ] — Q)S[,扌 A t)
• Double exponential smoothing
s
i = Ti 龄=art + (1 -
— /3(5f — -5(_1) + (1 —
+ b—i)
5 = J 1 —工0
Ft^m = St + mbt
• Triple exponential smoothing
死=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 =广 工 \\Yt - ■ftl
1
t=i
预测结果如下:
n
■MRMSE =、 MSE = ^(yf - 2 t=l
rp
MAPE = 100n_1 工 IM 一 ftl/lytl
t=i
r>
t=L
Forecasts from Holt-Winters・ multiplieath/e method
结果还是很不错的。 (6)
ARIMA : AutoRegressive Integrated Moving Average
ARIMA是两个算法的结合：AR和MA。其公式如下：
•血=c +^iQ-i
•是白噪声，均值为 0, C是常数。 '_ '
—
■
'
_ Vs?
ARIMA 的前半部分就是 Autoregressive :
，后半部分是 moving
y _ j, t ,
average :
n .
。 AR实际上就是一个无限脉冲响应滤波器
_
(infinite impulse resopnse ) , MA 是一 -个有限脉冲响应( finite impulse resopnse )，
输入是白噪声。
ARIMA 里面的I指Integrated （差分）。 ARIMA （ p,d,q ）就表示p阶AR, d次差分，
ARIMA的前提是数据是stationary
的，也就是说统
q阶MA。 为什么要进行差分呢？
计特性（mean，varianee ，correlation 等）不会随着时间窗口的不同而变化。用数学表 示就是联合分布相同：
+T ? ■ * -1
^tk ） — Fx （赴 i s * i 龙起开
airli ne passe nger
数据。有很多方式对原
当然很多时候并不符合这个要求，例如这里的 始数据进行变换可以使之 statio nary : （1）差分，即Integrated
。例如一阶差分是把原数列每一项减去前一项的值。二阶差
分是一阶差分基础上再来一次差分。这是最推荐的做法 （2 ）先用某种函数大致拟合原始数据，再用 合airline passenger
ARIMA处理剩余量。例如，先用一条直线拟
的趋势，于是原始数据就变成了每个数据点离这条直线的偏移。再
用ARIMA去拟合这些偏移量。
（3）对原始数据取log或者开根号。这对 varianee 不是常数的很有效。 如何看数据是不是 stationary 呢？这里就要用到两个很常用的量了： eorrelation function non-stationary
）禾口 PACF（patial auto eorrelation function）
ACF （auto 。对于
下面是三张ACF
的数据，ACF图不会趋向于0 ,或者趋向0的速度很慢。
图，分别对应原始数据，一阶差分原始数据，去除周期性的一阶差分数据：
00
05
1 0
15
Lag
diff act
diff act and『emave seasonality
Lag
确保stationary
之后，下面就要确定 p和q的值了。定这两个值还是要看
ACF和PACF :
MF形状 Model
AR modek用PACF團去确认p的倩 MA model, q等于知:F中大于0的个数 AR和MA模型的结合
扌旨数寰减至0
正员交替，衰减至0
1个或者多个倩大于6其余大致为0 在某个点之后开始指数衰减至0
所有点都罡0 ________________________ 固定点上有大于。的值 不衰减至0
原始数据完全随机 ___________________ 原始数据有周期性 — ■■原始数据不罡 占tmtiongryl
确定好p和q之后，就可以调用 R里面的arime函数了。 值得一提的是，R里面有两个 很强大的函数：ets和auto.arima 。用户什么都不需要做，这两个函数会自动挑选一个 最恰当的算法去分析数据。 在R中各个算法的效果如下:
；p rrn
meanf(train, h=12} naive(tfainf h12)
二
226.2657 102.9765 50.70832 92.66636 89.77035 76.86677 16.36156 24390252 22.07215 23.538735 1835567
snaiveftrain, h-12) rwfftrairt, h-12, drift=T)
sesftrain^ ti=12,initial=lsimple,f alpha=0.2)
holtftrain, h=12/damped=F/ initiak^siimpfe\ hw(train, h=12,seasonal-multiplicative1) ets(train) stlf(train|
auto.arimaftrain)
arimpftrain, order = c(0? 1, 3), seasonal=list(order=c(0,lT3), period=12))
代码如下:
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 . 2657
pred_naive<-naive(train,h= rmse(pred_naive$mean,
12) test )# 102 . 9765
pred_snaive<-snaive(train,h= rmse(pred_snaive$mean,
12)
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.86677 without beta= 12 ,seasonal= 'multiplicative',
test )# 16. 36156
0. 65 it would be 84 . 41239
pred_hw<-hw(train,h= rmse(pred_hw$mean
fitv-ets(train)
accuracy(predict(fit.
12), test ) # 24. 390252
pred_stlf<-stlf(train) rmse(pred_stlf$mean,
test )# 22.07215
plot(stl(train,s.window= by Loess fitv-auto.arima(train) accuracy(forecast(fit,h=
\"periodic\" )) #Seasonal Decomposition of Time Series
12), test ) # 23.538735
ma = arima(train, order = c( period= 12)) pv-predict(ma, 12)
0, 1, 3), seasonal=list(order=c( 0,1,3),
accuracy(p$pred, test ) # 18. 55567 BT = Box. test (ma$residuals, lag= \"Ljung-Box\" , fitdf= 2)
30, type =
因篇幅问题不能全部显示，请点此查看更多更全内容