《R语言在时间序列中的应用要点.pdf》由会员分享，可在线阅读，更多相关《R语言在时间序列中的应用要点.pdf（13页珍藏版）》请在三一文库上搜索。
1、- 1 - 时间序列分析在人口预测问题中的应用 摘要 时间序列分析是研究动态数据的动态结构和发展变化规律的统 计方法。以 1949 年至 2004 年中国大陆人口自然增长率为例, 用时 间序列分析和统计学软件R 建立模型 , 并对人口进行预测 , 取得较 好的效果。说明时间序列分析在人口预测问题上是有效的。 关键词 : ARMA 模型; R 软件; 平稳性; 可逆性 - 2 - Application of time series analysis in population prediction Abstract Time series analysis is a statistic met
2、hod studying dynamic structure of dynamic data and the law of de-velopment and change. Based on the example of population growth rate between 1949 and 2004 in the mainlandof China, mathematic models were established with time series analysis method and statistic software R,and population was predict
3、ed with it. It received a good result. Therefore the application of time series analysisis effective in population prediction. Key words:ARMA model; R software; stability; invertibility - 3 - 一时间序列概述 1. 概念 所谓时间序列就是按照时间的顺序记录的一列有序数据。对时间序列进行观 察、研究，找寻它变化发展的规律， 预测它将来的走势就是时间序列分析时间序 列分析有着非常广泛的应用领域。 2. 定义 在
4、统计研究中，常用按时间序列排列的一组随机变量 , 12 , t XXX, 来表示一个随机事件的时间序列，简记为, t X tT 或 t X。 3. 主要分析方法 时间序列分析方法主要有描述性时序分析和统计时序分析。 描述性时序分析主要通过直观数据比较或绘图测绘，统计时序分析主要有频域分 析方法以及时域分析方法。 常用的是时域分析法， 时域分析法的基本思想是源于 事件的发展通常具有一定的惯性， 这种惯性用统计语言来描述就是序列值之间存 在一定的相关关系， 而这种关系具有某种统计规律。 我们分析的重点就是找寻这 种规律，选取合适的数学模型拟合，进而预测该事件发展走向。 4. 研究意义 事件序列分析
5、具有现实意义，在金融经济、气象水文、信号处理、机械振动 等众多领域具有广泛的应用。 二时间序列的预处理 通常得到一个观察值序列后首先要对其进行平稳性以及纯随机性进行检验。 根据检验结果的不同我们有不同的处理方法 1. 平稳性 时间序列的平稳性分为严平稳与宽平稳 (1) 严平稳定义 设 t X一时间序列。对任意整数m，任取 12 , m t tt, T ，对任意整数， 有 1212 ,12,12 (,)(,) mm t ttmtttm Fx xxFx xx，则称序列 t X为严稳序列。 - 4 - 其中 F 为分布函数。 (2) 宽平稳定义 如果 t X满足：任取tT，有 2 t EX； 任取t
6、T，有 t EX，为常数； 任取, ,t s kT，且ktsT，有( , )( ,)t sk kst； 则称 t X为宽平稳序列。其中( , )t s表示 t X 与 s X 的自相关系数。 (3) 平稳性的检验 平稳性检验主要有时序图检验以及自相关图检验。 2. 纯随机性 (1) 纯随机性定义 如果时间序列 t X 满足以下性质： 任取tT，有 t EX，为常数； 任取, t sT，有 2, ( , ) 0, ts t s ts 则称序列为纯随机序列，也称为白噪声(white noise) 序列。 (2) 纯随机性检验 构造检验统计量，主要是Q统计量以及 LB统计量。 三时间序列分析的主要方
7、法及模型 1. 平稳时间序列分析的模型 (1) AR 模型( auto regression model ) 具有如下结构的模型称为p 阶自回归模型，记为AR(p)： 011 2 0 ()0,(),()0, 0, ttptpt p ttts st xxx EVarEst Exst (2) MA 模型(moving average) - 5 - 具有如下结构的模型称为q 阶移动平均模型，记为MA(q) ： 1122 2 0 ()0,(),()0, ttttqt q q ttts x EVarEst (3)ARMA 模型(auto regression moving average) 具有如下结构
8、的模型称为自回归移动平均模型，记为ARMA(p,q) ： 01111 2 0,0 ()0,(),()0, 0, ttptpttqtq pq ttts st xxx EVarEst Exst 若 0 0，该模型称为中心化ARMA(p,q) 模型。 2. 非平稳序列分析 事实上在自然界中绝大部分序列都是非平稳的，因而对非平稳序列的分析更 普遍更重要。 对非平稳时间序列的分析法通常分为确定性时序分析和随机时序分析。这里 简要介绍常用确定性时序分析方法。 (1) 趋势分析 有些时间序列具有非常显著的趋势，我们分析的目的就是要找到序列中的这 种趋势，并利用这种趋势对序列对序列的发展做出合理的预测。 (2
9、) 季节效应分析 在日常生活中我们可以看到许多有季节效应的时间序列，如四季气温等等。 凡是呈现出固定的周期性变化的时间，我们都称其有季节效应。 (3) 综合分析 既有趋势起伏变动又有季节效应的复杂序列的分析方法，常用模型有： 加法模型 tttt xTSI 乘积模型 tttt xT S I 混合模型 - 6 - a. tttt xT SI b.() tttt xTSI 式中， t T代表序列的长期趋势波动； t S 代表序列的季节性（周期性）变化； t I代 表随机波动。 3. 非平稳序列的模型 事实上，许多非平稳序列差分后会显示出平稳序列的性质，称之为差分平稳 序列。对差分平稳序列可以用ARI
10、MA 模型拟合。 具有如下结构的模型称为求和自回归移动平均(autoregressive integrated moving average) 模型，简记为 ARIMA(p,d,q) 模型： 2 ( )( ) ()0,(),()0, ()0, d tt ttts st BxB EVarEst E xst 式中： (1) dd B； 1 ( )1 p p BBB ，为平稳可逆ARMA(p,q) 模型的自回归系 数多项式； 1 ( )1 q q BBB ，为平稳可逆 ARMA(p,q) 模型的移动平滑系 数多项式； t 为零均值白噪声序列。 由上式可知 ARIMA 模型的实质就是差分运算与ARMA
11、 模型的组合。 当序列具有非常显著的确定性趋势或季节效应时，人们会怀念确定性因素分解方 法对各种确定性效应的解释， 但又因为它对残差信息的浪费而不敢轻易使用。为 了解决这个问题人们构造了残差自回归(auto-regressive) 模型。 Auto-Regressive 模型的构造思想是首先通过确定性因素分解方法提取序列中主 要的确定性信息： tttt xTS 式中， t T为趋势效应拟合， t S 为季节效应拟合。 考虑到因素分解方法对确定性信息的提取可能不够充分，因而需要进一步检 验残差序列 t 的相关性。如果检验结果显示残差序列自相关性不显著，说明确 定性回归模型对信息提取比较充分，可以
12、停止分析。 如果检验结果显示残差序列 - 7 - 自相关性显著，这时可以考虑对残差拟合自回归模型，进一步提取相关信息： 11ttptpt a 这样构造的模型： 11 2 ()0,(),(,)0,1 tttt ttptpt tttti xTS a E aVar aCov a ai 称为残差自回归模型。 四实例 本文以中国大陆人口自然增长率( 19492004 年)为样本进行分析 , 数据(数 据来自 2005 年统计年鉴 ) 见表 1。 表 1 中国大陆人口自然增长率 年份人口自然增长率 % 1940 16 1950 19 20 20 23 24 20.32 20.5 23.23 17.24 1
13、0.19 1960 -4.74 3.78 26.99 33.33 27.64 28.38 26.22 25.53 27.38 26.08 1970 25.83 23.33 22.16 20.89 17.48 15.69 12.66 12.06 12 11.61 1980 11.87 14.55 15.68 13.29 13.08 14.26 15.57 16.61 15.73 15.04 1990 14.39 12.98 11.6 11.45 11.21 10.5 10.42 10.06 9.14 8.18 2000 7.58 6.95 6.45 6.01 5.87 第一步：原数据的平稳性检验
14、 为判断一个序列是否平稳， 我们主要通过时序图以及自相关图进行检验。因 为用到 ARIMA 模型的拟合和检验，所以在程序的开头会载入tseries 。 首先绘出时序图、自相关图、偏自相关图（如下）： - 8 - Time V 1 01020304050 0 1 0 2 0 3 0 051015 - 0 .2 0 . 2 0 . 6 1 . 0 Lag A C F V1 51015 - 0 .4 0 .0 0 . 4 0 .8 Lag P a r t ia l A C F Series pr 图 1 原数据的时序图、自相关图、偏自相关图 通过观察时序图，序列有递减趋势，所以我们基本可以判断该序列
15、非平稳。 第二步：差分并检验新序列的平稳性，完成序列的定阶 Time V 1 1020304050 - 1 0 0 1 0 2 0 051015 - 0 . 5 0 . 0 0 .5 1 .0 Lag A C F V1 51015 - 0 .4 0 .0 0 . 2 Lag P a r ti a l A C F Series e - 9 - 图 2 1 阶差分后的时序图、自相关图、偏自相关图 为了增加说服力， 我们再次使用单位根检验对1 阶差分后的的序列进行验证。 结 果如下： Augmented Dickey-Fuller Test data: d Dickey-Fuller = -4.89
16、92, Lag order = 3, p-value = 0.01 alternative hypothesis: stationary Warning message: In adf.test(d) : p-value smaller than printed p-value 从上可以看出， p 值等于 0.01，远远小于 0.05，因此拒绝原假设。所以可以 认定 1 阶差分后的序列基本平稳，模型定阶结束。 第三步：拟合模型并通过AIC 准则检验 在 ARIMA 模型拟合过程中中，模型为一阶差分，p 和 q 分别取 0 到 3 并依 次验证，取 AIC 值最小者。 P Q AIC 0 1 3
17、08.89 0 2 309.62 0 3 303.79 1 0 319.71 1 1 310.5 1 2 304.96 1 3 305.31 2 0 305.87 2 1 307.67 2 2 305.57 2 3 307.14 3 0 307.83 3 1 308.85 - 10 - 3 2 307.41 3 3 309.12 从上面的表格可以看出，当p 取 0，q 取 3 的时候， AIC 函数的值最小。因 此，模型初步定为ARIMA （0,1,3） 。 第四步：残差检验 Time r 01020304050 - 1 5 - 1 0 - 5 0 5 1 0 图 3 残差的时序图 Box-P
18、ierce test data: r X-squared = 0.0017, df = 1, p-value = 0.9673 由于 p 值远大于 0.05， 则该序列为白燥声， 故模型最终定为ARIMA（0,1,3） 。 第五步：画出 qq图 - 11 - -2-1012 - 1 5 - 1 0 - 5 0 5 1 0 Normal Q-Q Plot Theoretical Quantiles S a m p le Q u a n ti le s 从 qq 图的效果可以看出，二者还是比较吻合的，但是右侧头部和左侧尾部 偏离期望的正态分布。 第六步：预测并比较 结果如下： $pred Time
19、 Series: Start = 57 End = 58 Frequency = 1 1 6.646174 7.601378 $se Time Series: Start = 57 End = 58 Frequency = 1 - 12 - 1 3.521358 6.133766 因此我的未来 2 年预测结果分别是6.646174， 7.601378。 原论文的预测结果为： 两个预测结果进行比较， 略有差距， 其中的原因可能是拟合的问题，也可能 是软件的差异。 参考文献 1王燕 应用时间序列分析 . 北京：中国人民大学出版社，2005 2 何书元 应用时间序列分析 . 北京：北京大学出版社，2
20、003 3 Paul Teetor R语言经典实例 . 北京：机械工业出版社，2013 年份2003 2004 实际值6.01 5.87 预测值6.16626 5.8761 - 13 - 完整的程序如下： library(tseries) pr-read.table(D:1.txt) par(mfrow=c(2,2);plot.ts(pr);acf(pr);pacf(pr) d=diff(ts(pr) par(mfrow=c(2,2);plot.ts(d);acf(d);pacf(d) adf.test(d) nihe=arima(pr,order=c(0,1,1),method=“ML“) nihe r=nihe$residuals plot.ts(r) Box.test(r) qqnorm(r) qqline(r) pr.fore = predict(arima(pr, order = c(0,1,3), n.ahead =2) pr.fore