2013統(tǒng)計學實驗第15章時間序列分析

2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-1統(tǒng)計學實驗

—SPSS和R軟件應用與實例

主編：費宇2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-2第15章時間序列分析2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-3一、實驗目的1.運用R中plot()函數(shù)繪制時間序列圖，acf()和pacf()函數(shù)作時間序列數(shù)據(jù)的自相關和偏自相關圖檢驗序列的平穩(wěn)性，HoltWinters()濾波函數(shù)實現(xiàn)時間序列數(shù)據(jù)的指數(shù)平滑；一、實驗目的運用arima()函數(shù)進行ARIMA模型的識別、估計和預測，用pose()函數(shù)作時間序列的季節(jié)分解；掌握運用R軟件中的函數(shù)對時間序列數(shù)據(jù)進行分析的基本操作過程，并能讀懂R輸出的結果。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-42023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-5二、實驗環(huán)境系統(tǒng)軟件Windows2000或WindowsXP或Windows7；統(tǒng)計軟件R2.13.2或更高版本。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-6三、實驗內(nèi)容檢驗時間序列的平穩(wěn)性根據(jù)自相關系數(shù)和偏自相關系數(shù)來初步判斷ARMA模型的階數(shù)p和q對ARIMA模型進行估計、預測和診斷識別時間序列數(shù)據(jù)的季節(jié)變動，能看出其季節(jié)波動趨勢，學會剔除季節(jié)因素的方法2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-7第15章時間序列分析15.1時間序列的圖形化觀察15.2指數(shù)平滑15.3時間序列的Box-Jenkin

模型15.4時間序列的季節(jié)分解2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-815.1時間序列的圖形化觀察【例15.1】(數(shù)據(jù)文件為li15.1.txt)選取我國從1978年至2008年國內(nèi)生產(chǎn)總值(GDP，億元)數(shù)值，總31個觀測值，數(shù)據(jù)來源“國泰安宏觀經(jīng)濟數(shù)據(jù)庫”。試對該數(shù)據(jù)進行分析，判斷其是否具有平穩(wěn)性。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-9【統(tǒng)計理論】時間序列分析中一個重要的概念是平穩(wěn)性，它是時間序列建模的一個重要基礎。如果一個隨機過程的均值和方差在時間過程上都是常數(shù)，并且在任何兩時期的協(xié)方差僅依賴于這兩個時期間的間隔，而不依賴于計算這個協(xié)方差的實際時間，就稱它是寬平穩(wěn)的。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-10【統(tǒng)計理論】時間序列平穩(wěn)性檢驗的常用方法有：(1)時序圖依據(jù)均值和方差的統(tǒng)計意義，一個平穩(wěn)序列時間序列的實現(xiàn)大致是由在某一水平線附近等幅波動的點構成的。平穩(wěn)序列的時序圖應該顯示出序列始終圍繞一個常數(shù)值波動，且波動的范圍不大。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-11【統(tǒng)計理論】實踐中真正具有平穩(wěn)性的時間序列并不多見，通常都會表現(xiàn)出非平穩(wěn)性。非平穩(wěn)性的表現(xiàn)形式多種多多樣，主要特征有：趨勢性、異方差性、周期性等。在許多時間序列數(shù)據(jù)中，差分可以使時間序列達到平穩(wěn)性目的。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-12假設

為一時間序列，則該序列的一階差分定義為

，其中

表示一階差分算子。為了達到平穩(wěn)性的目的，一階差分是最常用的方法，有時我們甚至需要對其進行多次一階差分?！窘y(tǒng)計理論】2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-13【統(tǒng)計理論】

(2)自相關和偏自相關系數(shù)所謂自相關是指序列與其自身經(jīng)過某些階數(shù)滯后形成的序列之間存在某種程度的相關性。隨著

