R语言建模和预测

R语言下时间序列建模和预测研究摘要：时间序列是指将某种现象某一个统计指标在不同时间上的各个数值，按时间先后顺序排列而形成的序列。研究时间序列的领域多的难以罗列，但是时间序列分析的目的一般有两个方面：一是认识产生观测序列的随机机制，即建立数据生成模型；二是基于序列的已知数据，对序列未来的可能取值给出预测或预报。本文在分析和理解平稳时间序列和非平稳时间序列的基础上，在R语言环境下，根据某一地震观测数据的时间序列，实现了对该时间序列的平稳性分析，模型识别，参数估计，模型诊断，预测和更新。关键字：时间序列，平稳性分析，模型识别，参数估计，预测更新1.引言时间序列是指将某种现象某一个统计指标在不同时间上的各个数值，按时间先后顺序排列而形成的序列。随着现代社会的发展，数据的产量越来越大，对于数据时间序列的分析研究越来越多，时间序列的研究应用于很多领域，经济上用于统计和预测经济增长等各项指标，在天气预报上，也通过对云层的观测序列，预测未来的天气情况，在测量上的应用更是比比皆是。时间序列模型，可以根据其平稳性可以分为平稳时间序列模型和非平稳时间序列模型。而平稳时间序列模型又可以分为一般线性模型、滑动平均模型、自回归模型和自回归滑动混合模型;而非平稳模型则是多种多样。2.模型概述2.1平稳时间序列模型平稳时间序列满足期望为零，均值在所有时间上均为常数，且任意两个时刻的相关函数与时间t无关，仅与两个时刻的时间差相关，平稳时间序列模型包括一般线性过程、滑动平均过程、自回归过程、自回归滑动平均过程[1]。2.1.1一般线性过程一般线性过程*𝑌𝑡+可以表示成现在和过去白噪声变量的加权线性组合：𝑌𝑡=𝑒𝑡+𝜑1𝑒𝑡−1+𝜑2𝑒𝑡−2+∙∙∙(2.1)如果表达式的右边事实上是一个无穷级数，那么需要给权数𝜑加上一定条件，这样右边的表达式才在数学上有意义，为此，只需假设∑𝜑𝑖2∞∞𝑖=1(2.2)2.1.2滑动平均模型当有限个系数𝜑不为零时，就得到了所谓的滑动平均过程，改变正负号，如下式所示：𝑌𝑡=𝑒𝑡−𝜃1𝑒𝑡−1−𝜃2𝑒𝑡−2−∙∙∙−𝜃𝑞𝑒𝑡−𝑞(2.3)(2.3)式为q阶的滑动平均过程，记作MA(q)。2.1.3自回归模型自回归过程就是用自身做回归变量，具体来说，p阶自回归过程*Yt+满足方程：𝑌𝑡=∅1𝑌𝑡−1+∅2𝑌𝑡−2+∙∙∙+∅𝑝𝑌𝑡−𝑝+𝑒𝑡(2.4)(2.4)式为p阶的自回归模型，记为AR(p),，序列Yt的当期值是自身最近p滞后项和新信息项𝑒𝑡的组合，其中包含了序列在t期无法用过去值来解释的所有新信息，因此对于每一个t，假设𝑒𝑡独立于Yt−1，Yt−2，Yt−3∙∙∙2.1.4自回归滑动平均混合模型如果假定序列中部分是自回归过程，部分是滑动平均过程，就可以得到一个相当于普遍的时间序列模型：𝑌𝑡=∅1𝑌𝑡−1+∅2𝑌𝑡−2+∙∙∙+∅𝑝𝑌𝑡−𝑝+𝑒𝑡−𝜃1𝑒𝑡−1−𝜃2𝑒𝑡−2−∙∙∙−𝜃𝑞𝑒𝑡−𝑞(2.5)(2.5)式是𝑌𝑡的自回归滑动平均模型，滑动平均部分的阶数为q,自回归部分的阶数为p,记为ARMA(p,q)。2.2非平稳时间序列模型具有时变均值的时间序列都是非平稳的，形如：Y𝑡=𝜇𝑡+𝑋𝑡(2.6)其中：μt为时变均值函数，Xt为平稳序列。该序列因为均值μt是不确定是不断变化的，因此该模型的时间序列趋势就是随机的，可能一直向上，一直向下，或是断断续续的上下颠簸，这样的模型不具备平稳性，平稳性模型是不适合的，必须有一个含有随机趋势的非平稳模型才是适当的，这就提出了ARIMA模型。如果一个时间序列*Yt+的d次差分W∇𝑑=Yt,Yt是一个ARMA的平稳模型，那么我们称*Yt+为自回归滑动平均求和模型，记为𝐴𝑅𝐼𝑀𝐴(𝑝,𝑑,𝑞),其中d为差分次数。可以看出，𝐴𝑅𝐼𝑀𝐴(𝑝,𝑑,𝑞)，可以表示平稳模型和非平稳模型，比如说：𝐴𝑅𝐼𝑀𝐴(𝑝,0,𝑞)模型就是自回归滑动平均模型𝐴𝑅𝑀𝐴(𝑝,𝑞)，𝐴𝑅𝐼𝑀𝐴(𝑝,0,0)就是自回归模型AR(p)，𝐴𝑅𝐼𝑀𝐴(0,0,𝑞)就是滑动平均模型𝑀𝐴(𝑞)。3.模型识别3.1基本知识在模型识别的过程中，具有两个非常重要的参数，分别为自相关函数𝜌𝑘=𝛾𝑘/𝛾0和偏自相关函数∅𝑘𝑘，如果某一时间序列的自相关函数随着滞后k的增加而很快地下降为0，那么我们就认为该序列为平稳序列；如果自相关函数不随着k的增加而迅速下降为0，就表明该序列不平稳。如果一个时间序列的自相关和偏相关图没有任何模式，而且数值很小,那么该序列可能就是一些互相独立的无关的随机变量。其中自相关函数与偏相关函数又具有两个非常重要的性质，分别为拖尾性和截尾性，其中拖尾性是指随着滞后k的增加，函数渐变为0，图像像是拖了一条尾巴；截尾性是指当滞后k＞p或者k＞q处，函数突然下降为0图像像是截断了尾巴一样。函数的拖尾性和截尾性对时间序列模型的识别具有重要作用，见表3-1：图表3-1ARMA模型ACF和PACF的一般特征AR(p)MA(q)ARMA(p,q),p0,q0ACF拖尾滞后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和PACFmydata=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-410次差分后的ACF和PACFmydata=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和PACFmydata=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模型后，就开始从低阶到高阶来估计模型参数，并不断进行残差检验，最终确定出了模型参数为𝑝=5，𝑞=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$residualsBox.test(model,type=Ljung-Box,lag=6,fitdf=1)由图表4-1残差检验结果可以看出，首先从残差的时间序列图中可以看出：残差是均匀分布在0值上下的，没有任何趋势的散点图，满足残差检验的条件一。从残差的自相关ACF可以看出，除了第13阶处有一点超限，其他的都在限差之内，在Ljung-Box检验中，p=0.5348＞5%，因此可以说残差有很好的独立性，满足残差检验的第二个条件。下面我们将在R预言下利用qqnorm和qqline函数画