R語言在時間序列中的應用

1、時間序列分析在人口預測問題中的應用摘 要時間序列分析是研究動態(tài)數(shù)據(jù)的動態(tài)結構和發(fā)展變化規(guī)律的統(tǒng)計方法。以 1949 年至 2004 年中國大陸人口自然增長率為例, 用時間序列分析和統(tǒng)計學軟件R建立模型, 并對人口進行預測, 取得較好的效果。說明時間序列分析在人口預測問題上是有效的。關鍵詞: ARMA 模型; R軟件; 平穩(wěn)性; 可逆性Application of time series analysis in population predictionAbstract Time series analysis is a statistic method studying dynamic str

2、ucture 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 predicted with it. It received

3、a good result. Therefore the application of time series analysisis effective in population prediction.Key words: ARMA model; R software; stability; invertibility一時間序列概述1. 概念所謂時間序列就是按照時間的順序記錄的一列有序數(shù)據(jù)。對時間序列進行觀察、研究，找尋它變化發(fā)展的規(guī)律，預測它將來的走勢就是時間序列分析時間序列分析有著非常廣泛的應用領域。2. 定義在統(tǒng)計研究中，常用按時間序列排列的一組隨機變量,來表示一個隨機事件的時間序列，

4、簡記為或。3. 主要分析方法時間序列分析方法主要有描述性時序分析和統(tǒng)計時序分析。描述性時序分析主要通過直觀數(shù)據(jù)比較或繪圖測繪，統(tǒng)計時序分析主要有頻域分析方法以及時域分析方法。常用的是時域分析法，時域分析法的基本思想是源于事件的發(fā)展通常具有一定的慣性，這種慣性用統(tǒng)計語言來描述就是序列值之間存在一定的相關關系，而這種關系具有某種統(tǒng)計規(guī)律。我們分析的重點就是找尋這種規(guī)律，選取合適的數(shù)學模型擬合，進而預測該事件發(fā)展走向。4. 研究意義事件序列分析具有現(xiàn)實意義，在金融經(jīng)濟、氣象水文、信號處理、機械振動等眾多領域具有廣泛的應用。二時間序列的預處理通常得到一個觀察值序列后首先要對其進行平穩(wěn)性以及純隨機性進行

5、檢驗。根據(jù)檢驗結果的不同我們有不同的處理方法1. 平穩(wěn)性時間序列的平穩(wěn)性分為嚴平穩(wěn)與寬平穩(wěn)(1)嚴平穩(wěn)定義設一時間序列。對任意整數(shù)，任取，對任意整數(shù)，有，則稱序列為嚴穩(wěn)序列。其中為分布函數(shù)。(2)寬平穩(wěn)定義如果滿足：任取，有；任取，有，為常數(shù)；任取，且，有；則稱為寬平穩(wěn)序列。其中表示與的自相關系數(shù)。(3)平穩(wěn)性的檢驗平穩(wěn)性檢驗主要有時序圖檢驗以及自相關圖檢驗。2. 純隨機性(1) 純隨機性定義如果時間序列滿足以下性質(zhì)： 任取，有，為常數(shù)； 任取，有則稱序列為純隨機序列，也稱為白噪聲(white noise)序列。(2) 純隨機性檢驗構造檢驗統(tǒng)計量，主要是Q統(tǒng)計量以及LB統(tǒng)計量。三時間序列分析的

6、主要方法及模型1. 平穩(wěn)時間序列分析的模型(1) AR模型(auto regression model)具有如下結構的模型稱為p階自回歸模型，記為AR(p)：(2) MA模型(moving average)具有如下結構的模型稱為q階移動平均模型，記為MA(q)：(3)ARMA模型(auto regression moving average)具有如下結構的模型稱為自回歸移動平均模型，記為ARMA(p,q)：若，該模型稱為中心化ARMA(p,q)模型。2. 非平穩(wěn)序列分析事實上在自然界中絕大部分序列都是非平穩(wěn)的，因而對非平穩(wěn)序列的分析更普遍更重要。對非平穩(wěn)時間序列的分析法通常分為確定性時序分析和

7、隨機時序分析。這里簡要介紹常用確定性時序分析方法。(1) 趨勢分析有些時間序列具有非常顯著的趨勢，我們分析的目的就是要找到序列中的這種趨勢，并利用這種趨勢對序列對序列的發(fā)展做出合理的預測。(2) 季節(jié)效應分析在日常生活中我們可以看到許多有季節(jié)效應的時間序列，如四季氣溫等等。凡是呈現(xiàn)出固定的周期性變化的時間，我們都稱其有季節(jié)效應。(3) 綜合分析既有趨勢起伏變動又有季節(jié)效應的復雜序列的分析方法，常用模型有： 加法模型 乘積模型 混合模型a. b.式中，代表序列的長期趨勢波動；代表序列的季節(jié)性（周期性）變化；代表隨機波動。3. 非平穩(wěn)序列的模型事實上，許多非平穩(wěn)序列差分后會顯示出平穩(wěn)序列的性質(zhì)，稱