的增加，自相關函數(shù)下降且趨于0。從下降的速度看，平穩(wěn)序列要比非平穩(wěn)序列快得多。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-14【軟件操作】首先，用plot()函數(shù)繪制時間序列圖setwd("D:/R-Statistics/data/chap-15")#設定工作路徑x=read.table("li15.1.txt",header=T)#從li15.1.txt中讀入樣本數(shù)據(jù)xdx=diff(x[,2])#對DGP作一階差分變換【軟件操作】plot(x[,1],x[,2],type=“1",lty=1,col=1,xlab="時間",ylab="")lines(x[2:31,1],dx,lty=2,col=1)legend(1980,250000,c("GDP","GDP_1"),col=c(1,1),text.col="black",lty=c(1,2))2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-152023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-16【運行結果】圖15.1例15.1中GDP及GDP_1的時序圖【軟件操作】然后，采用acf()和pacf()函數(shù)作相關圖x=read.table("li15.1.txt",header=T)#從li15.1.txt中讀入樣本數(shù)據(jù)xGDP=x[,2]#從數(shù)據(jù)文件中讀取GDP數(shù)據(jù)acf(GDP)#GDP序列的自相關函數(shù)條形圖2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-17【運行結果】2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-18圖15.2例15.1中GDP的自相關函數(shù)的條形圖2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-1915.2指數(shù)平滑【例15.2】(數(shù)據(jù)文件為li15.2.txt)選取1990年至2004年我國月度社會消費品零售總額數(shù)據(jù)(億元)，利用指數(shù)平滑方法試對該數(shù)據(jù)進行分析。

2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-20【統(tǒng)計理論】當利用過去數(shù)據(jù)的加權平均預測將來的數(shù)值時(這一過程成為平滑)，對序列中較近的數(shù)據(jù)給予較大的權重，遠期的數(shù)據(jù)給于較小的權重。令

，則

的預測值為：即利用過去的觀測值進行加權平均預測時間序列在

時刻的數(shù)值。【軟件操作】使用HoltWinters()濾波函數(shù)可以作時間序列數(shù)據(jù)的指數(shù)平滑x=read.table("li15.2.txt",header=T)#從li15.2.txt中讀入樣本數(shù)據(jù)xtx=ts(x,frequency=12,start=c(1990,1),end=c(2004,8))#添加合適的時間格式2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-21【軟件操作】m=HoltWinters(tx,gamma=FALSE,beta=FALSE,seasonal="mult")#參數(shù)gamma=FALSE,濾波函數(shù)將作指數(shù)平滑，參數(shù)beta為Holt-Winters濾波函數(shù)的參數(shù)，參數(shù)seasonal為濾波函數(shù)中季節(jié)趨勢，seasonal="mult"表示季節(jié)效應為乘法模型。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-22【軟件操作】ts.plot(tx)#作時間序列tx的點圖lines(fitted(m)[,1],lty=2)#fitted(m)為擬合的序列l(wèi)egend(1990,4200,c("原序列","指數(shù)平滑后擬合值"),text.col="black",lty=c(1,2))2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-232023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-24【運行結果】圖15.3原序列和通過指數(shù)平滑后擬合值和觀測值的時間序列圖2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-2515.3Box-Jenkin模型

15.3.1ARMA模型的識別和估計【例15.3】(數(shù)據(jù)文件為li15.3.txt)該數(shù)據(jù)是由一個AR(1)模型模擬得到的數(shù)據(jù)，其中模擬參數(shù)，試對該數(shù)據(jù)進行模型識別和參數(shù)估計。

(數(shù)據(jù)來源：費宇等,《統(tǒng)計學》第7章,高等教育出版社,2010)2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-26【統(tǒng)計理論】(1)自回歸(AR)模型①模型形式假設

為觀測時間序列，

；對

的統(tǒng)計建模，一個自然的想法是用其過去某個時段上的數(shù)據(jù)對其進行回歸，寫成模型形式，即其中

為模型參數(shù)，

是相互獨立的均值為零，方差為

