時(shí)間序列建模培訓(xùn)

1、時(shí)間序列建模培訓(xùn)n背景及引例n平穩(wěn)序列n平穩(wěn)序列建模（結(jié)合例子來(lái)講）n非平穩(wěn)序列處理方法n時(shí)間序列建模的Matlab實(shí)現(xiàn)（例子）n自主練習(xí)概率統(tǒng)計(jì)學(xué)科中應(yīng)用性較強(qiáng)的一個(gè)分支概率統(tǒng)計(jì)學(xué)科中應(yīng)用性較強(qiáng)的一個(gè)分支廣泛的應(yīng)用領(lǐng)域：廣泛的應(yīng)用領(lǐng)域：金融經(jīng)濟(jì)金融經(jīng)濟(jì)氣象水文氣象水文信號(hào)處理信號(hào)處理機(jī)械振動(dòng)機(jī)械振動(dòng)一 背景及引例時(shí)時(shí)間序列的定義間序列的定義 n時(shí)間序列時(shí)間序列:按時(shí)間次序排列的隨機(jī)變量序列按時(shí)間次序排列的隨機(jī)變量序列n 個(gè)觀測(cè)樣本個(gè)觀測(cè)樣本:隨機(jī)序列的隨機(jī)序列的 個(gè)有序觀個(gè)有序觀測(cè)測(cè)值值 n稱(chēng)序列稱(chēng)序列 是時(shí)間序列是時(shí)間序列(1.1)(1.1)的一次實(shí)現(xiàn)或一條軌道的一次實(shí)現(xiàn)或一條軌道)1 .

2、1 (,21XX)2 . 1 (,21nxxxnn) 3 . 1 (,21xxn按照時(shí)間的順序把隨機(jī)事件變化發(fā)展的過(guò)按照時(shí)間的順序把隨機(jī)事件變化發(fā)展的過(guò)程記錄下來(lái)就構(gòu)成了一個(gè)時(shí)間序列。對(duì)時(shí)程記錄下來(lái)就構(gòu)成了一個(gè)時(shí)間序列。對(duì)時(shí)間序列進(jìn)行觀察、研究，找尋它變化發(fā)展間序列進(jìn)行觀察、研究，找尋它變化發(fā)展的規(guī)律，預(yù)測(cè)它將來(lái)的走勢(shì)就是時(shí)間序列的規(guī)律，預(yù)測(cè)它將來(lái)的走勢(shì)就是時(shí)間序列分析。分析。太陽(yáng)黑子對(duì)地球的影響n會(huì)出現(xiàn)磁暴現(xiàn)象會(huì)出現(xiàn)磁暴現(xiàn)象 n會(huì)引起地球上氣候的變化會(huì)引起地球上氣候的變化 n會(huì)影響地球上的地震會(huì)影響地球上的地震n會(huì)影響樹(shù)木生長(zhǎng)會(huì)影響樹(shù)木生長(zhǎng)n會(huì)影響到我們的身體會(huì)影響到我們的身體n Wolfe

3、r記錄的300年的太陽(yáng)黑子數(shù)房地產(chǎn)業(yè)、房?jī)r(jià)房地產(chǎn)業(yè)、房?jī)r(jià)n關(guān)乎國(guó)計(jì)民生的支柱產(chǎn)業(yè)關(guān)乎國(guó)計(jì)民生的支柱產(chǎn)業(yè)n影響著城鎮(zhèn)居民的住房消費(fèi)影響著城鎮(zhèn)居民的住房消費(fèi)n影響著水泥，鋼鐵，建材，冶金等相關(guān)影響著水泥，鋼鐵，建材，冶金等相關(guān)行業(yè)的發(fā)展行業(yè)的發(fā)展n影響著地方政府財(cái)政收入影響著地方政府財(cái)政收入n.杭州近三年房?jī)r(jià)走勢(shì)n股市是經(jīng)濟(jì)的晴雨表股市是經(jīng)濟(jì)的晴雨表 n從股市本身看，我國(guó)股市的確有自己的從股市本身看，我國(guó)股市的確有自己的特點(diǎn)特點(diǎn) n股票是一種高風(fēng)險(xiǎn)的資本投資股票是一種高風(fēng)險(xiǎn)的資本投資n1985至2000年廣州月平均氣溫國(guó)際航空公司月旅客數(shù)050100150100200300400500600700