8、之為差分平穩(wěn)序列。對差分平穩(wěn)序列可以用ARIMA模型擬合。具有如下結構的模型稱為求和自回歸移動平均(autoregressive integrated moving average)模型，簡記為ARIMA(p,d,q)模型：式中：；，為平穩(wěn)可逆ARMA(p,q)模型的自回歸系數(shù)多項式；，為平穩(wěn)可逆ARMA(p,q)模型的移動平滑系數(shù)多項式；為零均值白噪聲序列。由上式可知ARIMA模型的實質(zhì)就是差分運算與ARMA模型的組合。當序列具有非常顯著的確定性趨勢或季節(jié)效應時，人們會懷念確定性因素分解方法對各種確定性效應的解釋，但又因為它對殘差信息的浪費而不敢輕易使用。為了解決這個問題人們構造了殘差自回歸

9、(auto-regressive)模型。Auto-Regressive模型的構造思想是首先通過確定性因素分解方法提取序列中主要的確定性信息：式中，為趨勢效應擬合，為季節(jié)效應擬合?？紤]到因素分解方法對確定性信息的提取可能不夠充分，因而需要進一步檢驗殘差序列的相關性。如果檢驗結果顯示殘差序列自相關性不顯著，說明確定性回歸模型對信息提取比較充分，可以停止分析。如果檢驗結果顯示殘差序列自相關性顯著，這時可以考慮對殘差擬合自回歸模型，進一步提取相關信息：這樣構造的模型：稱為殘差自回歸模型。四實例 本文以中國大陸人口自然增長率(19492004年)為樣本進行分析,數(shù)據(jù)(數(shù)據(jù)來自2005年統(tǒng)計年鑒)見表1。

10、表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.04199014.3912.9811.611.4511.2110.510.4210.069.148.1820007.586.956.456.015.87

11、第一步：原數(shù)據(jù)的平穩(wěn)性檢驗為判斷一個序列是否平穩(wěn)，我們主要通過時序圖以及自相關圖進行檢驗。因為用到ARIMA模型的擬合和檢驗，所以在程序的開頭會載入tseries。首先繪出時序圖、自相關圖、偏自相關圖（如下）：圖1 原數(shù)據(jù)的時序圖、自相關圖、偏自相關圖通過觀察時序圖，序列有遞減趨勢，所以我們基本可以判斷該序列非平穩(wěn)。第二步：差分并檢驗新序列的平穩(wěn)性，完成序列的定階圖2 1階差分后的時序圖、自相關圖、偏自相關圖為了增加說服力，我們再次使用單位根檢驗對1階差分后的的序列進行驗證。結果如下： Augmented Dickey-Fuller Test data: d Dickey-Fuller = -

12、4.8992, 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階差分后的序列基本平穩(wěn)，模型定階結束。第三步：擬合模型并通過AIC準則檢驗在ARIMA模型擬合過程中中，模型為一階差分，p和q分別取0到3并依次驗證，取AIC值最小者。PQAIC01308.8902309.6203303.79103

13、19.7111310.512304.9613305.3120305.8721307.6722305.5723307.1430307.8331308.8532307.4133309.12從上面的表格可以看出，當p取0，q取3的時候，AIC函數(shù)的值最小。因此，模型初步定為ARIMA（0,1,3）。第四步：殘差檢驗圖3 殘差的時序圖 Box-Pierce test data: r X-squared = 0.0017, df = 1, p-value = 0.9673由于p值遠大于0.05，則該序列為白燥聲，故模型最終定為ARIMA（0,1,3）。第五步：畫出qq圖從qq圖的效果可以看出，二者還是比

14、較吻合的，但是右側頭部和左側尾部偏離期望的正態(tài)分布。第六步：預測并比較 結果如下： $pred Time Series: Start = 57 End = 58 Frequency = 1 1 6.646174 7.601378 $se Time Series: Start = 57 End = 58 Frequency = 1 1 3.521358 6.133766 因此我的未來2年預測結果分別是6.646174， 7.601378。 原論文的預測結果為：年份20032004實際值6.015.87預測值6.166265.8761兩個預測結果進行比較，略有差距，其中的原因可能是擬合的問題，也可能是軟件的差異。參考文獻1 王燕 應用時間序列分析.北京：中國人民大學出版社，20052 何書元 應用時間序列分析. 北京：北京大學出版社，20033 Paul Teetor R語言經(jīng)典實例. 北京：機械工業(yè)出版社，2013完整的程序如下： 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=