ARIMA模型預(yù)測(cè)

1、1 arima模型預(yù)測(cè)一、模型選擇預(yù)測(cè)是重要的統(tǒng)計(jì)技術(shù)， 對(duì)于領(lǐng)導(dǎo)層進(jìn)行科學(xué)決策具有不可替代的支撐作用。常用的預(yù)測(cè)方法包括定性預(yù)測(cè)法、傳統(tǒng)時(shí)間序列預(yù)測(cè) （如移動(dòng)平均預(yù)測(cè)、 指數(shù)平滑預(yù)測(cè)）、現(xiàn)代時(shí)間序列預(yù)測(cè)（如arima 模型） 、灰色預(yù)測(cè)（ gm） 、線性回歸預(yù)測(cè)、非線性曲線預(yù)測(cè)、馬爾可夫預(yù)測(cè)等方法。綜合考量方法簡(jiǎn)捷性、科學(xué)性原則，我選擇arima模型預(yù)測(cè)、 gm(1,1)模型預(yù)測(cè)兩種方法進(jìn)行預(yù)測(cè)， 并將結(jié)果相互比對(duì)， 權(quán)衡取舍， 從而選擇最佳的預(yù)測(cè)結(jié)果。二 arima模型預(yù)測(cè)（一）預(yù)測(cè)軟件選擇 -r軟件arima模型預(yù)測(cè)，可實(shí)現(xiàn)的軟件較多，如spss 、sas 、eviews 、r等。使用

2、r軟件建模預(yù)測(cè)的優(yōu)點(diǎn)是：第一，r是世最強(qiáng)大、最有前景的軟件，已經(jīng)成為美國(guó)的主流。第二， r是免費(fèi)軟件。而spss 、sas 、eviews正版軟件極為昂貴，盜版存在侵權(quán)問(wèn)題，可以引起法律糾紛。第三、r軟件可以將程序保存為一個(gè)程序文件，略加修改便可用于其它數(shù)據(jù)的建模預(yù)測(cè)，便于方法的推廣。（二）指標(biāo)和數(shù)據(jù)指標(biāo)是銷(xiāo)售量 (x)，樣本區(qū)間是 1964-2013年，保存文本文件data.txt 中。（三）預(yù)測(cè)的具體步驟1、準(zhǔn)備工作（1）下載安裝 r軟件目 前 最 新 版 本 是r3.1.2， 發(fā) 布 日 期 是2014-10-31， 下 載 地 址 是/ 。

3、我使用的是 r3.1.1。（2）把數(shù)據(jù)文件 data.txt 文件復(fù)制“我的文檔”。（3）把 data.txt 文件讀入 r軟件，并起個(gè)名字。 具體操作是： 打開(kāi) r軟件，輸入（輸入每一行后，回車(chē)） ：data=read.table(data.txt,header=t) 我的文檔是默認(rèn)的工作目錄，也可以修改自定義工作目錄。2 data #查看數(shù)據(jù)回車(chē)表示執(zhí)行。完成上面操作后，r窗口會(huì)顯示：（4）把銷(xiāo)售額（ x）轉(zhuǎn)化為時(shí)間序列格式x=ts(x,start=1964) x 結(jié)果：2、對(duì) x進(jìn)行平穩(wěn)性檢驗(yàn)arma模型的一個(gè)前提條件是，要求數(shù)列是平穩(wěn)時(shí)間序列。所以，要先對(duì)數(shù)列 x 進(jìn)行平穩(wěn)性檢驗(yàn)。先做