4、化學(xué)反應(yīng)過(guò)程中溶液濃度數(shù)據(jù)化學(xué)反應(yīng)過(guò)程中溶液濃度數(shù)據(jù)0204060801001201401601802001616.51717.51818.5二、平穩(wěn)序列自協(xié)方差函數(shù)的性質(zhì)自協(xié)方差函數(shù)的性質(zhì)自相關(guān)系數(shù)白噪聲、白噪聲模擬Poisson白噪聲的60樣本的產(chǎn)生1. 隨機(jī)產(chǎn)生服從隨機(jī)產(chǎn)生服從(0,1)上均勻的上均勻的200個(gè)樣本個(gè)樣本:2. 給出服從參數(shù)為給出服從參數(shù)為1的指數(shù)分布的的指數(shù)分布的200個(gè)獨(dú)立個(gè)獨(dú)立樣本樣本;3. 給出參數(shù)為給出參數(shù)為1的的Poisson過(guò)程一條樣本軌道過(guò)程一條樣本軌道在在i=1,61上的取值；上的取值；參數(shù)為參數(shù)為1的的Poisson白噪聲的白噪聲的60個(gè)個(gè)樣本樣本I0

5、102030405060-1-0.500.511.522.533.54樣本樣本II0102030405060-1-0.500.511.52例：布朗運(yùn)動(dòng)標(biāo)準(zhǔn)正態(tài)白噪聲的標(biāo)準(zhǔn)正態(tài)白噪聲的60個(gè)樣本個(gè)樣本: A=randn(1,60)；plot(A)在經(jīng)濟(jì)工作中我們遇到的時(shí)間序列通常未必是平在經(jīng)濟(jì)工作中我們遇到的時(shí)間序列通常未必是平穩(wěn)的，那么如何檢驗(yàn)與鑒別呢？首先我們從圖穩(wěn)的，那么如何檢驗(yàn)與鑒別呢？首先我們從圖形上直觀地理解一下平穩(wěn)性與非平穩(wěn)性：形上直觀地理解一下平穩(wěn)性與非平穩(wěn)性：例例 某市某市1985-1994年各月的工業(yè)產(chǎn)值年各月的工業(yè)產(chǎn)值： a=10.93,9.34,11,10.98,11.2

6、9,11.84,10.62,10.9,12.77,12.15,12.24,12.3, 9.91,10.24,10.41,10.47,11.51,12.45,11.32,11.73,12.61,13.04,13.14,14.15, 10.85,10.3,12.74,12.73,13.08,14.27,13.18,13.75,14.42,13.95,14.53,14.91, 12.94,11.43,14.36,14.57,14.25,15.86,15.18,15.94,16.54,16.9,16.88,18.1, 13.7,10.88,15.79,16.36,17.22,17.75,16.62,1

7、6.96,17.69,16.4,17.51,19.73, 13.73,12.85,15.68,16.79,17.59,18.51,16.8,17.27,20.83,19.18,21.4,23.76, 15.73,13.14,17.24,17.93,18.82,19.12,17.7,19.87,21.17,21.44,22.14,22.45, 17.88,16,20.29,21.03,21.78,22.51,21.55,22.01,22.68,23.02,24.55,24.67, 19.61,17.15,22.46,23.19,23.40,26.26,22.91,24.03,23.94,24.1

