1、r語言：時間序列ARMA基礎學習26 二月, 2013,oldlee11,R語言與數據挖掘,時間序列分析,?010203040506070809101112131415161718192021222324252627282930313233343536# 術語 #白噪音：其均值=0，并且獨立分布的（同時間無關）#穩定時間序列：任意j-i時間段的序列：其均值相等#自相關系數acf圖：研究yt同yt-l序列之間的相關性# 在純的ma(q)序列下，acf圖形表現為q+1以后的自相關系數約為0（虛線內）#偏相關系數pacf圖：在yt同yt-l之間的序列固定的情況下，研究研究yt同yt-l序列之間的相關

2、性# 在純的ar(p)序列下，pacf圖形表現為p+1以后的偏相關系數約為0（虛線內）#擴展相關系數圖eacf：如果yt不是純的ar或ma，而是arma（混合體），無法通過acf確定q，也不能通過pacf確認p，需要通過eacf確認p和q# 模擬產生ma ar arma 序列 # MA 時間序列的模擬試驗：產生一個ma時間序列y.ma-function(a1,a2,a3=0,a4=0,num=200,pic=TRUE)#MA滑動平均時間序列的模擬（也可以使用filter函數）e-rnorm(num,0,1)#模擬白噪聲,均值=0result-0result1-e1result2-e2-a1*e

3、1result3-e3-a1*e2-a2*e1result4-e4-a1*e3-a2*e2-a3*e1for(t in 5:num) resultt-et-a1*et-1-a2*et-2-a3*et-3-a4*et-4 #構造一個ma型時間序列if(pic=TRUE)#畫圖形dev.new()ts.plot(result,main=paste(y.mat=et-,a1,*et-1-,a2,*et-2-,a3,*et-3-,a4,*et-4的時間序列散點圖)dev.new()lag.plot(result, 9, do.lines=FALSE)dev.new()par(mfrow=c(2,1)a