4、時(shí)間序列圖：#后的提示語(yǔ)句是給自己看的，并不影響r 運(yùn)行3 從時(shí)間序列圖可以看出， 銷(xiāo)售量 x 不具有上升的趨勢(shì)， 也不具有起降的趨勢(shì)，初步判斷，銷(xiāo)售量 x 是平穩(wěn)時(shí)間序列。但觀察時(shí)間序列圖是不精確的，更嚴(yán)格的辦法是進(jìn)行單位根檢驗(yàn)。單位根檢驗(yàn)是通行的檢驗(yàn)數(shù)列平穩(wěn)性的工具，常用的有adf 單位根檢驗(yàn)、pp單位根檢驗(yàn)和 kpss 單位根檢驗(yàn)三種方法。單位根檢驗(yàn)的準(zhǔn)備工作是， 安裝 tseries 程序包。安裝方法： 在聯(lián)網(wǎng)狀態(tài)下，點(diǎn)菜單“ packages install packages ” ，在彈出的對(duì)話框中，選擇一個(gè)鏡像，如china(beijing1) ，確定。然后彈出附加包列表，選擇ts

5、eries，確定即可。安裝完附加包后，執(zhí)行下面操作：library(tseries) #加載 tseries 包adf.test(x) #adf檢驗(yàn)pp.test(x) #pp檢驗(yàn)kpss.test(x) #kpss 檢驗(yàn)結(jié)果：timex197019801990200020100200004000060000800004 上面分別給出了 adf檢驗(yàn)、pp檢驗(yàn)和 kpss 檢驗(yàn)的結(jié)果。其中， adf檢驗(yàn)顯示 x 是不平穩(wěn)的（ p值=0.990.05） ，而 pp檢驗(yàn)和 kpss 檢驗(yàn)則表明 x 是平穩(wěn)時(shí)間序列。再結(jié)合時(shí)間序列圖的判斷，我們認(rèn)為x是平穩(wěn)時(shí)間序列，因而符合建立arma模型的前提條件。

6、3、選擇模型做 x的自相關(guān)圖（左圖）和偏自相關(guān)圖（右圖）：acf(x) #做自相關(guān)圖pacf(x) #做偏自相關(guān)圖無(wú)論是自相關(guān)系數(shù)圖（左） ，還是偏自相關(guān)系數(shù)圖（右），都顯著第 4 階的系數(shù)突破了虛線，表明相關(guān)性顯著。因此，我們建立4 階 ar模型，寫(xiě)作 ar(4)。4、估計(jì)模型參數(shù)fit=arima(xse,order=c(4,0,0) #把估計(jì)結(jié)果取名為fit fit #查看 fit pp檢驗(yàn)的原假設(shè)是不平穩(wěn)，p 值=0.01，小于 0.05，拒絕原假設(shè)，表明序列是平穩(wěn)的。kpss 檢驗(yàn)與 pp檢驗(yàn)和 adf檢驗(yàn)不同，它的原假設(shè)是平穩(wěn)的。p值=0.1，大于 0.05，接受原假設(shè)，表明序列是

7、平穩(wěn)序列。5 上面給出了 ar 模型的回歸系數(shù)的估計(jì)值，其中，截距為44079.31，1 到 4階自回歸系數(shù)分別是0.0344，-0.0174，-0.2002和 0.4560。5、模型效果的檢驗(yàn)?zāi)Ｐ托Ч臋z驗(yàn)非常重要， 因?yàn)橹挥型ㄟ^(guò)檢驗(yàn)，才證明是可靠、有效的模型，才能進(jìn)行后續(xù)的預(yù)測(cè)分析。主要的檢驗(yàn)工具有兩個(gè)， 一是對(duì)回歸系數(shù)的顯著性檢驗(yàn)。 四個(gè)自回歸系數(shù)中，第 4 個(gè)回歸系數(shù)的 t 統(tǒng)計(jì)值=0.4560/0.1241=3.67，大于 2，因此，通過(guò)了顯著性檢驗(yàn)，表明確實(shí)存在四階自相關(guān)。 這與前面看自相關(guān)圖和偏自相關(guān)系數(shù)圖的結(jié)論相吻合。第二個(gè)檢驗(yàn)是殘差的白噪聲檢驗(yàn)（ljung-box 檢驗(yàn)） ，