8、2,25.87,28.25;作出其原始數(shù)據(jù)圖形與一階差分后的圖形作出其原始數(shù)據(jù)圖形與一階差分后的圖形解：解：c=a;b=c(:);plot(b)02040608010012051015202530原 始 數(shù) 據(jù)原始數(shù)據(jù)圖原始數(shù)據(jù)圖從圖從圖1.2可以發(fā)現(xiàn)該時(shí)間序列具有以下特征：可以發(fā)現(xiàn)該時(shí)間序列具有以下特征：隨著時(shí)間的推移，有逐漸增大的趨勢(shì)；隨著時(shí)間的推移，有逐漸增大的趨勢(shì)；具有一定的周期性；當(dāng)然還有隨機(jī)性。這是非具有一定的周期性；當(dāng)然還有隨機(jī)性。這是非平穩(wěn)時(shí)間序列平穩(wěn)時(shí)間序列. 在時(shí)間序列建模的方法中很多都是對(duì)平在時(shí)間序列建模的方法中很多都是對(duì)平穩(wěn)序列作出的，因此對(duì)于非平穩(wěn)時(shí)間序列穩(wěn)序列作出的

9、，因此對(duì)于非平穩(wěn)時(shí)間序列常常先將其化為平穩(wěn)序列，差分就是最常常常先將其化為平穩(wěn)序列，差分就是最常用的一種方法，對(duì)于圖用的一種方法，對(duì)于圖1.2中的數(shù)據(jù)我們利中的數(shù)據(jù)我們利用一次差分得到的數(shù)據(jù)圖形如圖用一次差分得到的數(shù)據(jù)圖形如圖1.3所示：所示：020406080100120-10-50510一 階 差 分 圖圖1-3 經(jīng)過(guò)一階差分處理后的圖經(jīng)過(guò)一階差分處理后的圖注意：注意： 時(shí)間趨勢(shì)可以通過(guò)差分解決，因?yàn)槎鄷r(shí)間趨勢(shì)可以通過(guò)差分解決，因?yàn)槎囗?xiàng)式差分一次，其階數(shù)就降低一次；項(xiàng)式差分一次，其階數(shù)就降低一次； 指數(shù)指數(shù)趨勢(shì)可以通過(guò)先取對(duì)數(shù)，然后再差分。趨勢(shì)可以通過(guò)先取對(duì)數(shù)，然后再差分。例如：例如：yt

10、=e2t，取對(duì)數(shù)變?yōu)?，取?duì)數(shù)變?yōu)閘nyt=2t，再做一次差分得，再做一次差分得lnyt+1-lnyt=2(t+1-t)=2此時(shí)為常數(shù)故平穩(wěn)此時(shí)為常數(shù)故平穩(wěn)三三 平穩(wěn)時(shí)間序列建模平穩(wěn)時(shí)間序列建模1.常見(jiàn)的三類(lèi)時(shí)間序列模型常見(jiàn)的三類(lèi)時(shí)間序列模型(1) P階自回歸模型：階自回歸模型：AR(p)t1t 12t 2ptptxxxx tt1t 12t 2qt qx 其中其中 t是白噪聲序列是白噪聲序列.其中其中 t是白噪聲序列，且與是白噪聲序列，且與xk不相關(guān)（不相關(guān)（kq時(shí)時(shí) q2kii 11 T 12/()是否檢驗(yàn)偏相關(guān)是否截尾，考察檢驗(yàn)偏相關(guān)是否截尾，考察kp時(shí)時(shí) kk2 1 T/ 是否其中其中T是

