R语言建模和预测

R 语言下时间序列建模和预测研究

摘要：

时间序列是指将某种现象某一个统计指标在不同时间上的各个数值，按时间先后顺序排列而形成的序列。研究时间序列的领域多的难以罗列，但是时间序列分析的目的一般有两个方面：一是认识产生观测序列的随机机制，即建立数据生成模型；二是基于序列的已知数据，对序列未来的可能取值给出预测或预报。本文在分析和理解平稳时间序列和非平稳时间序列的基础上，在R 语言环境下，根据某一地震观测数据的时间序列，实现了对该时间序列的平稳性分析，模型识别，参数估计，模型诊断，预测和更新。

关键字：

时间序列，平稳性分析，模型识别，参数估计，预测更新

1.引言

时间序列是指将某种现象某一个统计指标在不同时间上的各个数值，按时间先后顺序排列而形成的序列。随着现代社会的发展，数据的产量越来越大，对于数据时间序列的分析研究越来越多，时间序列的研究应用于很多领域，经济上用于统计和预测经济增长等各项指标，在天气预报上，也通过对云层的观测序列，预测未来的天气情况，在测量上的应用更是比比皆是。时间序列模型，可以根据其平稳性可以分为平稳时间序列模型和非平稳时间序列模型。而平稳时间序列模型又可以分为一般线性模型、滑动平均模型、自回归模型和自回归滑动混合模型;而非平稳模型则是多种多样。

2.模型概述

2.1 平稳时间序列模型

平稳时间序列满足期望为零，均值在所有时间上均为常数，且任意两个时刻的相关函数与时间t 无关，仅与两个时刻的时间差相关，平稳时间序列模型包括一般线性过程、滑动平均过程、自回归过程、自回归滑动平均过程[1]

。 2.1.1 一般线性过程

一般线性过程*Y t +可以表示成现在和过去白噪声变量的加权线性组合： Y t =e t +φ1e t?1+φ2e t?2+??? (2.1) 如果表达式的右边事实上是一个无穷级数，那么需要给权数φ加上一定条件，这样右边的表达式才在数学上有意义，为此，只需假设

∑φ

i 2

<∞∞i=1 (2.2)

2.1.2 滑动平均模型

当有限个系数φ不为零时，就得到了所谓的滑动平均过程，改变正负号，如下式所示：

Y t=e t?θ1e t?1?θ2e t?2?????θq e t?q(2.3)

(2.3)式为q阶的滑动平均过程，记作MA(q)。

2.1.3 自回归模型

自回归过程就是用自身做回归变量，具体来说，p阶自回归过程*Y t+满足方程：

Y t=?1Y t?1+?2Y t?2+???+?p Y t?p+e t(2.4) (2.4)式为p阶的自回归模型，记为AR(p),，序列Y t的当期值是自身最近p 滞后项和新信息项e t的组合，其中包含了序列在t 期无法用过去值来解释的所有新信息，因此对于每一个t ，假设e t独立于Y t?1，Y t?2，Y t?3???

2.1.4 自回归滑动平均混合模型

如果假定序列中部分是自回归过程，部分是滑动平均过程，就可以得到一个相当于普遍的时间序列模型：

Y t=?1Y t?1+?2Y t?2+???+?p Y t?p+e t?θ1e t?1?θ2e t?2?????θq e t?q(2.5)

(2.5) 式是Y t的自回归滑动平均模型，滑动平均部分的阶数为q,自回归部分的阶数为p,记为ARMA(p,q)。

2.2 非平稳时间序列模型

具有时变均值的时间序列都是非平稳的，形如：

Y t=μt+X t(2.6)

其中：μt为时变均值函数，X t为平稳序列。

该序列因为均值μt是不确定是不断变化的，因此该模型的时间序列趋势就是随机的，可能一直向上，一直向下，或是断断续续的上下颠簸，这样的模型不具备平稳性，平稳性模型是不适合的，必须有一个含有随机趋势的非平稳模型才是适当的，这就提出了ARIMA 模型。

如果一个时间序列*Y t+的d次差分W?d=Y t, Y t是一个ARMA的平稳模型，那么我们称*Y t+为自回归滑动平均求和模型，记为ARIMA(p,d,q),其中d为差分次数。

可以看出，ARIMA(p,d,q)，可以表示平稳模型和非平稳模型，比如说：

ARIMA(p,0,q)模型就是自回归滑动平均模型ARMA(p,q)，ARIMA(p,0,0)就是自回归模型AR(p)，ARIMA(0,0,q)就是滑动平均模型MA(q)。

3.模型识别

3.1 基本知识

