R语言在时间序列中的应用

1、时间序列分析在人口预测问题中的应用摘要时间序列分析是研究动态数据的动态结构和发展变化规律的统计方法。以1949年至2004年中国大陆人口自然增长率为例，用时 间序列分析和统计学软件 R建立模型，并对人口进行预测，取得较 好的效果。说明时间序列分析在人口预测问题上是有效的。关键词：ARMA模型;R软件；平稳性；可逆性-1 -Application of time series analysis in populationpredictionAbstractTime series an alysis is a statistic method study ing dyn amic structur

2、e of dyn amic data and the law of de-velopme nt and cha nge. Based on the example of populati on growth rate betwee n 1949 and 2004 in the mainlan dof China, mathematic models were established with time series an alysis method and statistic software R,a nd populatio n was predicted with it. It recei

3、ved a good result. Therefore the applicati on of time series an alysisis effective in populati on predict ion.Key words: ARMA model; R software; stability; invertibility-2 -.时间序列概述1. 概念所谓时间序列就是按照时间的顺序记录的一列有序数据。对时间序列进行观 察、研究，找寻它变化发展的规律，预测它将来的走势就是时间序列分析时间序 列分析有着非常广泛的应用领域。2. 定义在统计研究中，常用按时间序列排列的一组随机变量X1

4、 , X2，,Xt ,来表示一个随机事件的时间序列，简记为Xt，tT或Xt。3. 主要分析方法时间序列分析方法主要有描述性时序分析和统计时序分析。描述性时序分析主要通过直观数据比较或绘图测绘，统计时序分析主要有频域分 析方法以及时域分析方法。常用的是时域分析法，时域分析法的基本思想是源于 事件的发展通常具有一定的惯性，这种惯性用统计语言来描述就是序列值之间存 在一定的相关关系，而这种关系具有某种统计规律。我们分析的重点就是找寻这 种规律，选取合适的数学模型拟合，进而预测该事件发展走向。4. 研究意义事件序列分析具有现实意义，在金融经济、气象水文、信号处理、机械振动等众多领域具有广泛的应用。二.

5、 时间序列的预处理通常得到一个观察值序列后首先要对其进行平稳性以及纯随机性进行检验。根据检验结果的不同我们有不同的处理方法1. 平稳性时间序列的平稳性分为严平稳与宽平稳（1）严平稳定义设Xt 一时间序列。对任意整数 m，任取 H，嚅T，对任意整数 ， 有Fg,tm（X1，X2，Xm）" ；,tm .（为必，,Xm），则称序列Xt为严稳序列。其中F为分布函数。(2) 宽平稳定义如果Xt满足：任取tT ,有EX：:; 任取r T ,有EXt 二，为常数； 任取 t,s,k T，且 k ts T，有(t,s)二(k,k s t)；则称Xt为宽平稳序列。其中(t,s)表示Xt与Xs的自相关系

