第10章时间序列分析

第10章时间序列分析7/14/2019常用时间序列模型ARIMA模型时间序列对象在R软件中，使用时间序列建模前需要先将数据存储到一个时间序列对象中。我们可以使用函数ts()将数值类型的观测对象存储为时间序列对象。使用格式：ts(data

=

NA,

start

=

1,

end

=

numeric(),

frequency

=

1,…

)其中，data是时间序列观测值对象，必须为数值类型的向量、矩阵或数据框；

start是用来指定时间序列观测值对象的第一个时间点，比如2000年1月，则设置

start=c(2000,1)；end用来指定时间序列的终止时间点；frequency用来指定数据在一年中的频数。还可以通过函数as.ts()可以将对象转换成时间序列；通过函数is.ts()可以判断象是否为时间序列对象。ARIMA模型绘制时间序列图R软件中，可以使用plot.ts()函数来画出时间序列的时序图。plot.ts()用法同plot。根据平稳时间序列的均值和方差都为常数的性质，平稳序列的时序图显示该序列值始终在一个常数附近随机波动，而且波动的范围有界；如果有明显的趋势性或者周期性那它通常不是平稳序列。ARIMA模型时间序列检验分析（1）自相关性检验自相关图中的两条虚线标示置信界限是自相关系数的上下界。如果自相关系迅速衰减落入置信区间内，就可能是白噪声；如果自相关系数超出置信区间，那么表示存在相关关系，而且从哪一阶落在置信区间内，就表示自相关的阶数是几阶。使用格式：acf(x,

lag.max

=

NULL,type

=

c(correlation,

covariance,

partial),plot

=

TRUE,

na.action

=

na.fail,

demean

=

TRUE,

…)pacf(x,

lag.max,

plot,

na.action,

…)acf()函数中，参数x为观测值序列，acf()为观测值序列自相关函数，lag.max为与acf对应

的最大延迟，type为计算acf的形式，默认为correlation。当没有输出，即为acf(Series)时，画观测值序列的自相关系数图。pacf()函数中的输入参数与输出参数的含义同acf()函数的类似。在acf()和pacf()中设定