8、這個(gè)最主要、最關(guān)鍵。一般來(lái)說(shuō)，只要通過(guò)了殘差的白噪聲檢驗(yàn)，則表明模型是有效的。殘差白噪聲檢驗(yàn)的r代碼：tsdiag(fit) 結(jié)果：上邊是殘差的自相關(guān)圖， 圖形顯示，除了 0 階以外， 各階自相關(guān)系數(shù)都很小，6 基本在 0 左右。表明殘差中已經(jīng)沒(méi)有多少有用的信息，殘差是純隨機(jī)序列， 即白噪聲。換個(gè)角度說(shuō)， 時(shí)間序列的有價(jià)值信息絕大部分都已經(jīng)被模型提取了，建模獲得了成功。下邊是更為精確的ljung-box檢驗(yàn)結(jié)果，所有小圈都在虛線之上（虛線值為0.05） ，表明在 0.05 的顯著性水平上，各階自相關(guān)系數(shù)和零的差別不顯著，殘差為白噪聲序列，模型效果優(yōu)良。這與上面的殘差的自相關(guān)圖相吻合。6、ari

9、ma模型預(yù)測(cè)r軟件代碼：predict(fit,n.ahead=3) #預(yù)測(cè)下三年（ 2014-2016）的數(shù)值若想預(yù)測(cè)后五年，就把3 改成 5，依此類(lèi)推。結(jié)果：pred 即 predict(預(yù)測(cè))的前四個(gè)字母，下面是時(shí)間2014-2016，表明要預(yù)測(cè)2014-2016年三年的。 結(jié)果在最后一行， 2014 年銷(xiāo)售額預(yù)測(cè)值為61768.02，2015年為 36563.83，2016 年為 45464.87。（四）模型的再檢驗(yàn) 用 aic準(zhǔn)則尋找更優(yōu)上面建模預(yù)測(cè)， 通過(guò)的顯著性檢驗(yàn)和殘差的白噪聲檢驗(yàn)，證明模型優(yōu)良， 可以進(jìn)行預(yù)測(cè)。一般的預(yù)測(cè)報(bào)告就到此結(jié)束了。但考慮到預(yù)測(cè)對(duì)于企業(yè)家的決策重要，而決

10、策的失誤將會(huì)產(chǎn)生很大的不良后果。因此，更嚴(yán)謹(jǐn)起見(jiàn)，我們建立了24 個(gè)可能的 arma模型，一個(gè)一個(gè)比較，想看一看還有沒(méi)有比前面我們建立的模型擬合效果更好的。挑選標(biāo)準(zhǔn)是國(guó)際通行的 aic準(zhǔn)則。aic是日本統(tǒng)計(jì)學(xué)家akaike 于 1973 年提出的。其基本思想是，變量越多，一般來(lái)說(shuō)模型的擬合優(yōu)度會(huì)越高。 但是我們又不能單純地以擬合的準(zhǔn)確度的衡量模型的好壞， 因?yàn)樽宰兞康脑龆鄷?huì)導(dǎo)致未知參數(shù)的增多，而參數(shù)越多， 參數(shù)估計(jì)的難度就越大， 估計(jì)的精度也越差。 因此，應(yīng)該尋求在擬合優(yōu)度和參數(shù)個(gè)數(shù)之間7 的一個(gè)平衡， aic達(dá)到最小時(shí)的模型被認(rèn)為是最優(yōu)的模型。這就是說(shuō)，我們所建立的預(yù)測(cè)模型，經(jīng)過(guò)ljung-