6、数。(3) 平稳性的检验平稳性检验主要有时序图检验以及自相关图检验。2. 纯随机性(1) 纯随机性定义如果时间序列Xt满足以下性质： 任取t T，有EXt = i 为常数;任取t,sT ,有卢2丄a ,t=s_0,t 式 s则称序列为纯随机序列，也称为白噪声 (white noise)序列。(2) 纯随机性检验构造检验统计量，主要是 Q统计量以及LB统计量。三. 时间序列分析的主要方法及模型1. 平稳时间序列分析的模型(1) AR 模型(auto regression mode)l具有如下结构的模型称为p阶自回归模型，记为AR(p):Xt = J 必4 pXy；tp = 0EG)=0,Var(

7、q)M；, E(吓s)=0,stlEXst =0,%彳 MA 模型(moving average)具有如下结构的模型称为q阶移动平均模型，记为MA(q):Xt =卩 + 日1 ®t 二一日2 % _2 一 _q“02E(eJ =0,Var(st)"器,(翳民)=0,s芒t(3) ARMA 模型(auto regression moving average)具有如下结构的模型称为自回归移动平均模型，记为ARMA(p,q):Xt = °0 + °1人二 + °P Xt _p + 色一日1 色 4 一一日q 色 _q% H0,8q 芒 02Eg) =

8、0,Var(Et)"器,(可是)=0,sKtExsEt =0,7s <t若0 =0，该模型称为中心化 ARMA(p,q)模型。2. 非平稳序列分析事实上在自然界中绝大部分序列都是非平稳的，因而对非平稳序列的分析更 普遍更重要。对非平稳时间序列的分析法通常分为确定性时序分析和随机时序分析。这里 简要介绍常用确定性时序分析方法。(1) 趋势分析有些时间序列具有非常显著的趋势，我们分析的目的就是要找到序列中的这 种趋势，并利用这种趋势对序列对序列的发展做出合理的预测。(2) 季节效应分析在日常生活中我们可以看到许多有季节效应的时间序列，如四季气温等等。 凡是呈现出固定的周期性变化的时

9、间，我们都称其有季节效应。(3) 综合分析既有趋势起伏变动又有季节效应的复杂序列的分析方法，常用模型有： 加法模型 乘积模型X 二T S h混合模型a.xt =Tt Sth b.xt =T (St It)式中，Tt代表序列的长期趋势波动；St代表序列的季节性(周期性)变化；h代 表随机波动。3. 非平稳序列的模型事实上，许多非平稳序列差分后会显示出平稳序列的性质，称之为差分平稳序列。对差分平稳序列可以用 ARIMA模型拟合。具有如下结构的模型称为求和自回归移动平均(autoregressive in tegrated movi ngaverage模型，简记为 ARIMA(p,d,q)模型：=O

10、(B)勺2E( ；t) =0,Var( ；t) 乂 E( ；t ；s) =0,s"E(Xs4) =0,Ws vt式中：id =(1-B)d(B) =1-止-BP，为平稳可逆ARMA(p,q)模型的自回归系 数多项式；0(B) =17B-，qBq，为平稳可逆ARMA(p,q)模型的移动平滑系 数多项式； .为零均值白噪声序列。由上式可知ARIMA模型的实质就是差分运算与 ARMA模型的组合。当序列具有非常显著的确定性趋势或季节效应时，人们会怀念确定性因素分解方 法对各种确定性效应的解释，但又因为它对残差信息的浪费而不敢轻易使用。为 了解决这个问题人们构造了残差自回归(auto-regr

11、essive模型。Auto-Regressive模型的构造思想是首先通过确定性因素分解方法提取序列中主要的确定性信息：Xt 珂 St t式中，Tt为趋势效应拟合，St为季节效应拟合。考虑到因素分解方法对确定性信息的提取可能不够充分，因而需要进一步检验残差序列 ；t的相关性。如果检验结果显示残差序列自相关性不显著，说明确 定性回归模型对信息提取比较充分， 可以停止分析。如果检验结果显示残差序列自相关性显著，这时可以考虑对残差拟合自回归模型，进一步提取相关信息：；t = itat这样构造的模型：人=Tt +S十勺；t = 1 2 t_p at2口耳）=O,Var（at） = ,Cov（at,耳丄）

12、=0,E称为残差自回归模型。四. 实例本文以中国大陆人口自然增长率（19492004年）为样本进行分析,数据（数 据来自2005年统计年鉴）见表1。表1中国大陆人口自然增长率年份人口自然增长率%1940161950192020232420.3220.523.2317.2410.191960-4.743.7826.9933.3327.6428.3826.2225.5327.3826.08197025.8323.3322.1620.8917.4815.6912.6612.061211.61198011.8714.5515.6813.2913.0814.2615.5716.6115.7315.041

13、99014.3912.9811.611.4511.2110.510.4210.069.148.1820007.586.956.456.015.87第一步：原数据的平稳性检验为判断一个序列是否平稳，我们主要通过时序图以及自相关图进行检验。因为用到ARIMA模型的拟合和检验，所以在程序的开头会载入 tseries。首先绘出时序图、自相关图、偏自相关图（如下）：-8 -V1TimeSeries prFC Allrtrap-9 -# -图1原数据的时序图、自相关图、偏自相关图通过观察时序图，序列有递减趋势，所以我们基本可以判断该序列非平稳第二步：差分并检验新序列的平稳性，完成序列的定阶-# -# -T

14、ox- 50 oo 5OFATimeSeries eFCAPLaLrap-10 -图21阶差分后的时序图、自相关图、偏自相关图为了增加说服力，我们再次使用单位根检验对1阶差分后的的序列进行验证。结 果如下：Augme nted Dickey-Fuller Testdata: dDickey-Fuller = -4.8992, Lag order = 3, p-value = 0.01alter native hypothesis: stati onaryWarning message:In adf.test(d) : p-value smaller tha n prin ted p-value

15、从上可以看出，p值等于0.01,远远小于0.05，因此拒绝原假设。所以可以 认定1阶差分后的序列基本平稳，模型定阶结束。第三步：拟合模型并通过AIC准则检验在ARIMA模型拟合过程中中，模型为一阶差分，p和q分别取0到3并依 次验证，取AIC值最小者。PQAIC01308.8902309.6203303.7910319.7111310.512304.9613305.3120305.8721307.6722305.5723307.1430307.8331308.8532307.4133309.12从上面的表格可以看出，当p取0, q取3的时候，AIC函数的值最小。因 此，模型初步定为ARIMA

16、（0,1,3）。第四步：残差检验图3残差的时序图Box-Pierce testdata: rX-squared = 0.0017, df = 1, p-value = 0.9673由于p值远大于0.05,则该序列为白燥声，故模型最终定为ARIMA（ 0,1,3）第五步：画出qq图se.tnauQ orpmasNormal Q-Q Plot-2-1012Theoretical QuantilesnV5 o5-HV_从qq图的效果可以看出，二者还是比较吻合的，但是右侧头部和左侧尾部 偏离期望的正态分布。第六步：预测并比较结果如下：$predTime Series:Start = 57End = 5

17、8Freque ncy = 11 6.646174 7.601378$seTime Series:Start = 57End = 58Freque ncy = 11 3.521358 6.133766因此我的未来2年预测结果分别是6.646174,7.60137&原论文的预测结果为：年份20032004实际值6.015.87预测值6.166265.8761也可能两个预测结果进行比较，略有差距，其中的原因可能是拟合的问题, 是软件的差异。参考文献1 王燕 应用时间序列分析.北京：中国人民大学出版社，20052 何书元 应用时间序列分析.北京：北京大学出版社，2003Paul Teetor R语言经典实例.北京：机械工业出版社，2013-14 -完整的程序如下:library(tseries)pr<-read.table('D:1.txt')par(mfrow=c(2,2);plot