應(yīng)用時(shí)間序列分析實(shí)驗(yàn)報(bào)告

1、 應(yīng)用時(shí)間序列分析實(shí)驗(yàn)報(bào)告學(xué)院名稱 理學(xué)院 專業(yè)班級 應(yīng)用統(tǒng)計(jì)學(xué)14-2 學(xué)生姓名 張艷雪 學(xué) 號 201411081051 齊魯工業(yè)大學(xué)實(shí)驗(yàn)報(bào)告 成績 課程名稱 應(yīng)用時(shí)間序列分析實(shí)驗(yàn) 指導(dǎo)教師 黃玉林 實(shí)驗(yàn)日期 2017.6.30 院（系） 理學(xué)院 專業(yè)班級 統(tǒng)計(jì)14-2 實(shí)驗(yàn)地點(diǎn) 機(jī)電樓C428 學(xué)生姓名 張艷雪 學(xué)號 201411081051 同組人 無 實(shí)驗(yàn)項(xiàng)目名稱 ARIMA模型、確定性分析法,多元時(shí)間序列建模 一、 實(shí)驗(yàn)?zāi)康暮鸵?.熟悉非平穩(wěn)序列的確定性分析法：趨勢分析、季節(jié)效應(yīng)分析、綜合分析2.熟悉差分平穩(wěn)序列的建模步驟。3.掌握單位根檢驗(yàn)、協(xié)整檢驗(yàn)、動(dòng)態(tài)回歸模型的建立。二、

2、實(shí)驗(yàn)原理1. 序列的各種變化都?xì)w結(jié)于四大因素的綜合影響：長期趨勢（Trend）,循環(huán)波動(dòng)（Circle）,季節(jié)性變化（Season）,機(jī)波動(dòng)（Immediate）.常假設(shè)它們有如下的相互模型：加法模型 乘法模型 混合模型 模型結(jié)構(gòu)不唯一 2.非平穩(wěn)序列如果能通過適當(dāng)階數(shù)的差分后實(shí)現(xiàn)平穩(wěn)，就可以對差分后序列進(jìn)行ARMA模型擬合了，所以ARIMA模型是差分運(yùn)算與ARMA模型的組合 3.單位根檢驗(yàn)：（1）DF檢驗(yàn)；（2）ADF檢驗(yàn)； (3)PP檢驗(yàn)； 4.動(dòng)態(tài)回歸模型ARIMAX 如果兩個(gè)非平穩(wěn)序列之間具有協(xié)整關(guān)系，則先建立它們的回歸模型，再對平穩(wěn)的殘差序列建立ARMA模型。三、實(shí)驗(yàn)內(nèi)容1、P202

3、頁：第7 題（X11因素分解法）2、P155頁：第3題（乘積季節(jié)模型）3、P240頁：第4題出口為，進(jìn)口為，回答以下問題（1）畫出，的時(shí)序圖，用單位根檢驗(yàn)序列它們的平穩(wěn)性；（2）對分別擬合模型（提示：建立ARIMA模型）；(3)考察的協(xié)整關(guān)系，建立的協(xié)整模型，同時(shí)建立誤差修正模型。四、實(shí)驗(yàn)過程(一) P202頁：第7 題（X11因素分解法）1.繪制序列時(shí)序圖。（程序見附錄）由上圖可得季節(jié)序列的振幅隨序列水平的變化而變化，所以季節(jié)效應(yīng)與趨勢效應(yīng)不獨(dú)立，采用乘法模型: xt=Tt×St×It2.進(jìn)入x-11季節(jié)調(diào)整模型經(jīng)過三個(gè)階段共十步的重復(fù)迭代后，得到如下的擬合效果圖：顯然，