在模型识别的过程中，具有两个非常重要的参数，分别为自相关函数ρk=γk/γ0和偏自相关函数?kk，如果某一时间序列的自相关函数随着滞后k 的增加而很快地下降为0，那么我们就认为该序列为平稳序列；如果自相关函数不随着k 的增加而迅速下降为0，就表明该序列不平稳。如果一个时间序列的自相关和偏相关图没有任何模式，而且数值很小, 那么该序列可能就是一些互相独立的无关的随机变量。其中自相关函数与偏相关函数又具有两个非常重要的性质，分别为拖尾性和截尾性，其中拖尾性是指随着滞后k 的增加，函数渐变为0，图像像是拖了一条尾巴；截尾性是指当滞后k＞p或者k＞q处，函数突然下降为 0 图像像是截断了尾巴一样。函数的拖尾性和截尾性对时间序列模型的识别具有重要作用，见表3-1：

图表3-1 ARMA模型ACF和PACF的一般特征

AR(p) MA(q) ARMA(p,q),p>0,q>0 ACF 拖尾滞后q阶后截尾拖尾

PACF 滞后p阶后截尾拖尾拖尾

3.2 模型的平稳性分析

3.2.1 原始数据的平稳性分析

本文中运用的原始数据是某一地震的Z观测数据，总数据有2400个历元，运用前2300个数据进行建模和预测更新，图表3-2 就是该数据前2300个历元的时间序列图。

图表3-2 原始数据的时间序列

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资\\1-2300.txt",header=F)

prop=ts(mydata)

plot(prop,xlab="time",ylab="Observation: Z")

图表3-3 原始数据的时间序列的ACF和PACF

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资\\1-2300.txt",header=F)

prop=ts(mydata)

acf(prop)

pacf(prop)

从图表3-2可以看出，该时间序列的是在0值上下小范围波动，到1300历元时波动范围变大，从10的三次方慢慢跃升到10的六次方量级，但是总体上时间序列没有表现出整体上升或是下降的趋势，始终在0值上下震荡，很有可能是平稳模型，但是从该时间序列的ACF和PACF可以看出，ACF和PACF并不满足平稳模型的拖尾或是截尾的特点，又有可能是非平稳模型。假设是非平稳模型，我们可以对原始的时间序列进行差分处理后画出其ACF和PACF进行分析，如图表3-4所示，画出的是10次差分后ACF和PACF。从图表3-4可以看出，在10次差分后还是没有满足拖尾或是截尾的特征，进而做了更多此时的差分，但是还是没有满足。这说明整个序列也并不满足ARIMA模型的条件。从数据的时间序列上看，

因为前半段的数据很平稳，突然发生地震，数据做了几个数量级的跨越并维持了很长时间，这样就会造成数据自相关函数和偏自相关函数的不正常，不能正常的表现出拖尾和截尾的特征，很难建立一个已有的模型和之对应，因此得出结论，原始数据很难整体建模。因此，本文就对原始数据进行了划分，对每一段分别建模。

图表3-4 10次差分后的ACF和PACF

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资\\1-2300.txt",header=F)

prop=ts(mydata)

prop1=diff(prop,difference=10)

acf(prop1)

pacf(prop1)

3.2.2 第一段数据的平稳性分析

从原始数据的时间序列可以看出，前1200个点的序列表现出了相对平稳的性质，在0值上下做小范围的波动，在第1208个历元的时候，数据的绝对值出现了跨数量级的变化，可以把第1208历元作为一个节点，为了探讨模型的预测和更新，因此第一段数据就是原始数据的前1107个历元。图表3-5画出了第一段数据的时间序列，图表给出了第一段时间序列的ACF和PACF。

图表3-6 第一段的时间序列

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107.txt",header=F) prop=ts(mydata)

plot(prop,xlab="time",ylab="Observation: Z")

图表 3-7 第一段时间序列的ACF和PACF

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107.txt",header=F) prop=ts(mydata)

acf(prop)

pacf(prop)

从图表3-6 可以看出，第一段的时间序列整体表现出了均值稳定的性质，都是在0值上下波动，没有明显的上升和下降趋势。从图表3-7，可以看出PACF 表现出了拖尾的特征，而ACF并没有。因为这是真实的观测数据，跟理论上的模型数据有一点的差异。根据第一段时间序列的均值稳定性和PACF的拖尾特征，暂且识别为ARMA模型，在后面的模型诊断中去检验该假设的正确性。

4.模型参数的估计

虽然第3章中第一段的模型暂定为ARMA模型，但是还不知道其ARMA模型的相应模型参数，本章将简述模型参数估计的方法。