4、cf(result, 30,main=paste(y.ma自相關圖,y.mat=et-,a1,*et-1-,a2,*et-2-,a3,*et-3-,a4,*et-4)pacf(result, 30,main=paste(y.ma偏自相關圖,y.mat=et-,a1,*et-1-,a2,*et-2-,a3,*et-3-,a4,*et-4)resulty.ma-y.ma(0.92,0.65)結果一：繪制散點圖結果二：繪制出yt同yt-1(延遲1)、yt-2(延遲2).yt-9(延遲9)的2維散點圖，用以觀察yt同yt-i的相關性（i=19）結果三：繪制自相關和偏自相關圖（在虛線外的表示有相關性）可

5、以看到：1）在純的ma(q)序列下，acf圖形表現為q+1以后的自相關系數約為0（虛線內）2）在純的ma(q)序列下，pacf則不規則的會大于1/-1可以通過acf圖確定q的數值?01020304050607080910111213141516171819202122# AR 時間序列的模擬試驗：產生一個ar時間序列y.ar-function(b1,b2,b3=0,b4=0,num=200,pic=TRUE)#AR自回歸型時間序列的模擬（也可以使用filter函數）e-rnorm(num,0,1)#模擬白噪聲,均值=0result-0result1-e1result2-e2+b1*result

6、1result3-e3+b1*result2+b2*result1result4-e4+b1*result3+b2*result2+b3*result1for(t in 5:num) resultt-et+b1*resultt-1+b2*resultt-2+b3*resultt-3+b4*resultt-4 #構造一個ar型時間序列if(pic=TRUE)#畫圖形dev.new()ts.plot(result,main=paste(y.art=et+,b1,*y.art-1+,b2,*y.art-2+,b3,*y.art-3+,b4,*y.art-4的時間序列散點圖)dev.new()lag.

7、plot(result, 9, do.lines=FALSE)dev.new()par(mfrow=c(2,1)acf(result, 30,main=paste(y.ar自相關圖,y.art=et+,b1,*y.art-1+,b2,*y.art-2+,b3,*y.art-3+,b4,*y.art-4)pacf(result, 30,main=paste(y.ar偏自相關圖,y.art=et+,b1,*y.art-1+,b2,*y.art-2+,b3,*y.art-3+,b4,*y.art-4)resulty.ar-y.ar(0.92,-0.65)結果一：繪制散點圖結果二：繪制出yt同yt-1

8、(延遲1)、yt-2(延遲2).yt-9(延遲9)的2維散點圖，用以觀察yt同yt-i的相關性（i=19）結果三：繪制自相關和偏自相關圖（在虛線外的表示有相關性）1）在純的ar(q)序列下，pacf圖形表現為q+1以后的自相關系數約為0（虛線內）2）在純的ar(q)序列下，acf則不規則的會大于1/-1可以通過pacf圖確定p的數值?0102030405060708091011121314151617181920# ARMA 時間序列的模擬試驗：產生一個arma時間序列library(TSA)y.arma-function(a1,a2,a3=0,a4=0,b1,b2,b3=0,b4=0,num

9、=200,pic=TRUE)result-y.ma(a1=a1,a2=a2,a3=a3,a4=a4,pic=F,num=num)+y.ar(b1=b1,b2=b2,b3=b3,b4=b4,pic=F,num=num)#產生自回歸滑動平均時間序列ARIMAexp.str=q+1的自相關系數全部約小于0 | 逐步將接近0# AR(p) 逐步將接近0 | lag=p+1的自相關系數全部約小于0 # ARMA(p,q) 逐步將接近0 | 逐步將接近0# step2 序列的類型確定后：求p q# 如果是MA(q) 類型的，看acf圖，看前幾個不為0# 如果是AR(p) 類型的，看pacf圖，看前幾個不為

10、0# 如果是ARMA(p,q) 類型的，看eacf圖，看x三角形的頂端#代碼：data是向量數據library(TSA)plot.cf-function(data)dev.new()ts.plot(data,main=時間序列散點圖)dev.new()lag.plot(data, 9, do.lines=FALSE)dev.new()par(mfrow=c(2,1)acf(data, 30,main=自相關圖)pacf(data, 30,main=偏自相關圖)eacf(data)data-y.arma(a1=0.92,a2=0.65,b1=0.92,b2=-0.65,num=100,pic=F

11、)plot.cf(data)# 產生ariam模型 # p AR的參數# q MA的參數# d 移動取差的階數# arima(data.ts,order=c(p,d,q)#產生模型。arima只處理時間序列，不處理向量等p=2;d=0;q=2data-y.arma(a1=0.92,a2=0.65,b1=0.92,b2=-0.65,num=100,pic=F)data.ts-ts(data,freq=1,star=1)#把向量化為時間序列sol-arima(data.ts,order=c(p,d,q)#產生模型。arima只處理時間序列，不處理向量等dev.new()plot(sol,n.ahe

12、ad=30)#預測30個點# 檢驗預測模型的性能 # 序列的正態分布驗證函數norm.test-function(input.data,alpha=0.05,pic=TRUE)if(pic=TRUE)#畫圖形dev.new()par(mfrow=c(2,1)qqnorm(input.data,main=qq圖)qqline(input.data)hist(input.data,freq=F,main=直方圖和密度估計曲線)lines(density(input.data),col=blue)#密度估計曲線x-c(round(min(input.data):round(max(input.dat

13、a)lines(x,dnorm(x,mean(input.data),sd(input.data),col=red)#正態分布曲線solalpha)print(paste(success:服從正態分布，p.value=,sol$p.value,alpha)elseprint(paste(error:不服從正態分布，p.value=,sol$p.value,=,alpha) sol#殘差分析#方案一（自己編寫的，不一定好）# ariam.sol由ariam()產生的模型aya.res-function(ariam.sol)#step1 ：檢驗正太性，如果符合正太，表示均值相同：穩定的dev.ne

14、w()norm.test(ariam.sol$residuals)#畫出QQ圖,并給出檢驗正態分布的測試結果par(mfrow=c(3,1)#step2 ：查看時序，如果有趨勢變化，則mean不為相同，即不是穩定的plot(ariam.sol$residuals,main=殘差序列圖,type=o)#畫出時序圖abline(h=0)#step3 ：查看acf圖，如果全部在虛線內，則表示殘差和時間無關系acf(ariam.sol$residuals, 30,main=殘差的自相關圖:如果所有lag的系數均處于虛線內，則原始模型較好)#step4 ：使用Ljung來進行白噪聲驗證：表示均值全部為0

15、，即殘差很小test.p-0for(i in 1:20)test.pi-Box.test(sol$residuals,i,fitdf=2,type=Ljung)$p.value#選取殘差中lag=1-i的自相關系數來建立Ljung量。print(test.p)plot(c(NA,NA,test.p-c(1:2),main=殘差的Ljung檢驗p.value:如果lag=3-以后的所有點均大于0.05，則是白噪聲,ylim=c(0,1)abline(h=0.05,col=red) aya.res(sol)#方案二(好，書上的)#step1 ：檢驗正太性，如果符合正太，表示均值相同：穩定的dev.

16、new()norm.test(sol$residuals)#畫出QQ圖,并給出檢驗正態分布的測試結果dev.new()#使用對arima產生的模型，使用tsdiag()函數，直接畫出：殘差的時序圖，殘差的acf圖，Ljung產生的p.value（和0.05比較）tsdiag(sol,gof=15,omit.initial=F)# sol是由ariam()產生的模型# 使用ariam模型進行預測 # arima.pred函數：# ariam.sol由ariam()產生的模型# nahead預測的點數# alpha計算預測區間時的置信度：1-alpha# 輸出的數據為數據框：$pred是預測的估值# $max是預測的估值的最大值（1-alpha置信度，雙側估計）# $min是預測的估值的最小值（1-alpha置信度，雙側估計）arima.pred-function(ariam.sol,nahead=5,alpha=0.05)pred-predict(ariam.sol,n.ahead=nahead)#產生nahead個預測數據，包括預測的平均值$pred.和預測標準誤$sek-qnorm(1-alpha/2),0,1)pred.data-data.frame(pred=pred$pred,max