4、該地區(qū)奶牛的月度產(chǎn)奶量序列具有顯著的季節(jié)變動(dòng)特征。(二) P155頁：第3題（乘積季節(jié)模型） 1.繪制序列時(shí)序圖。 繪制時(shí)序圖，如圖1所示（程序見附錄1）。圖1 美國月度事故死亡人數(shù)序列時(shí)序圖時(shí)序圖顯示該序列具有以年為周期的季節(jié)效應(yīng)。2.差分平穩(wěn)化：對原序列作1階12步差分，希望提取原序列季節(jié)效應(yīng)，差分后序列時(shí)序圖如圖2所示。圖2 美國月度事故死亡人數(shù)1階12步差分后序列時(shí)序圖時(shí)序圖顯示差分后序列類似平穩(wěn)。3.模型定階：考察差分后序列自相關(guān)圖，如圖3，進(jìn)一步確定平穩(wěn)性判斷，并估計(jì)擬合模型的階數(shù)。圖3 美國月度事故死亡人數(shù)1階12步差分后序列自相關(guān)圖 自相關(guān)圖顯示延遲12階自相關(guān)系數(shù)顯著大于2倍

5、標(biāo)準(zhǔn)差范圍，這說明差分后序列中仍蘊(yùn)含著非常顯著的季節(jié)效應(yīng)。延遲1階的自相關(guān)系數(shù)也大于2倍的標(biāo)準(zhǔn)差，這說明差分后序列還具有短期相關(guān)性。觀察偏自相關(guān)圖，如圖4，得到的結(jié)論和上面的結(jié)論一致。圖4 美國月度事故死亡人數(shù)1階12步差分后序列偏自相關(guān)圖圖5 序列白噪聲檢驗(yàn)圖5顯示，原序列延遲各階LB統(tǒng)計(jì)量的P值小于顯著性水平0.05，所以拒絕原假設(shè)，序列不通過白噪聲檢驗(yàn)。 根據(jù)差分后序列的自相關(guān)圖和偏自相關(guān)圖的性質(zhì)，擬合乘積季節(jié)模型。 自相關(guān)圖顯示，12階以內(nèi)的自相關(guān)系數(shù)1階截尾，偏自相關(guān)圖顯示，12階以內(nèi)的偏自相關(guān)系數(shù)1階截尾，所以嘗試使用ARMA(1，0)模型提取差分后序列的短期自相關(guān)信息。再考慮季節(jié)

6、自相關(guān)特征，這時(shí)考察延遲12階、24階等以周期長度為單位的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)的特征。自相關(guān)圖顯示延遲12階自相關(guān)系數(shù)顯著非零，而偏自相關(guān)圖顯示延遲12階偏自相關(guān)系數(shù)顯著非零，這時(shí)用以12步為周期的模型提取差分后序列的季節(jié)自相關(guān)信息。4. 參數(shù)估計(jì)：圖6 擬合模型綜合前面的差分信息，我們要擬合的乘積季節(jié)模型為。使用條件最小二乘估計(jì)方法，確定該模型的口徑為：5.模型檢驗(yàn)：對序列擬合模型，模型及模型參數(shù)的顯著性檢驗(yàn)如圖7、8所示。圖7 模型參數(shù)的顯著性 由圖7知，擬合效果顯示模型參數(shù)顯著。圖8 殘差白噪聲檢驗(yàn) 對擬合模型進(jìn)行白噪聲檢驗(yàn)，結(jié)果顯示P值都大于顯著性水平0.05.接受原假設(shè)，殘差序列

7、通過白噪聲檢驗(yàn)，模型顯著，說明模型擬合良好，對序列相關(guān)信息提取充分。將序列擬合值和序列觀察值聯(lián)合作圖，如圖9所示。圖9 美國月度事故死亡人數(shù)擬合效果圖說明：圖中，點(diǎn)為序列觀察值；曲線為序列擬合值。從圖9可以直觀地看出該乘積季節(jié)模型對原序列的擬合效果良好。（三）P240頁：第4題1畫出，的時(shí)序圖，用單位根檢驗(yàn)序列的平穩(wěn)性； 輸出時(shí)序圖如圖1所示（程序見附錄2）。圖1 我國出口總額Xt、進(jìn)口總額yt時(shí)序圖圖1中，黑色為出口總額xt序列時(shí)序圖，紅色為進(jìn)口總額yt序列時(shí)序圖。從圖1中可以看出出口總額xt序列、進(jìn)口總額yt序列均顯著非平穩(wěn)，這個(gè)直觀判斷還可以通過單位根檢驗(yàn)驗(yàn)證。同時(shí)時(shí)序圖顯示這兩個(gè)序列具