4.1 模型参数估计的方法

4.1.1参数估计思路

在ARMA模型的参数估计中，由于参数p，q的取值一般不大，参数p，q 的获取方法是逐阶尝试，(p,q)的取值从低阶(1,1)、(1,2)、 (2,1)、(2,2)…向高阶尝试，在每一次尝试中定出模型，然后经过对估计模型的残差分析确定该模型是否被接受，即是否与原始的序列符合的很好；如若符合不好，调整(p,q)的值，继续尝试，重新进行残差分析，直至出现符合要求的模型参数为止[3]。

4.1.2模型参数的检验

如果模型被正确识别，且估计的模型参数合适，那么该模型应该与真实时间序列拟合很好，也就是说残差就应该是近似于具有白噪声性质，即残差均匀出现在0值附近，并且没有上下趋势，即表现为独立同分布的具有零均值和相同标准差的正态变量。残差就是真实值与模型值得差值。残差分析法，检验模型估计参数的正确性的步骤如下：

?观察残差的时间序列图，如果模型识别正确，且模型参数合适，那么残差的时间序列应该表现为在0值附近上下分布均匀且无任何趋势的长方形散

点图。

?检验残差的自相关性，如果模型识别正确，且模型参数合适，那么残差应该是相互独立的，具有良好的独立性，采用Ljung-Box检验统计量的值，若，p＞5%则认为残差有很好的独立性。

?检验残差的正态性，主要利用分位数-分位数图来检验残差的正态性，在R语言中，主要是用残差的分位数-分位数图，即qqnorm和qqline函数图。

4.2 模型参数估计的实现

第一段数据时间序列在暂定为ARMA模型后，就开始从低阶到高阶来估计模型参数，并不断进行残差检验，最终确定出了模型参数为p=5，q=2，残差检验的结果见图表4-1残差检验结果和图表4-2 残差正态性检验qq图。

图表4-1 残差检验结果

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,2),method="ML")

tsdiag(model)

r= model $residuals

Box.test(model,type="Ljung-Box",lag=6,fitdf=1)

由图表4-1 残差检验结果可以看出，首先从残差的时间序列图中可以看出：残差是均匀分布在0值上下的，没有任何趋势的散点图，满足残差检验的条件一。从残差的自相关ACF可以看出，除了第13阶处有一点超限，其他的都在限差之内，在Ljung-Box检验中，p=0.5348＞5%，因此可以说残差有很好的独立性，满足残差检验的第二个条件。

下面我们将在R预言下利用qqnorm和qqline函数画出残差的分位图，进而检验残差的正态性，检验结果如图表4-2所示。

图表4-2 残差模型的正态QQ图

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107.txt",header=F)

prop=ts(mydata)

model =arima(prop,order=c(5,0,2),method="ML")

qqnorm(model $residuals)

qqline(model $residuals)

从图表4-2中，可以看出QQ图中的点的分布非常接近一条直线，只有在两段的时候有一点偏差，大部分数据分布很好。该QQ图说明残差有很好的正态性。满足残差检验的第三个条件，因此我们接受模型ARMA(5,2)。

4.3 模型参数估计结果

图表4-3 模型参数估计结果

AR1 AR2 AR3 AR4 AR5 MA1 MA2 Intercept Value 1.2186 -1.2379 0.8243 -0.2611 0.2719 -1.1352 0.8301 -345.2400 s.e. 0.0365 0.0507 0.0552 0.0488 0.0336 0.0239 0.0225 33.4983 mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107.txt",header=F)

prop=ts(mydata)

model =arima(prop,order=c(5,0,2),method="ML")

model

5.预测与更新

5.1预测

时间序列建模的目的有两个：一是认识产生观测序列的随机机制；二是基于序列的已知数据，对序列未来的可能取值给出预测或预报。模型的建立我们已经完成，下面我们将根据模型和已知的数据，对未来数据进行预测。由于前1207个历元中我们预留了100个历元，因此我们在此预测未来100个历元的数据，预测结果如图表5-1所示。

图表5-1 预测未来100步

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,2),method="ML")

prop.fore=predict(model,n.ahead=100)

Up=prop.fore$pred+1.96*prop.fore$se

Low=prop.fore$pred-1.96*prop.fore$se

ts.plot(prop,prop.fore$pred,col=1:2,xlab="time",ylab="Observation: Z")

lines(Up,col="green",lty="dashed")

lines(Low,col="green",lty="dashed")