11、時(shí)間序列的長(zhǎng)度（建模數(shù)據(jù)的個(gè)數(shù)）是時(shí)間序列的長(zhǎng)度（建模數(shù)據(jù)的個(gè)數(shù)） ARMA(p,q)由于由于ARMA模型的自相關(guān)、偏相關(guān)系數(shù)均不截尾，模型的自相關(guān)、偏相關(guān)系數(shù)均不截尾，不能用上述方法，所以采用以下準(zhǔn)則：不能用上述方法，所以采用以下準(zhǔn)則：AIC準(zhǔn)則（赤池準(zhǔn)則），準(zhǔn)則（赤池準(zhǔn)則），SIC準(zhǔn)則準(zhǔn)則 AIC準(zhǔn)則（赤池準(zhǔn)則）準(zhǔn)則（赤池準(zhǔn)則）2AIC p qp q2(p+q) T(, )ln(, )/ T2t2t 1p qTpq (, ) 其中其中00p qAIC p qAIC p q,(,)min(, )在在MATLAB中有如下命令：中有如下命令：am = aic(Model1,Model2,.)SI

12、C準(zhǔn)則：準(zhǔn)則： T2tt 1SIC nTnTT()ln|/|ln/ 其中其中n是參數(shù)的個(gè)數(shù)，是參數(shù)的個(gè)數(shù)， 為殘差。為殘差。4.模型的檢驗(yàn)?zāi)Ｐ偷臋z驗(yàn)通常有通常有異方差檢驗(yàn)、殘差自相關(guān)檢驗(yàn)、正態(tài)性檢驗(yàn)異方差檢驗(yàn)、殘差自相關(guān)檢驗(yàn)、正態(tài)性檢驗(yàn)?zāi)Ｐ蜋z驗(yàn)的目的是檢驗(yàn)殘差序列是否為白噪聲模型檢驗(yàn)的目的是檢驗(yàn)殘差序列是否為白噪聲 四四 非平穩(wěn)序列處理方法非平穩(wěn)序列處理方法(2) 作對(duì)數(shù)差分（先取自然對(duì)數(shù)，再求差分）作對(duì)數(shù)差分（先取自然對(duì)數(shù)，再求差分）(3) 對(duì)于長(zhǎng)記憶時(shí)間序列做分?jǐn)?shù)差分對(duì)于長(zhǎng)記憶時(shí)間序列做分?jǐn)?shù)差分常見(jiàn)模型為：常見(jiàn)模型為：ARIMA(p,d,q), d為正整數(shù)為正整數(shù)ARFIMA(p,d,q)