8、有某種同變關(guān)系。對我國出口總額序列xt進(jìn)行ADF檢驗(yàn)，單位根檢驗(yàn)結(jié)果如圖2所示。圖2 出口總額xt白噪聲、單位根檢驗(yàn)檢驗(yàn)結(jié)果顯示，無論考慮何種類型的模型，檢驗(yàn)統(tǒng)計(jì)量的P值均顯著大于0.05的顯著性水平，所以可以認(rèn)為中國我國出口總額序列xt顯著非平穩(wěn)，且這六種處理均不能實(shí)現(xiàn)殘差序列平穩(wěn)。 對我國進(jìn)口總額序列yt進(jìn)行ADF檢驗(yàn)，單位根檢驗(yàn)結(jié)果如圖3所示。圖3 進(jìn)口總額yt白噪聲、單位根檢驗(yàn)同出口序列xt的檢驗(yàn)結(jié)果一樣，在顯著性水平取為0.05時(shí)，可以認(rèn)為我國進(jìn)口序列yt非平穩(wěn)，且這六種處理均不能實(shí)現(xiàn)殘差序列平穩(wěn)。顯然，這兩個(gè)序列的ADF檢驗(yàn)結(jié)果與根據(jù)時(shí)序圖得到的直觀判斷完全一致2對分別擬合模型（提

9、示：建立ARIMA模型）； 對我國出口對數(shù)序列l(wèi)nxt和進(jìn)口對數(shù)序列l(wèi)nyt繪制時(shí)序圖，如圖4所示。圖4 我國出口總額Xt、進(jìn)口總額yt取對數(shù)時(shí)序圖 圖4中，黑色線代表我國出口對數(shù)序列l(wèi)nxt，紅色線代表我國進(jìn)口對數(shù)序列l(wèi)nyt。時(shí)序圖顯示這兩個(gè)對數(shù)序列有顯著的上升趨勢，為典型的非平穩(wěn)序列。同時(shí)時(shí)序圖顯示這兩個(gè)序列具有某種同變關(guān)系。因?yàn)樾蛄谐尸F(xiàn)出近似線性趨勢，所以選擇1階差分。1階差分后出口對數(shù)序列l(wèi)nxt時(shí)序圖如圖5所示。圖5 對數(shù)序列Lnx差分時(shí)序圖 時(shí)序圖顯示，lnxt差分后序列在均值附近比較穩(wěn)定地波動(dòng)。為了進(jìn)一步確定平穩(wěn)性，考察差分后序列的自相關(guān)圖，如圖6所示。圖6 對數(shù)序列Lnxt差

10、分后自相關(guān)圖自相關(guān)圖顯示序列有很強(qiáng)的短期相關(guān)性，所以可以初步認(rèn)為lnxt1階差分后序列平穩(wěn)。對平穩(wěn)的1階差分序列進(jìn)行白噪聲檢驗(yàn)，白噪聲檢驗(yàn)結(jié)果如圖7所示。圖7 lnxt一階差分后序列白噪聲檢驗(yàn)在檢驗(yàn)的顯著性水平取為0.05的條件下，由于延遲6階、12階的P值均小于0.05，所以lnxt差分后的序列不能視為白噪聲序列，即差分后序列還蘊(yùn)含著不容忽視的相關(guān)信息可以提取。對平穩(wěn)非白噪聲差分序列擬合ARMA模型，1階差分后序列的自相關(guān)圖(見圖6）已經(jīng)顯示該序列有不截尾的性質(zhì)。再考察其偏自相關(guān)系數(shù)的性質(zhì)，如圖8所示。圖8 對數(shù)序列Lnxt差分后偏自相關(guān)圖偏自相關(guān)圖顯示出1階截尾性，所以考慮用AR(1)模型

11、擬合lnxt1階差分后序列?？紤]到前面已經(jīng)進(jìn)行的1階差分運(yùn)算，實(shí)際上是用模型擬合原序列。對序列擬合模型，模型參數(shù)及模型的顯著性檢驗(yàn)如圖9、10所示。圖9 模型參數(shù)顯著性檢驗(yàn) 由圖9知，系數(shù)顯著性檢驗(yàn)顯示兩參數(shù)均顯著。 對殘差序列進(jìn)行白噪聲檢驗(yàn)，檢驗(yàn)結(jié)果如圖10所示。圖10 殘差白噪聲檢驗(yàn) 顯然，擬合檢驗(yàn)統(tǒng)計(jì)量的P值都顯著大于顯著性檢驗(yàn)水平0.05，可以認(rèn)為殘差序列即為白噪聲序列，模型顯著，這說明模型對lnxt序列建模成功。圖11 模型在條件最小二乘估計(jì)原理下，擬合結(jié)果為：2將對數(shù)序列擬合值lnxt和對數(shù)序列觀察值lnxt聯(lián)合作圖，如圖12所示。圖12 對數(shù)序列Lnxt擬合效果圖說明：圖中，星號