plot=FALSE可以得到自相关和偏自相关的真实值。ARIMA模型原始序列时序图 原始序列的自相关图时序图显示该序列具有明显的单调递增趋势，可以判断为是非平稳序列；自相关图显示自相关系数长期大于零，说明序列间具有很强的长期相关性，可以判断为非平稳序列；ARIMA模型时间序列检验分析（2）单位根检验对时间序列的平稳性检验通常使用单位根检验的方法。在R软件中，单位根检验使用fUnitRoots包中的unitrootTest()函数可以实现使用格式：unitrootTest(x,lags=1,type=c("nc","c","ct"),title=NULL,description=其中，输入参数x为观测值序列，lags为用于校正误差项的最大滞后项，type为单位根的回归类型，返回的参数p值，p值小于0.05表示满足单位根检验。单位根检验统计量对应的p值显著大于0.05，判断该序列为非平稳序列（非平稳序列一定不是白噪声序列）。ARIMA模型ARIMA建模分析（1）非平稳时间序列差分对于非平稳时间序列，首先需要对其进行差分直到得到一个平稳时间序列。在R软件中，可以使用diff()函数对时间序列进行差分运算。diff()函数的使用格式：diff(x,lag=1,differences=1,...)其中，输入参数"x"代表观测值序列；"lag"代表差分运算的步数，缺省值代表一步差分；"differences"代表差分运算的阶数，缺省值代表一阶差分。对一阶差分后的序列再次做平稳性判断过程同上。一阶差分之后序列的时序图一阶差分之后序列的自相关图结果显示，一阶差分之后的序列的时序图在均值附近比较平稳的波动、自相关图有很强的短期相关性、单位根检验p值小于0.05，所以一阶差分之后的序列是平稳序列。ARIMA模型（2）时间序列模型识别定阶使用R软件中的acf()和pacf()函数来分别给出时间序列的自相关图和偏自相关图。可根据自相关图和偏自相关图对时间序列模型进行定阶。根据自相关图和偏自相关图对时间序列模型进行定阶：

若平稳序列的偏相关系数是截尾的，而自相关系数是拖尾的，则序列适合AR模型；

若平稳序列的偏相关系数是拖尾的，而自相关系数是截尾的，则序列适合MA模型；若平稳序列的偏相关系数与自相关系数都是拖尾的，则序列适合ARMA

模型；ARIMA模型pacf(difsales,lag.max=30)一阶差分之后序列的自相关图 一阶差分后序列的偏自相关图一阶差分后自相关图中，ACF值在一阶后迅速跌入置信区间，并且数值徘徊在置信区间，没有收敛趋势，显示出拖尾性。偏自相关图中，PACF值在一阶后迅速跌入置信区间，并且有向零收敛的趋势，显示出截尾性，所以可以考虑用AR模型拟合1阶差分后的序列，即对原始序列建立ARIMA(1,1,0)模型。acf(difsales,lag.max=30)ARIMA模型另外，模型还可以通过BIC进行定阶。此处计算ARMA(p,q)当p和q均小于等于5的所有组合的BIC信息量，取其中BIC信息量达到最小的模型阶数。一阶差分后序列的BIC图显示BIC值从下往上，依次递减。模型选用变量的单元格用阴影表示。较好的模型（具有较低的BIC值）处于较高的行中。第一行中，test-lag1被选入模型，error-lag1到error-lag5均未被选取，取零阶。即在test-lag1，error-lag0处，取得最小BIC值。因此p值为1、q值为0时，BIC值最小。p、q定阶完成,，即对原始序列建立ARIMA(1,1,0)模型。ARIMA模型（3）ARIMA模型R软件中，arima()函数设置时序模式的建模参数，创建ARIMA时序模型或者把一个回归时序模型转换为ARIMA模型。使用格式：Arima(Series,

order,

seasonal,

period,method,...)其中，Series为观测值序列，order为构建的ARIMA(p,d,q)模型的参数，seasonal为模型的季节性参数，period为观测值序列的周期，method为估计模型参数所使用的方法。ARIMA模型（4）白噪声检验在使用ARIMA模型定阶完成后，还要对模型进行假设检验，检验残差序列是否为白噪声序列。R软件中，Box.test()函数用于检测序列是否符合白噪声检验。使用格式：Box.test(x,lag=1,type=c("Box-Pierce","Ljung-Box"),fitdf=0)其中，x为Arima()函数返回的结果对象的residuals残差。ARIMA模型p-value>0.05，说明残差为白噪声序列，通过白噪声检验。ARIMA模型预测（1）模型预测根据参数估计得到的时间序列模型对未来的序列值进行预测，通过forecast包中的函数forecast.Arima()来完成。使用格式：forecast.Arima(object,

h=10,

level=c(80,95),…)object是函数arima()返回的对象；h指定预测的时间点；level指定预测区间的置信水平，默认情况为80%和95%置信水平下的预测区间。Forecast的输入结果是指：通过arima模型，置信水平为80%和95%时，预测出未来五天的营业额及营业额区间。ARIMA模型预测（2）绘制原始值与预测值图形为查看原始值与预测值的图形效果，可以使用plot.forecast()函数绘制原始值与预测值的图形。使用格式：plot.forecast(x,…)其中，x为forecast.Arima()函数返回的结果对象。plot.forecast(forecast)使用函数plot.forecast()可以查看原始值和预测值的图形效果，例如从图中可以直观地看到营业额的发展趋势。其他时间序列模型（1）组合模型时间序列的变化主要受到长期趋势、季节变动、周期变动和噪声变动这四个因素的影响。根据序列的特点，可以构建加法模型和乘法模型。decompose()函数、stl()函数可以估计出时间序列中趋势的、季节性的和不规则的部分，而此时间序列需是可以用相加模型描述的。函数使用格式：decompose(x,

type

=

c("additive",

"multiplicative"),

filter

=

NULL)在decompose()函数中，x时间序列对象；type指定分解为加法模型还是乘法模型；filter是滤波系数。stl(x,

s.window,

s.degree=0,…)在stl()函数中，x同样为时间序列对；s.window因为没有默认值，所以必须手动设置，可以采用’periodic’或Loess方法提取季节跨度，若采用Loess方法，

s.window的值必须为大于7的奇数；s.degree可取1或0，为局部多项式拟合季节性提取的程度。sales<-ts(Data,start=c(2013,1),frequency=12)

##示例数据sales.de

=

decompose(sales,type="additive")sales.de$xJan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec2013

3023

3039

3056

3138

3188

3224

3226

3029

2859

2870

2910

30122014

3142

3252

3342

3365

3339

3345

3421

3443

3428

3554

3615

36462015

3614

3574

3635

3738

3707

3827

4039

4210

4493

4560

4637

47552016

4817$seasonalJan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec2013

40

33

50

44

-56

-65

79

22

-94

-50

-23

202014

40

33

50

44

-56

-65

79

22

-94

-50

-23

202015

40

33

50

44

-56

-65

79

22

-94

-50

-23

202016

40$trendJan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec2013

NA

NA

NA

NA

NA

NA

3053

3067

3087

3109

3125

31362014

3149

3174

3215

3268

3325

3381

3427

3460

3486

3514

3545

35802015

3626

3684

3760

3846

3931

4020

4116

NA

NA

NA

NA

NA2016

NA$randomJan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Dec2013

NA

NA

NA

NA

NA

NA

94

-59

-134

-189

-192

-1442014

-47

44

77

54

69

29

-86

-39

36

90

93

462015

-52

-143

-175

-152

-168

-127

-156

NA

NA

NA

NA

NA2016

NA$figure[1]

40

33

50

44

-56

-65

79

22

-94

-50

-23

20$type[1]

"additive"attr(,"class")[1]

"decomposed.ts"plot(sales.de)图中包含了四部分，自上而下依次为：原始时间序列的观测值、时间序列分解趋势图、时间序列分解季节变动图、时间序列分解噪声图。其他时间序列模型sales.stl

=

stl(sales,s.window

=

"periodic")sales.stlCall:stl(x

=

sales,

s.window

=

"periodic")Componentsseasonal

trend

remainderJan

20135.73

3109-91.9Feb

2013-53.40

3103-10.3Mar

2013-28.08

3096-12.2Apr

20130.99

309145.7May

2013-41.61

3086143.4Jun

2013-26.07

3085165.6…… ……

……

……Mar

2015

-28.08

3752Apr

2015 0.99

3843-89.2-105.6May

2015-41.613933-184.3Jun

2015-26.074036-182.5Jul

201532.134138-131.4Aug

2015-4.784243-28.0Sep

2015-7.684347153.5Oct

201519.36445685.1Nov

201537.73456435.4Dec

201565.68467513.8Jan

20165.73

478724.2plot(sales.stl)stl季节性时间序列分解图图中包含了四部分，自上而下依次为：原始时间序列的观测值、时间序列分解趋势图、时间序列分解季节变动图、残差自相关图。采用decompose()与stl()对时间序列进行分解后，建模、预测的工作一般通过指数平滑算法实现。其他时间序列模型（2）指数平滑法在R软件中，简单指数平滑法、Holt双参数线性指数平滑法和Winters线性和季节性指数平滑法可以通过HoltWinters()函数实现。HoltWinters()函数会给出时间序列自身时间段内的预测，生成的预测值存在一个变量列表fitted中。可以用plot()函数绘制原始数据和预测值的图形查看预测效果。使用格式：HoltWinters(x,

alpha

=

NULL,

beta

=

NULL,

gamma

=

NULL,seasonal

=

c("additive",

"multiplicative"),start.periods

=

2,

l.start

=

NULL,

b.start

=

NULL,s.start

=

NULL,optim.start

=

c(alpha

=

0.3,

beta

=

0.1,

gamma

=

0.1),optim.control

=

list())其中，x为时间序列对象；alpha、beta、gama都是HoltWinters算法的过滤系数。beta和gama是Holt指数平滑法或Winters指数平滑法的参数。如果beta设置为FALSE，该函数将做指数平滑；如果gama指定为

FALSE，那么该函数适用于拟合非季节性模型；如果gama与beta同时指定为FALSE，那么该函数适用于简单指数平滑模型；通过seasonal参数选择”additive"或"multiplicative"季节性模型，默认选择”additive”而仅在gamma非零时生效；start.period是用于x对象的frequency自动检测，不能小于2；l.start、b.start、s.start分别表示启动值、趋势值和季节分量的初始值；optim.start设置向量命名的组件alpha、beta、gamma包含优化的初始值，必须指定唯一需要的值，在只使用alpha、beta、gamma中的一个参数时忽略本参数。HoltWinters()函数仅得到预测模型，如果要对未来的时间做预测，需要调用forecast()函数。其他时间序列模型eta

=

TRUE,

gamma

=

TRUE)s1110.763889s2159.805556s3130.597222s446.722222s5-3.069444s6206.388889s7-4.444444s8

-195.236111s9

-205.611111s10

-181.361111s11

-90.694444s12

26.138889plot(hw.sales)该餐馆的营业数据存在明显的趋势性与季节性。因此可以采用Winters线性和季节性指数平滑法进行建模预测。hw.sales

=

HoltWinters(sales,alpha

=

TRUE,beta

=

TRUE,

gamma

=TRUE)hw.salesHolt-Winters

exponential

smoothing

with

trend

and

additive

seasonalcomponent.Call:HoltWinters(x

=

sales,

alpha

=

TRUE,

bSmoothing

parameters:alpha:

TRUEbeta

:

TRUEgamma:

TRUECoefficients:[,1]a

4790.861111b

-54.833333餐饮营业额Winters模型拟合图其他时间序列模型建模并对模型的残差进行自相关检验与白噪声检验。>library(forecast