的隨機誤差。由于它是對自身變量的回歸，模型(15.4)稱為p階自回歸模型，記為。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-27【統(tǒng)計理論】②識別條件當使用ARMA模型擬合時間序列數(shù)據(jù)時，關鍵的問題是如何確定其階數(shù)，即即AR(p)中的p，MA(q)中的q，以及ARMA(p,q)中的p，q。通過某種方法確定這些階數(shù)的過稱成為ARMA模型的模型識別。識別模型階數(shù)的兩個重要統(tǒng)計量是自相關函數(shù)(autocorrelationfunction，簡記為ACF)和偏自相關函數(shù)(partialautocorrelationfunction，簡記為PACF)。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-28【統(tǒng)計理論】對于AR(p)模型，自相關函數(shù)

隨著p的增加呈現(xiàn)指數(shù)衰減或震蕩性衰減；而偏相關函數(shù)

是

階截尾的，即當時，

等于零。以AR(1)為例，此時

，由于，當

時，

呈現(xiàn)指數(shù)衰減，而當

時，

呈或震蕩性衰減。偏相關函數(shù)

，因此是一階截尾的。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-29【統(tǒng)計理論】假設

為后移(滯后)算子，即，

，令則(15.4)可以寫為：

2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-30【統(tǒng)計理論】有幾種信息準則也可用來決定AR過程的階p，它們都基于似然函數(shù)。例如，著名的AIC準則定義為：其中

是殘差平方和的估計。在實際應用中，事先給定一個正整數(shù)

，對

計算

，然后選擇階，使

達最小值。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-31【統(tǒng)計理論】③平穩(wěn)條件假設

表示

的一個

階多項式，則AR(p)模型平穩(wěn)的充分必要條件是：的根全部落在單位圓之外。最簡單的AR模型是AR(1)，模型為：

。此時，平穩(wěn)模型的充要條件是

。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-32【統(tǒng)計理論】④模型意義僅通過時間序列變量的自身歷史觀測值來反映有關因素對預測目標的影響和作用，影響序列變化的主要因素是時間序列在不同時期的取值。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-33【統(tǒng)計理論】(2)移動平均(MA)模型①模型形式對時間序列

，一個滑動平均模型假設

可以由該時刻和該時刻之前的隨機誤差線性表示，即

我們稱(15.6)為

階滑動平均模型，記為MA(q)。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-34【統(tǒng)計理論】滑動平均模型

(15.6)也可以表為，其中

。由于

只有有限項，因此式(15.6)總滿足平穩(wěn)性條件，是一個平穩(wěn)過程。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-35【統(tǒng)計理論】②識別條件對于MA(p)模型，自相關函數(shù)

是

階截尾的，即當

時，

接近于零。而偏相關系數(shù)

隨著的增加呈現(xiàn)指數(shù)衰減或震蕩性衰減，并逐漸趨于零。以MA(1)為例，此時

，

，即一階截尾；而偏相關系數(shù)

，因為

，偏相關系數(shù)將出現(xiàn)指數(shù)衰減。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-36【統(tǒng)計理論】③可逆條件MA模型中的一個重要內(nèi)容是研究它的可逆性。令

表示z的q階多項式，則MA(q)模型滿足可逆性的充分必要條件是：

的根全部落在單位圓之外。一階滑動平均模型MA(1)的可逆條件為：

。當滿足可逆條件時，MA(q)模型往往可以轉化為AR(p)模型。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-37【統(tǒng)計理論】④模型意義用過去各個時期的隨機干擾或預測誤差的線性組合來表達當前預測值。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-38【統(tǒng)計理論】(3)自回歸滑動平均(ARMA)模型①模型形式當模型既有自回歸項，又有滑動平均項時，就得到一個混合模型，稱為自回歸滑動平均模型：

ARMA(p,q)模型的平穩(wěn)條件為：的根全部落在單位圓之外。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-39【統(tǒng)計理論】②識別條件偏相關函數(shù)