12、為序列觀察值；曲線為擬合值。從圖可以直觀地看出該模型對原序列的擬合效果良好。因?yàn)閷?shù)序列l(wèi)nyt呈現(xiàn)出近似線性趨勢，所以選擇1階差分。1階差分后進(jìn)口對數(shù)序列l(wèi)nyt時(shí)序圖如圖13所示。圖13 對數(shù)序列Lny差分時(shí)序圖 時(shí)序圖顯示，lnyt差分后序列在均值附近比較穩(wěn)定地波動(dòng)。為了進(jìn)一步確定平穩(wěn)性，考察差分后序列的自相關(guān)圖，如圖14所示。圖14 對數(shù)序列Lnyt差分后自相關(guān)圖自相關(guān)圖顯示序列有很強(qiáng)的短期相關(guān)性，所以可以初步認(rèn)為lnyt1階差分后序列平穩(wěn)。對平穩(wěn)的1階差分序列進(jìn)行白噪聲檢驗(yàn)，白噪聲檢驗(yàn)結(jié)果如圖15所示。圖15 lnyt一階差分后序列白噪聲檢驗(yàn)在檢驗(yàn)的顯著性水平取為0.05的條件下，由

13、于延遲6階的P值小于0.05，所以lnyt差分后的序列不能視為白噪聲序列，即差分后序列還蘊(yùn)含著不容忽視的相關(guān)信息可以提取。對平穩(wěn)非白噪聲差分序列擬合ARMA模型，1階差分后序列的自相關(guān)圖(見圖14）已經(jīng)顯示該序列有1階截尾的性質(zhì)。再考察其偏自相關(guān)系數(shù)的性質(zhì)，如圖16所示。圖16 對數(shù)序列Lnyt差分后偏自相關(guān)圖偏自相關(guān)圖顯示該序列1階截尾的性質(zhì)，所以考慮用AR(1)模型擬合lnyt1階差分后序列?？紤]到前面已經(jīng)進(jìn)行的1階差分運(yùn)算，實(shí)際上是用模型擬合原序列。對序列擬合模型，模型參數(shù)及模型的顯著性檢驗(yàn)如圖17、18所示。圖17 模型參數(shù)顯著性檢驗(yàn) 由圖17知，系數(shù)顯著性檢驗(yàn)顯示兩參數(shù)均顯著。 對殘

14、差序列進(jìn)行白噪聲檢驗(yàn)，檢驗(yàn)結(jié)果如圖18所示。圖18 殘差白噪聲檢驗(yàn) 顯然，擬合檢驗(yàn)統(tǒng)計(jì)量的P值都顯著大于顯著性檢驗(yàn)水平0.05，可以認(rèn)為殘差序列即為白噪聲序列，模型顯著。這說明模型對該序列建模成功。圖19 模型在條件最小二乘估計(jì)原理下，擬合結(jié)果為：將對數(shù)序列擬合值lnyt和對數(shù)序列觀察值lnyt聯(lián)合作圖，如圖20所示。圖20 對數(shù)序列Lnyt擬合效果圖說明：圖中，星號為序列觀察值；曲線為擬合值。從圖20可以直觀地看出該模型對原序列的擬合效果良好。3考察的協(xié)整關(guān)系，建立的協(xié)整模型，同時(shí)建立誤差修正模型。 對我國出口對數(shù)序列l(wèi)nxt和進(jìn)口對數(shù)序lnyt繪制時(shí)序圖，如圖4所示?？梢园l(fā)現(xiàn)時(shí)序圖顯示這兩