由图表5-1是对该1107个时间历元序列进行的往前100步，预测极限为0.95的预测结果，可以看出，预测的结果（红色部分）在预测极限（绿色部分）之内，且随着预测的步数变大，其与已知数据的相关性减弱，就给出了平稳模型的均值。可以从图中看出，红线部分为预测的100步数据，在预测30步之后，预测值给出的数据基本是该ARMA模型的均值。

5.2模型更新

5.1节中我们根据模型预测了100部的数据，倘若在实际运用中，随着时间的流逝，我们又有了新的观测数据，比如说又知道了100步中前50 步数据，那么我们如果想要更向前预测，因为随着预测步数的增加，相关性减弱，我们就无法获得更合理的预测数据。因此，我们应该将新得到的50步真实观测数据，更新到我们的模型数据中去，即将我们新得到的50步数据加到我们数据序列的末端。与此同时，我们要删除掉原始序列中的前50个历元数据，因为如果不进行删除老数据，随着时间的流逝，数据会越来越多，如果数据过多不仅会造成模型的不稳定，而且还会降低预测的速度。然后再根据模型进行下一次的预测，不断循环，不断预测。数据的更新过程，在原始数据中手动改动。更新后预测结果见图表5-2.

图表5-2 模型更新50步后重新预测50步

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1-1107_new.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,2),method="ML")

prop.fore=predict(model,n.ahead=50)

Up=prop.fore$pred+1.96*prop.fore$se

Low=prop.fore$pred-1.96*prop.fore$se

ts.plot(prop,prop.fore$pred,col=1:2,xlab="time",ylab="Observation: Z")

lines(Up,col="green",lty="dashed")

lines(Low,col="green",lty="dashed")

6.剩余数据的建模分析

由于整个原始数据中前1200个历元是震前的正常数据，中间1207-1700历元则是地震时候的波峰数据，1701-2400历元时地震减弱，但还未恢复正常的数据，因为数值的震荡区间差距太大，所以进行了分段建模分析。

6.1 中间地震数据的建模

由于中间段1200-1550历元在震荡范围上还是相差很大，见图表6-1，在经过多次尝试之后还是没能找到合适的模型，因此中间段的地震数据就没有建模成功。如果继续尝试，个人觉得还的进行分段建模。

图表 6-1 中间地震数据时间序列

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\1207-1700.txt",header=F)

prop=ts(mydata)

plot(prop,xlab="time",ylab="Observation: Z")

6.2 第三段数据的建模和预测

6.2.1 建模阶段

在第一段数据建模的过程中我们做出了详细的介绍，因此在接下来建模的过程中就不在详细论述，流程是一致的。

第三段数据是1701-2300历元，2300-2400历元做后期的模型更新使用，故作保留。第三段数据，是地震后恢复到正常的阶段，数据的振荡幅度已经远低于波峰时段，且整个时段表现出了均值平稳，没有明显趋势的特征，其时间序列见图表6-2。

图表6-2 第三段数据时间序列

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300.txt",header=F) prop=ts(mydata)

plot(prop,xlab="time",ylab="Observation: Z")

图表6-3 第三段数据的ACF和PACF

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300.txt",header=F) prop=ts(mydata)

acf(prop)

pacf(prop)

从图表6-3中可以看出，第三段数据的ACF和PACF都表现出了很好的拖尾性，再加上图表6-2中，第三段数据的时间序列没有明显的上升或是下降的趋势、均值稳定，我们暂定为ARMA模型。

然后模型参数（p,q）从低阶到高阶尝试，并不断用残差检验的方法来检验模型参数的合适与否，最终得出模型参数p=5,q=4。残差检验结果见图表6-4.

图表6-4 第三段模型残差检验结果

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,4),method="ML")

tsdiag(model)

r= model $residuals

Box.test(r,type="Ljung-Box",lag=6,fitdf=1)

从图表6-4可以看出，当模型暂定ARMA(5,4)时，残差的时间序列呈现没有任何趋势，均匀分布在0值上下，符合模型参数检验的第一个条件；从残差的自相关函数ACF和box_Ljung检验的结果p=0.4835>0.05，说明残差有很好的独立性，符合模型参数检验的第二个条件。

同时根据图表6-5，QQ图中的点的分布非常接近一条直线，只有在两段的时候有一些偏差，大部分数据分布很好。该QQ图说明残差有很好的正态性。满足残差检验的第三个条件，因此我们接受模型ARMA(5,4)。

图表6-5 模型ARMA(5,4)的残差分位数-分位数图

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,4),method="ML")

r= model $residuals

qqnorm(model $residuals)

qqline(model $residuals)

模型参数估计结果：

我们接受了ARMA(5,4)模型，其模型参数的估计结果如图表6-6所示。