和自相關函數(shù)

均不截尾，但較快收斂到0，則該時間序列可能是ARMA（p,q）。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-40【統(tǒng)計理論】

③模型階數(shù)ARMA模型階數(shù)的確定很難從ACF和PACF獲得充分的信息，TsayandTiao(1984)提出了一種稱為EACF(extendedautocorrelationfunction)方法來識別ARMA(p,q)模型階數(shù)。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-41【統(tǒng)計理論】在模型識別過程中，ACF及PACF能提供ARMA模型階數(shù)的重要信息，但在建模過程中，我們還應該結合模型擬合統(tǒng)計量，如

、AIC及BIC，進行綜合分析。其中一種選擇模型階數(shù)的方法是計算對應各種不同階數(shù)的AIC或BIC統(tǒng)計量，其中對應AIC或BIC數(shù)值較小的模型即可選為較好的模型。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-42【軟件操作】首先，通過acf()或pacf()函數(shù)確定模型的階數(shù)x=read.table("li15.3.txt",header=T)#從li15.3.txt中讀入樣本數(shù)據(jù)xacf(x)#作自相關圖，拖尾pacf(x)#作偏自相關圖，1階截尾2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-43【運行結果】圖15.4（a）

例15.3中數(shù)據(jù)的ACF圖【運行結果】2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-44圖15.4（b）

例15.3中數(shù)據(jù)的PACF圖【軟件操作】然后，用arima()函數(shù)對模型進行參數(shù)估計arima(x,order=c(1,0,0),method="ML",include.mean=F)#用AR(1)模型擬合參數(shù)order=控制模型的階數(shù)，按順序分別為AR的階數(shù)、差分的步數(shù)、MA的階數(shù)，參數(shù)method控制模型參數(shù)估計方法，method="ML"表示估計方法為最大似然估計，參數(shù)include.mean=F表示模型不含截距項2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-45【軟件操作】a=arima(x,order=c(1,0,0),method="ML",include.mean=F)#把估計的結果賦于變量a，以便于通過輸出對象a進行殘差分析r=a$residuals#用r來保存殘差Box.test(r,type="Ljung-Box")#檢驗模型擬合殘差是否存在序列相關,參數(shù)type控制檢驗的方法，還可以選擇"Box-Pierce"檢驗統(tǒng)計量2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-46【運行結果】Call:arima(x=x,order=c(1,0,0),include.mean=F,method="ML")Coefficients:ar1

0.6783s.e.0.0327sigma^2estimatedas0.9788:loglikelihood=-704.41,aic=1412.822023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-47【運行結果】殘差檢驗的Ljung-Box檢驗結果為：

Box-Ljungtestdata:rX-squared=0.119,df=1,p-value=0.73012023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-482023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-4915.3.2ARIMA模型的識別和估計【例7.4】(數(shù)據(jù)文件為li7.4.txt)使用例7.1中1978年至2005年國內(nèi)生產(chǎn)總值(GDP，億元)數(shù)據(jù)，2006-2008年用于做預測。使用ARIMA模型進行分析，給出2006-2008年的預測值。

2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-50【統(tǒng)計理論】在許多實際問題中，時間序列可能是非平穩(wěn)的，此時不能簡單的用ARMA模型進行擬合。一個想法是將序列進行差分，使其達到平穩(wěn)性要求，再用ARMA模型進行分析。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-51【統(tǒng)計理論】假設零均值時間序列

經(jīng)過

次差分之后可以轉化為一個平穩(wěn)時間序列

，則擬合ARMA(p,q)模型后，有ARIMA(p,d,q)中

的平穩(wěn)性條件和可逆性條件與ARMA模型是一致的。當預先確定好差分階數(shù)