11、box檢驗(yàn)，表明是優(yōu)良的模型。但是，如果有好多模型通過(guò)檢驗(yàn)，證明優(yōu)良呢？這時(shí)，就可以比較aic的大小，達(dá)到優(yōu)中選優(yōu)的目的。統(tǒng)計(jì)界的經(jīng)驗(yàn)表明， arma模型最常見(jiàn)的 4 階之內(nèi)。這樣會(huì)產(chǎn)生 24 個(gè) arma模型。 分別是 ar(1) 、 ar(2)、 ar(3) 、 ar(4) 、 ma(1)、 ma(2)、 ma(3)、 ma(4)、 arma(1,1) 、arma(2,1) 、 arma3,1) 、 arma(4,1)、 arma(1,2) 、 arma(2,2)、 arma(3,2) 、 arma(4,2)、arma(1,3) 、 arma(2,3) 、 arma(3,3)、 arma(

12、4,3)、 arma(1,4)、 arma(2,4)、 arma(3,4)、arma(4,4) 。r代碼：fit=arima(x,c(1,0,0);fit #ar(1) fit=arima(x,c(2,0,0);fit #ar(2) fit=arima(x,c(3,0,0);fit #ar(3) fit=arima(x,c(4,0,0);fit #ar(4) fit=arima(x,c(0,0,1);fit #ma(1) fit=arima(x,c(0,0,2);fit #ma(2) fit=arima(x,c(0,0,3);fit #ma(3) fit=arima(x,c(0,0,4);fi

13、t #ma(4) fit=arima(x,c(1,0,1);fit #arma(1,1) fit=arima(x,c(2,0,1);fit #arma(2,1) fit=arima(x,c(3,0,1);fit #arma(3,1) fit=arima(x,c(4,0,1);fit #arma(4,1) fit=arima(x,c(1,0,2);fit #arma(1,2) fit=arima(x,c(2,0,2);fit #arma(2,2) fit=arima(x,c(3,0,2);fit #arma(3,2) fit=arima(x,c(4,0,2);fit #arma(4,2) fit

14、=arima(x,c(1,0,3);fit #arma(1,3) fit=arima(x,c(2,0,3);fit #arma(2,3) fit=arima(x,c(3,0,3);fit #arma(3,3) fit=arima(x,c(4,0,3);fit #arma(4,3) 8 fit=arima(x,c(1,0,4);fit #arma(1,4) fit=arima(x,c(2,0,4);fit #arma(2,4) fit=arima(x,c(3,0,4);fit #arma(3,4) fit=arima(x,c(4,0,4);fit #arma(4,4) 結(jié)果計(jì)算出每個(gè)模型的aic

15、值，如下：ar(4):aic = 1151.31 ar （1） ：aic = 1159.14 ar （2） ：aic = 1161.08 ar （3） ：aic = 1160.73 ma(1):aic = 1159.13 ma(2):aic = 1161.05 ma(3):aic =1160.06 ma(4):aic =1151.49 arma(1,1):aic =1155.99 arma(2,1):aic = 1157.59 arma(3,1):aic = 1154.97 arma(4,1):aic= 1152.97 arma(1,2):aic = 1156.66 arma(2,2):aic

16、 = 1156.67 arma(3,2):aic = 1155.4 arma(4,2):aic = 1154.38 arma(1,3):aic = 1155.68 arma(2,3):aic = 1157.67 arma(3,3):aic = 1156.19 arma(4,3):aic = 1153.55 arma(1,4):aic =1153.09 arma(2,4):aic = 1155.29 arma(3,4):aic = 1152.25 arma(4,4):aic = 1153.63 9 第一行是我們前面建立的四階自回歸模型ar(4) 的 aic值，為 1151.31。經(jīng)一一對(duì)比可以發(fā)

17、現(xiàn)，所有其它模型的aic值都大于 115.31。就是說(shuō)，我們前面建立的模型是所有可能模型中的最優(yōu)模型。三、gm(1,1)預(yù)測(cè)（一）方法簡(jiǎn)介灰色系統(tǒng)理論是由華中理工大學(xué)鄧聚龍教授于1982 年提出并加以發(fā)展的。二十幾年來(lái)，引起了不少?lài)?guó)內(nèi)外學(xué)者的關(guān)注，得到了長(zhǎng)足的發(fā)展。目前，在我國(guó)已經(jīng)成為社會(huì)、經(jīng)濟(jì)、科學(xué)技術(shù)在等諸多領(lǐng)域進(jìn)行預(yù)測(cè)、 決策、評(píng)估、規(guī)劃控制、系統(tǒng)分析與建模的重要方法之一。而其精華在于 gm(1,1)模型，gm(1,1)被廣泛應(yīng)用于預(yù)測(cè)，并且預(yù)測(cè)效果很好，其使用限制條件是：原始數(shù)據(jù)單調(diào)，預(yù)測(cè)背景呈現(xiàn)穩(wěn)定發(fā)展趨勢(shì)；其優(yōu)勢(shì)是：適用于觀測(cè)數(shù)據(jù)較少的預(yù)測(cè)問(wèn)題，算法簡(jiǎn)單易行，預(yù)測(cè)精度相對(duì)較高。鑒于