13、,d為正分?jǐn)?shù)為正分?jǐn)?shù)(1) 作差分（或季節(jié)差分）作差分（或季節(jié)差分）例、化學(xué)溶液濃度變化數(shù)據(jù)例、化學(xué)溶液濃度變化數(shù)據(jù)0204060801001201401601802001616.51717.51818.5020406080100120140160180200-1-0.500.511.5例 、Canadian lynx data（猞猁）例五、滬深例五、滬深1209（股指期貨）（股指期貨）0204060801001201401602200230024002500260027002800050100150-0.04-0.03-0.02-0.0100.010.020.03)(log) 1(log)(

14、tXtXtY例、國(guó)際航空公司的月客數(shù)例、國(guó)際航空公司的月客數(shù) 050100150100200300400500600700y2=log(y1); plot(y2);0501001504.64.855.25.45.65.866.26.46.6y3=diff(y2); y=y3(13:143)-y3(1:131);020406080100120140-0.2-0.15-0.1-0.0500.050.10.15數(shù) 據(jù) 處 理 后 的 圖五五 時(shí)間序列分析的時(shí)間序列分析的MATLAB實(shí)現(xiàn)實(shí)現(xiàn)1. AR(p)模型參數(shù)估計(jì)的最小二乘法模型參數(shù)估計(jì)的最小二乘法在在Matlab中實(shí)現(xiàn)自回歸模型參數(shù)估計(jì)方法如下

15、：中實(shí)現(xiàn)自回歸模型參數(shù)估計(jì)方法如下：% 使用使用Burg方法計(jì)算方法計(jì)算AR模型的參數(shù)估計(jì)模型的參數(shù)估計(jì) a,e = arburg(x,p)其中：其中：x是時(shí)間序列，是時(shí)間序列，p是選定的階數(shù)；是選定的階數(shù)；a是輸出的回歸是輸出的回歸模型的系數(shù)，模型的系數(shù)，e是是AR模型白噪聲的方差估計(jì)。模型白噪聲的方差估計(jì)。n如果選用不同的如果選用不同的p值，可以利用循環(huán)語(yǔ)句實(shí)現(xiàn)值，可以利用循環(huán)語(yǔ)句實(shí)現(xiàn)nfor i=1:p,a,e(i)=arburg(x,i);endn% Compute an estimate of AR model parameters using the covariance meth

16、odn a,e = arcov(x,p)n% Compute an estimate of AR model parameters using the Yule-Walker methodn n a,e = aryule(x,p) ARMA模型參數(shù)估計(jì)在模型參數(shù)估計(jì)在Matlab中的實(shí)現(xiàn)中的實(shí)現(xiàn)m = armax(data,orders)其中其中data為時(shí)間序列數(shù)據(jù)，為時(shí)間序列數(shù)據(jù)，orders=p,q即即ARMA模型的階數(shù)模型的階數(shù)m ：回歸模型的參數(shù)：回歸模型的參數(shù)對(duì)于上述的一個(gè)例題，如果選用對(duì)于上述的一個(gè)例題，如果選用ARMA(2,1)模型模型，則有，則有 m = armax(b,2,

17、1)輸出結(jié)果如下：輸出結(jié)果如下：A(q)y(t) = C(q)e(t)A(q) = 1 - 0.4188 q-1 - 0.1523 q-2 C(q) = 1 - 0.02786 q-1 Estimated using ARMAX from data set b Loss function 0.797581 and FPE 0.85135 1+d/nFPE=V1-d/n2tt(1-0.4188L-0.1523L )y =(1 - 0.02786L) 注：注：其中其中V是是Loss function值，值，d是模型參數(shù)個(gè)數(shù)，是模型參數(shù)個(gè)數(shù)，n是是建模數(shù)據(jù)個(gè)數(shù)。建模數(shù)據(jù)個(gè)數(shù)。 最后模型為：最后模型為

18、：即即tt-1t-2tt-1y =0.4188y +0.1523y + - 0.02786n在在MATLAB中有如下命令：中有如下命令：nam = aic(Model1,Model2,.)MATLAB中計(jì)算中計(jì)算AIC的公式為：的公式為：AIC=log(V)+2d/N其中，其中，V是模型輸出的是模型輸出的loss functionam輸出各種模型的輸出各種模型的AIC統(tǒng)計(jì)量的數(shù)值，然后統(tǒng)計(jì)量的數(shù)值，然后從中選擇最小者從中選擇最小者例如：例如：m1=armax(b,2,1) % 建立建立ARMA(2,1)模型模型m2=armax(b,3,2) % 建立建立ARMA(3,2)模型模型m3=arma

19、x(b,4,3) % 建立建立ARMA(4,3)模型模型am = aic(m1,m2,m3)n4.模型的檢驗(yàn)?zāi)Ｐ偷臋z驗(yàn)n模型檢驗(yàn)的目的是檢驗(yàn)殘差序列是否為白噪模型檢驗(yàn)的目的是檢驗(yàn)殘差序列是否為白噪聲，在聲，在MATLAB中檢驗(yàn)的步驟如下：中檢驗(yàn)的步驟如下： 計(jì)算殘差向量在計(jì)算殘差向量在MATLAB中計(jì)算模型殘差的命中計(jì)算模型殘差的命令如下：令如下：e=resid(m,data)其中其中m為所使用的模型（如上面的為所使用的模型（如上面的m1,m2,m3），data為時(shí)間序列數(shù)據(jù)為時(shí)間序列數(shù)據(jù) 進(jìn)行正態(tài)性檢驗(yàn)進(jìn)行正態(tài)性檢驗(yàn) h,p = jbtest(e)若若 h=0,接受正態(tài)分布，接受正態(tài)分布，h=1,拒絕正態(tài)分布拒絕正態(tài)分布