，ARIMA(p,d,q)模型的參數(shù)估計與ARMA(p,q)模型是一樣的。2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-52【軟件操作】對于ARIMA模型，應先計算該序列的0(即原始數(shù)據(jù))，1，2階差分后的ACF和PACF，以確定序列的差分階數(shù)x=read.table("li15.4.txt",header=T)tx=ts(x,frequency=1,start=c(1978),end=c(2005))x1=diff(tx,lag=1)x2=diff(tx,lag=2)【軟件操作】layout(matrix(1:3,1,3))#把原序列的0，1，2階差分后的ACF條形圖平鋪在同一張圖中acf(tx,main="GDP")#作原始數(shù)據(jù)的自相關圖acf(x1,main="GDP_1")#1階差分后的自相關圖acf(x2,main="GDP_2")#2階差分后的自相關圖2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-53【軟件操作】layout(matrix(1:3,1,3))pacf(tx,main="GDP")#作原始數(shù)據(jù)的偏自相關圖pacf(x1,main="GDP_1")#1階差分后的偏自相關圖pacf(x2,main="GDP_2")#2階差分后的偏自相關圖2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-542023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-55【運行結果】圖15.5GDP數(shù)據(jù)0,1,2階差分后的ACF圖【運行結果】2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-56圖15.6GDP數(shù)據(jù)0,1,2階差分后的PACF圖【軟件操作】對于未差分數(shù)據(jù)，從其ACF和PACF的結構來看似乎符合一個AR(1)模型。但原序列不平穩(wěn)，因此可以考慮MA(1)和ARMA(1,1)。如考慮一階差分后的數(shù)據(jù)，其PACF在一階截尾，因此我們選用ARIMA(1,1,0)、ARIMA(1,1,1)和ARIMA(0,1,2)進行擬合。接著我們對模型預設各種不同的階數(shù)，并進行參數(shù)估計和Ljun-Box檢驗，通過輸出結果中的統(tǒng)計量最終確定模型階數(shù)2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-57【軟件操作】model1=arima(tx,order=c(1,0,0),method="ML",include.mean=F)model2=arima(tx,order=c(0,0,1),method="ML",include.mean=F)model3=arima(tx,order=c(1,1,0),method="ML",include.mean=F)model4=arima(tx,order=c(1,1,1),method="ML",include.mean=F)2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-58【軟件操作】model5=arima(tx,order=c(0,1,2),method="ML",include.mean=F)model6=arima(tx,order=c(1,2,0),method="ML",include.mean=F)AIC(model1,model2,model3,model4,model5,model6)Box.test(model1$residuals,type="Ljung-Box")2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-59【軟件操作】Box.test(model2$residuals,type="Ljung-Box")Box.test(model3$residuals,type="Ljung-Box")Box.test(model4$residuals,type="Ljung-Box")Box.test(model5$residuals,type="Ljung-Box")Box.test(model6$residuals,type="Ljung-Box")2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-60【運行結果】各個預設模型的擬合殘差的Ljung-Box檢驗結果匯總如下：dfAICp-valueofBox-Ljungtestmodel12600.67270.000model22678.57040.000model32504.96670.0503model43501.58560.8861model53518.94160.041model62481.32570.73272023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-61【軟件操作】計算ARIMA(1,1,1)的參數(shù)估計arima(tx,order=c(1,1,1),method="ML",include.mean=F)2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-62【運行結果】Call:arima(x=tx,order=c(1,1,1),include.mean=F,method="ML")Coefficients:ar1ma10.95430.6364s.e.0.05160.1394sigma^2estimatedas4748365:loglikelihood=-247.79,aic=501.592023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-63【軟件操作】依據(jù)殘差的平穩(wěn)性特征再次驗證模型階數(shù)的合理性：layout(matrix(1:2,1,2))acf(model4$residuals,main="reisdual")pacf(model4$residuals,main="reisdual")2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-64【運行結果】2023/3/14《統(tǒng)計學實驗》第15章時間序列分析15-65圖15.7擬合ARIMA(1,1,1)模型殘差的ACF和PACF圖【軟件操作】最后，對2006-2008年的GDP進行預測：tx.fo