18、 gm(1,1)預(yù)測(cè)適合于數(shù)據(jù)較少情況， 我們選擇 2008-2013年的近期數(shù)據(jù)，進(jìn)行預(yù)測(cè)嘗試。（二） r軟件操作1、把腳本文件 gm11.r和 data2.txt 復(fù)制到“我的文檔”；2、r窗口輸入命令：source(gm11.r) data2=read.table(data2.txt,header=t) x=data2$x x gm11(x,length(x)+1) #預(yù)測(cè)后 1 期值預(yù)測(cè)結(jié)果： 2014 年銷(xiāo)售量分別為28385.3。在 r顯示結(jié)果中，是 x(0)的模擬值的最后 1 個(gè)數(shù)。10 (三)模型檢驗(yàn)利用 r軟件進(jìn)行操作處理結(jié)果顯示如下：平均相對(duì)誤差 = 22.11903% 相

19、對(duì)精度 = 77.88097% 可見(jiàn)， gm(1,1)模型預(yù)測(cè)的誤差較大，相對(duì)精度不高。再進(jìn)行后驗(yàn)差比值檢驗(yàn)。原始序列的方差：215970035633s =殘差序列 e 的方差：22192438750s =后驗(yàn)差比值為：c值=s2/s1=0.7108341 按照灰色預(yù)測(cè)理論， c值大于 0.65，表明 gm(1,1)預(yù)測(cè)精度等級(jí)為：不合格（標(biāo)準(zhǔn)見(jiàn)下表）。說(shuō)明該模型的預(yù)測(cè)結(jié)果不太可靠。123456730000400005000060000時(shí)間序列值預(yù)測(cè)序列原始序列11 小誤差概率檢驗(yàn)：r軟件給出了小誤差概率值，p= 4 / 6 = 0.6666667。屬于(0.7,0.8，gm(1,1)模型預(yù)測(cè)

20、精度等級(jí)為：勉強(qiáng)合格（標(biāo)準(zhǔn)見(jiàn)上表）。綜合來(lái)看，對(duì)于我公司的數(shù)據(jù)，gm(1,1)模型效果并不理想。四簡(jiǎn)要結(jié)論前面用 arima模型和 gm(1,1)模型，對(duì)我公司的銷(xiāo)售量進(jìn)行了預(yù)測(cè)。檢驗(yàn)表明，gm(1,1)模型的精度較低，而arima 模型則通過(guò) ljung-box 檢驗(yàn)，表明模型優(yōu)良。而且，我們進(jìn)行了全方位的24 個(gè)模型的一一比對(duì)，發(fā)現(xiàn)我們建立的模型在所有可能模型中最優(yōu)。因而，建立重點(diǎn)參照arima模型的預(yù)測(cè)結(jié)果。作為統(tǒng)計(jì)人員， 基于對(duì)統(tǒng)計(jì)工作的經(jīng)驗(yàn)認(rèn)識(shí)，明白任何預(yù)測(cè)， 都是根據(jù)過(guò)往的歷史信息進(jìn)行的。如果歷史條件不變，預(yù)測(cè)結(jié)果可靠。一旦歷史條件改變，則預(yù)測(cè)結(jié)果不一定準(zhǔn)確。而且，只要是預(yù)測(cè)，無(wú)論

21、方法多么科學(xué)，分析多么嚴(yán)謹(jǐn)，都沒(méi)有百分百的準(zhǔn)確，都會(huì)有預(yù)測(cè)誤差。五進(jìn)一步建議單變量時(shí)間預(yù)測(cè)， 是依據(jù)單一變量 銷(xiāo)售量的歷史信息進(jìn)行預(yù)測(cè)， 并沒(méi)有考慮其它因素。結(jié)果僅具有參考意義。建議領(lǐng)導(dǎo)層根據(jù)自己的經(jīng)驗(yàn)、 銷(xiāo)售部門(mén)根據(jù)其經(jīng)驗(yàn)， 尤其是對(duì)未來(lái)經(jīng)濟(jì)形勢(shì)的變化特點(diǎn)，修正前面的預(yù)測(cè)結(jié)果，使預(yù)測(cè)結(jié)果更加完善。附錄： gm（1,1）基本理論g 表示 grey(灰)， m 表示 model (模型 )，前一個(gè)“ 1”表示一階，后一個(gè)“1”表示一個(gè)變量， gm(1, 1)則是一階，一個(gè)變量的微分方程模型。給定等時(shí)間間隔的數(shù)據(jù)列，且設(shè)數(shù)據(jù)列單調(diào)：12 ),n(), 2(), 1 ()(,n21xxxkxk，k表

22、示時(shí)刻，kxkx)(表示 tk時(shí)刻某量的觀測(cè)值， 不妨設(shè)1kkxx，1, 2, 1nk，將數(shù)據(jù)列記成：0n030201)0(,xxxxx)0(x表示原始數(shù)據(jù)序列。比如：697.3,390.3,337. 3,278.3,874. 2)0(x。對(duì)原始數(shù)據(jù)作一次累加生成：即令nkxxkiik,2, 11)0()1(得一次累加生成數(shù)序列為：)1()1(21(1)1(,nxxxx）在此，)1(kx=2.874, 6.152, 9.489, 12.879, 16.558 給定的原始數(shù)據(jù)序列)0(kx已經(jīng)是單增序列，經(jīng)一次累加后生成的累加數(shù)序列具有更強(qiáng)烈的單調(diào)性。我們知道指數(shù)序列是單調(diào)的，但是，單調(diào)序列卻不

23、一定是指數(shù)型的，不過(guò)強(qiáng)烈的單調(diào)序列可近似看做是指數(shù)的，即可用指數(shù)型曲線進(jìn)行彌合。如果用指數(shù)曲線來(lái)彌合一次累加生成序列，那么，這條指數(shù)曲線一定是某個(gè)一階線性常系數(shù)微分方程uaxdtdx)1()1(6) 的滿(mǎn)足某個(gè)初始條件的一條積分曲線：aueauxxat)1()1()1(即aueauxxat)1()0()1(（7）其中 a， u 是待確定的未知參數(shù)，該微分方程中的導(dǎo)數(shù)dtdx)1(可用差商近似表示。tkxtkxdtdxt)1()1(0)1(limt為時(shí)間間隔，將時(shí)間間隔t看做是單位時(shí)間間隔，并且認(rèn)為時(shí)間被充分細(xì)化(秒，毫秒。微秒, 事實(shí)上只要單位時(shí)間內(nèi)函數(shù)的增量相對(duì)很小，這個(gè)單位時(shí)間間隔也可以是日，13 月，年等。 ) 此時(shí)有kxkxdtdx)1()1()1(1注意到一次累加生成數(shù))()1(tx在時(shí)刻 tk1 與 tk 時(shí)的差為：11)0()1()1(kxkxkx而)1(dtdx是在 k, k1上某一點(diǎn)取值，既然是近似，索性將dtdx)1(的值取在點(diǎn)k+1，即11)0()1()1(1)1(kxkxkxdtdxkt于是，一階線性常系數(shù)微分方程uaxdtdx)1()1(可近似化為：11)1()0(ktkut