15、個(gè)序列具有某種同變關(guān)系，可以考慮建立ARIMAX模型。對lnxt、lnyt、lnxt1階差分、序列分別進(jìn)行單位根檢驗(yàn)（ADF）。輸出結(jié)果如圖2124所示。圖21 對數(shù)序列l(wèi)nXt1階單位根檢驗(yàn)檢驗(yàn)結(jié)果顯示，無論考慮何種類型的模型，檢驗(yàn)統(tǒng)計(jì)量的P值均顯著大于0.05的顯著性水平，所以可以認(rèn)為中國我國出口總額對數(shù)序列l(wèi)nxt顯著非平穩(wěn)，且這六種處理均不能實(shí)現(xiàn)殘差序列平穩(wěn)。圖22 對數(shù)序列l(wèi)nyt1階單位根檢驗(yàn)同出口對數(shù)序列l(wèi)nxt的檢驗(yàn)結(jié)果一樣，在顯著性水平取為0.05時(shí)，可以認(rèn)為我國進(jìn)口對數(shù)序列l(wèi)nyt非平穩(wěn)，且這六種處理均不能實(shí)現(xiàn)殘差序列平穩(wěn)。圖23 1階單位根檢驗(yàn) 檢驗(yàn)結(jié)果顯示，無論考慮何種

16、類型的模型，檢驗(yàn)統(tǒng)計(jì)量的P值均顯著小于0.05的顯著性水平，拒絕原假設(shè)，所以可以認(rèn)為中國我國出口總額對序列顯著平穩(wěn)。圖24 1階單位根檢驗(yàn)同出口序列的檢驗(yàn)結(jié)果一樣，在顯著性水平取為0.05時(shí)，可以認(rèn)為我國進(jìn)口序列平穩(wěn)，且這六種處理均能實(shí)現(xiàn)殘差序列平穩(wěn)。利用最小二乘估計(jì)，回歸模型輸出結(jié)果如圖25所示。圖25 回歸模型結(jié)果構(gòu)造出的回歸模型如下：圖26 對數(shù)序列Lnyt與對數(shù)序列l(wèi)nxt之間的相關(guān)圖相關(guān)圖顯示對數(shù)序列l(wèi)nyt在延遲階數(shù)為零時(shí)與對數(shù)序列l(wèi)nxt相關(guān)關(guān)系最大。因此可以將對數(shù)序列l(wèi)nyt與對數(shù)序列l(wèi)nxt同期建模。圖27 殘差單位根檢驗(yàn) 殘差序列平穩(wěn)，說明對數(shù)序列l(wèi)nyt與對數(shù)序列l(wèi)nxt

17、之間具有協(xié)整的關(guān)系，我們可以大膽的在這兩個(gè)對數(shù)序列之間建立動(dòng)態(tài)回歸模型而不必?fù)?dān)心虛假回歸問題。圖28 模型參數(shù)顯著性檢驗(yàn)，無常數(shù)項(xiàng) 考察殘差序列白噪聲檢驗(yàn)結(jié)果，如圖29所示。圖29 殘差序列白噪聲檢驗(yàn)輸出結(jié)果顯示，延遲各階LB統(tǒng)計(jì)量的P值都大于顯著性水平0.05，可以認(rèn)為殘差序列為白噪聲序列，對序列相關(guān)信息提取充分。圖30 模型 根據(jù)輸出的模型擬合結(jié)果可知，最后的擬合模型口徑為： 將對數(shù)序列擬合值lnyt和對數(shù)序列觀察值lnyt聯(lián)合作圖，如圖31所示。圖31 對數(shù)序列Lnyt擬合效果圖說明：圖中，星號為序列觀察值；曲線為擬合值。從圖可以直觀地看出該ARIMAX模型對原序列的擬合效果良好。將序列

18、擬合值yt和序列觀察值yt聯(lián)合作圖，如圖32所示。圖32 yt擬合效果圖說明：圖中，星號為序列觀察值；曲線為擬合值。從圖可以直觀地看出該ARIMAX模型對原序列的擬合效果良好。構(gòu)造誤差修正模型: 為了研究我國的進(jìn)出口總額的短期波動(dòng)特征，我們利用差分序列和以及前期誤差序列，構(gòu)造ECM模型：輸出結(jié)果如圖33所。圖33誤差修正模型顯著性檢驗(yàn)由圖33輸出結(jié)果結(jié)果知ECM模型為：附錄程序1：data example4_7; input x; t=intnx('quarter','1jan1978'd,_n_-1); format t yyq4.; cards; 58956