图表6-6 模型参数估计结果

AR1 AR AR AR AR MA1 MA MA MA Intercept value 0.9311 -0.2109 -0.5594 0.4890 -0.3017 -0.4260 -1.224 0.3360 0.3676 -717.2480 s.e. 0.5297 0.9398 0.9060 0.4888 0.1784 0.5298 0.672 0.0992 0.2456 433.7555 mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,4),method="ML")

model

6.2.2 模型的预测和更新

首先，第三段数据的模型确定下来为ARMA(5,4)模型，接下来我们根据1701-2300历元的已知数据根据模型来预测接下来60步的数据，预测结果见图表6-7.

图表6-7 往后预测60步时间序列

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300.txt",header=F)

prop=ts(mydata)

model=arima(prop,order=c(5,0,4),method="ML")

prop.fore=predict(model,n.ahead=60)

Up=prop.fore$pred+1.96*prop.fore$se

Low=prop.fore$pred-1.96*prop.fore$se

ts.plot(prop,prop.fore$pred,col=1:2,xlab="time",ylab="Observation: Z")

lines(Up,col="green",lty="dashed")

lines(Low,col="green",lty="dashed")

可以从图表6-7往后预测60步的结果可以看出，红色部分为预测结果，绿色部分为0.95的预测极限，该模型预测的结果良好。最初几步很好体现了原始数据的特征，在均值上下来回振荡，但是随着预测步数的增加，相关性逐渐减弱，到一定程度之后，预测给出的结果就会是数据的均值。因此，根据原始数据和所确定的模型，只能预测未来很近的数据。因此，随着时间的流逝，我们会获得新的真实观测数据，因此我们需要对模型的原始数据进行更新，进行下一步的预测，这样才能保证预测结果与真实数据的贴近性，才能有效的利用预测数据去为决策者服务。

假如随着时间的流逝，我们获得了2300历元之后的30个真实观测数据，如果我们继续用我们预测得到的30步数据去进行下一步的预测，显然是不合适，因此我们就需要对原始数据进行更新，这样才能保证模型数据的真实性以及下一步预测结果的可靠性。我们需要将新获得的30步数据添加到原始数据的后面，同时将原始数据最开始的30部删去，这样才做到了数据的更新。原始数据的更新直接在数据中手动操作。图表6-8 给出了原始数据更新30步之后，又重新预测30步的结果。

图表6-8 更新30步再预测30步的结果

mydata=read.table("E:\\课程学习\\时间序列分析与应用\\建模资料\\ 1701-2300_new.txt",header=F) prop=ts(mydata)

model=arima(prop,order=c(5,0,4),method="ML")

prop.fore=predict(model,n.ahead=30)

Up=prop.fore$pred+1.96*prop.fore$se

Low=prop.fore$pred-1.96*prop.fore$se

ts.plot(prop,prop.fore$pred,col=1:2,xlab="time",ylab="Observation: Z")

lines(Up,col="green",lty="dashed")

lines(Low,col="green",lty="dashed")

比较图表6-8和图表6-7的预测结果，我们很明显得看出，更新之后再进行预测的30步数据比不更新数据直接预测60步中的后30步更能表现该时间序列的总体特征，更加可靠，因此随着时间的推迟，新数据的获取，我们要及时对数据进行更新。

7.总结

本文根据某一地震的观测数据，具体数为UEN-A-023.txt在R语言的编辑环境下，对该是数据时间序列的平稳性做出了分析。因为模型的建立一般是针对有规律的事件，根据已有的观测数据，用一个模型去拟合出它的发展规律，然后再利用建立的模型和已知的数据对未来进行一定程度的预测。但是经过分析，该数据是地震数据，前半段是正常的观测数据，中间突然发生了地震，持续一段时间后地震减弱，逐渐恢复的过程。对于一个正常的观测数据，地震的发生本来就是异常，并且持续的时间很长，占据了整个数据一半，这样以来，利用干预分析就很难去掉地震这个异常情况，经过数次尝试以后，得出结论，该数据不适合整体建模。因此，本文采用了分段建模的防范，对第一段和第三段时间序列进行了平稳性分析、模型的识别、参数的估计以及预测和模型的更新等过程。为以后时间序列建模分析提供参考。

参考文献

[1]Jonathan D.Cryer,Kung-Sik Chan.时间序列分析及应用R语言.第二版.北京。机械工业出版社.2011.1

[2].张翠娟，冯学军，盛敏.因子分析开发步骤及R语言代码实现.安庆师范学院计算机与信息学院.2013.2

[3]. 张美英,何杰.时间序列预测模型研究简介.江西科学.2009.5