19、164065672769764059956857755358260056665367374271666061758358756559862861868870577073667863960461159463465862270972278275670265361562160263567763573675581179873569766166764568871366776278483781776772268168766069871769677579685882678374070170667771173469078580587184580176472572369073475070780782488685

20、9819783740747711751 ; proc x11 data=example4_7; quarterly date=t; var x; output out=out b1=x d10=season d11=adjusted d12=trend d13=irr; data out;set out;estimate=trend*season/100; proc gplot data=out; plot season*t=2 adjusted*t=2 trend*t=2 irr*t=2; plot x*t=1 estimate*t=2/overlay; symbol1 c=black i=

21、join v=star; symbol2 c=red i=join v=none w=2; run;程序2：data ti5_3;input x;dif1_12=dif12(dif(x);time=intnx('month','1jan1973'd,_n_-1);format time year4.;cards;9007.008106.008928.009137.0010017.0010826.0011317.0010744.009713.009938.009161.008927.007750.006981.008038.008422.008714.009512

22、.0010120.009823.008743.009129.008710.008680.008162.007306.008124.007870.009387.009556.0010093.009620.008285.008433.008160.008034.007717.007461.007776.007925.008634.008945.0010078.009179.008037.008488.007874.008647.007792.006957.007726.008106.008890.009299.0010625.009302.008314.008850.008265.008796.0

23、07836.006892.007791.008129.009115.009434.0010484.009827.009110.009070.008633.009240.00;run;proc gplot;plot x*time=1 dif1_12*time=2;symbol1 c=coral v=circle i=join;symbol2 c=blue v=star i=join;run;proc arima;identify var=x(1,12);estimate p=1 q=(1)(12) ;forecast lead=0 id=time out=out;run;proc gplot d

24、ata=out;plot x*time=1 forecast*time=2 /overlay;symbol1 c=black i=none v=dot h=0.2;symbol2 c=red i=join v=none;run;程序3:data ti6_4; input year xt yt;lnxt=log(xt);lnyt=log(yt); diflnx=dif(lnxt);diflny=dif(lnyt);cards;19502021.3195124.235.3195227.137.5195334.846.119544044.7195548.761.1195655.753195754.5

25、5019586761.7195978.171.2196063.365.1196147.743196247.133.819635035.7196455.442.1196563.155.319666661.1196758.853.4196857.650.9196959.847.2197056.856.1197168.552.4197282.9641973116.9103.61974139.4152.81975143147.41976134.8129.31977139.7132.81978167.6187.41979211.7242.91980271.2298.81981367.6367.71982

26、413.8357.51983438.3421.81984580.5620.51985808.91257.819861082.11498.3198714701614.219881766.72055.1198919562199.919902985.82574.319913827.13398.719924676.34443.319935284.85986.2199410421.89960.1199512451.811048.1199612576.411557.4199715160.711806.5199815223.611626.1199916159.813736.5200020634.418638

27、.8200122024.420159.2200226947.924430.3200336287.934195.6200449103.346435.8200562648.154273.7200677594.663376.9200793455.673284.62008100394.979526.5;run;proc gplot data=ti6_4;plot xt*year=1 yt*year=2/overlay;symbol1 c=black i=join v=none;symbol2 c=red i=join v=none w=2 l=2;run;proc arima data=ti6_4;i

28、dentify var=xt stationarity=(adf=1); identify var=yt stationarity=(adf=1);run;proc gplot data=ti6_4;plot lnxt*year=1 lnyt*year=2/overlay;plot diflnx*year=1 diflny*year=2;symbol1 c=black i=join v=circle;symbol2 c=red i=join v=star;run;proc arima;identify var=lnxt(1) ;estimate p=1;forecast lead=0 id=year out=out1;identify var=lnyt(1);estimate p=1;forecast lead=0 id=year out=out2;run;proc gplot data=out1 ;plot lnxt*year=1 forecast*year=2 /overlay;symbol1 c=black i=none v=star;symbol2 c=red i=join v=none;run;proc gplot data