時(shí)間序列分析

1、第六章時(shí)間序列分析6.2自回歸模型(AR)自回歸模型中 最簡(jiǎn)單的是一階自回歸模型和二階自回歸模型。為節(jié)省篇幅，這里直接給出 階自回歸模型。6.2.1 功能求出p階自回歸方程的系數(shù)，從而得到 p階自回歸方程。6.2.2 方法說明6.2.3 子程序語句SUBROUTINE ARP(X,N,M,R,FAI)6.2.4 啞元說明X 輸入?yún)?shù)，一維實(shí)型數(shù)組，大小為N,存放觀測(cè)序列值。N 輸入?yún)?shù)，整型變量，為觀測(cè)序列的長(zhǎng)度。M輸入?yún)?shù)，整型變量，為自回歸的階數(shù)。R輸出參數(shù)，一維實(shí)型數(shù)組，存放自相關(guān)系數(shù)。FAI 輸出參數(shù)，二維實(shí)型數(shù)組，存放自回歸系數(shù)。6.2.5 子程序!落后時(shí)間!協(xié)方差!S2:方差，A1

2、,A2:中間變量SUBROUTINE ARP(X,N,M,R,FAI) INTEGER:TAO REAL(4),DIMENSION(N):X REAL(4),DIMENSION(M,M):FAI REAL(4),DIMENSION(M):R REAL(4),DIMENSION(M):S REAL(4):S2,A1,A2S=0DO TAO=1,MDO I=1,N-TAOS(TAO)=S(TAO)+X(I)*X(I+TAO)END DOS(TAO)=S(TAO)/(N-TAO)END DOS2=0DO I=1,NS2=S2+X(I)*X(I)END DOS2=S2/NDO TAO=1,MR(TAO

3、)=0DO I=1,N-TAOR(TAO)=R(TAO)+X(I)*X(I+TAO)/S2END DOR(TAO)=R(TAO)/(N-TAO)END DOFAI(1,1)=R(1)FAI(2,2)=(R(2)-R(1)*R(1)/(1-R(1)*R(1)FAI(1,2)=FAI(1,1)-FAI(2,2)*FAI(1,1)DO J=3,MA1=0A2=0DO K=1,J-1A1=A1+FAI(K,J-1)*R(J-K)A2=A2+FAI(K,J-1)*R(K)END DOFAI(J,J)=(R(J)-A1)/(1-A2)DO K=1,J-1FAI(K,J)=FAI(K,J-1)-FAI(J,

4、J)*FAI(J-K,J-1)END DOEND DOEND6.2.6 例以某海區(qū)的 22 年的逐月氣溫為例，計(jì)算出自回歸系數(shù)，并給出自回歸方程。PROGRAM MAININTEGER,PARAMETER:N=264INTEGER,PARAMETER:M=12REAL(4),DIMENSION(N):XREAL(4),DIMENSION(M,M):FAIREAL(4),DIMENSION(M):RREAL(4):XV !X 的平均值OPEN(10,FILE='AA2.DAT')DO I=1,NREAD(10,'(F8.2)')X(I)END DOCLOSE(10

5、)XV=0DO I=1,NXV=XV+X(I)END DOXV=XV/NX=X-XVCALL ARP(X,N,M,R,FAI)OPEN(12,FILE='ARP.DAT')WRITE(12,'(2X,"XV=",F8.4)')XVDO I=1,MWRITE(12,'("R(",I2,")=",F8.4," FAI(",I2,")=",F8.4)')I,R(I),I,FAI(I,M)END DOCLOSE(12)END輸出結(jié)果為：XV= 22.571

6、8R( 1)=.8383FAI( 1)=.6094R( 2)=.4648FAI( 2)=-.1669R( 3)=-.0148FAI( 3)=-.0701R( 4)=-.4776FAI( 4)=-.0564R( 5)=-.8080FAI( 5)=-.1197R( 6)=-.9222FAI( 6)=.0477R( 7)=-.8019FAI( 7)=-.0471R( 8)=-.4747FAI( 8)=-.1702R( 9)=-.0108FAI( 9)=.0053R(10)=.4665FAI(10)=.0977R(11)=.8211FAI(11)=.1246R(12)=.9508FAI(12)=.17

7、98從而得到自回歸方程為：x t .6094x t 1 .1669x t 2 .0701x t 3 .0564x t 4 .1197x t 5 .0477x t 6.0471x t 7 .1702x t 8 .0053x t 9 .0977x t 10 .1246x t 11 .1798x t 12 注意：以上是距平值，加上平均值即為實(shí)際值。63 滑動(dòng)平均模型( MA)6.3.1 功能求出q階滑動(dòng)平均模型方程的系數(shù)，從而得到q階滑動(dòng)平均方程。6.3.2 方法說明6.3.3 子程序語句SUBROUTINE MAQ(X,N,Q,EPS)6.3.4 啞元說明X 輸入?yún)?shù)，實(shí)型一維數(shù)組，大小為N ，存

8、放觀測(cè)序列值。N 輸入?yún)?shù)，整型變量，數(shù)組的長(zhǎng)度。Q輸入?yún)?shù)，整型變量，滑動(dòng)平均的階數(shù)。EPS輸入?yún)?shù)，實(shí)型變量，存放迭代精度。6.3.5 子程序SUBROUTINE MAQ(X,N,Q,EPS) INTEGER:TAO,Q REAL(8),DIMENSION(N):X REAL(8),DIMENSION(Q):THITA REAL(8),DIMENSION(Q):THIT REAL(8),DIMENSION(Q):R REAL(8),DIMENSION(Q):S REAL(8):S2,A1REAL(8):S2AREAL(8):EPS,EP1,EP2S=0DO TAO=1,Q!TAO:落后時(shí)間

9、；Q:滑動(dòng)平均的階數(shù)! 滑動(dòng)系數(shù)!迭代中用的滑動(dòng)系數(shù) ,中間變量! 相關(guān)系數(shù)!S 協(xié)方差!S2:方差,A1 :中間變量!S2A:序列a(t)的方差!EPS:迭代的精度DO I=1,N-TAOS(TAO)=S(TAO)+X(I)*X(I+TAO)END DOS(TAO)=S(TAO)/(N-TAO)END DOS2=0DO I=1,NS2=S2+X(I)*X(I)END DOS2=S2/NDO TAO=1,QR(TAO)=0DO I=1,N-TAOR(TAO)=R(TAO)+X(I)*X(I+TAO)END DOR(TAO)=R(TAO)/S2/(N-TAO)END DOTHIT=0S2B=S2

10、NN=0DONN=NN+1A1=1DO I=1,QA1=A1+THIT(I)*THIT(I)END DOS2A=S2/A1THITA=-R*S2/S2ADO K=1,Q-1DO I=1,Q-KTHITA(K)=THITA(K)+THIT(I)*THIT(K+I)END DOEND DOEP1=ABS(S2A-S2B)EP2=MAXV AL(ABS(THIT-THITA)IF(EP1<EPS.AND.EP2<EPS)EXITTHIT=THITAS2B=S2APRINT*,'NN=',NNEND DOOPEN(12,FILE='MAQ.DAT')WRIT

11、E(12,*)WRITE(12,'("S2=",D12.5)')S2WRITE(12,'("R=",<Q>D12.5)')RWRITE(12,'("S2A=",E12.5)')S2AWRITE(12,'("THITA=",<Q>D12.5)')THITACLOSE(12)END6.3.6 例計(jì)算北京 1951 年 1980 年 1 月的平均氣溫 2 階、3 階滑動(dòng)平均模型的系數(shù)(同時(shí)也算出 了 12 月、2 月的結(jié)果)PROGR

12、AM MAININTEGER,PARAMETER:N=30INTEGER,PARAMETER:Q=2REAL(8),DIMENSION(N):XREAL(8),PARAMETER:EPS=1.0E-5REAL(8):XV !X 的平均值OPEN(10,FILE='BEIJING .DAT')READ(10,*)XCLOSE(10)XV=0DO I=1,NXV=XV+X(I)END DOXV=XV/NX=X-XVCALL MAQ(X,N,Q,EPS)END 計(jì)算結(jié)果為： 2 階：S2= .11905D+01R= -.82189D-01 .65269D-01S2A= .11782E

13、+01THITA= .77908D-01 -.65949D-01 滑動(dòng)平均模型為：X t at 0.0779083at 1 0.065949at 23 階：S2= .11905D+01R= -.82189D-01 .65269D-01 .23275D-01S2A= .11770E+01THITA= .79343D-01 -.67884D-01 -.23542D-01 滑動(dòng)平均模型為：X t at 0.079343at 1 0.067884at 2 0.023542at 36.3.7 附注64自回歸滑動(dòng)平均模型( ARMA )6.4.1 功能求出(p, q)階自回歸一滑動(dòng)平均方程的系數(shù)，從而得到

14、(p, q)階自回歸一滑動(dòng)平均方程。6.4.2 方法說明6.4.3 子程序語句SUBROUTINE ARMA(X,N,P,Q,M,R,FAI,THITA,EPS)6.4.4 啞元說明X 輸入?yún)?shù)，一維實(shí)型數(shù)組，大小為N,存放觀測(cè)序列值。N 輸入?yún)?shù)，整型變量，為觀測(cè)序列的長(zhǎng)度。P輸入?yún)?shù)，整型變量，為自回歸的階數(shù)。Q輸入?yún)?shù)，整型變量，為滑動(dòng)平均的階數(shù)。M 輸入?yún)?shù)，整型變量，M=P+Q。R輸出參數(shù)，一維實(shí)型數(shù)組，存放自相關(guān)系數(shù)。FAI 輸出參數(shù)，一維實(shí)型數(shù)組，存放自回歸系數(shù)。THITA 輸出參數(shù)，一維實(shí)型數(shù)組，存放滑動(dòng)平均系數(shù)。EPS實(shí)型常量，存放迭代時(shí)要求的精度。6.4.5 子程序SUBR

15、OUTINE ARMA(X,N,P,Q,M,FAI,THITA,EPS)INTEGER:TAO!落后時(shí)間INTEGER:P!自回歸階數(shù)INTEGER:Q!滑動(dòng)平均階數(shù)INTEGER:M!M=P+QREAL(8),DIMENSION(N):XREAL(8),DIMENSION(0:P):FAIREAL(8),DIMENSION(P,P):AREAL(8),DIMENSION(P):BREAL(8),DIMENSION(0:M):SREAL(8),DIMENSION(0:Q):SCREAL(8),DIMENSION(Q):THITAREAL(8),DIMENSION(Q):THITREAL(8):

16、A1,A2,A3!A1,A2,A3:REAL(8):S2A!S2A:REAL(8):EPS,EP1,EP2 !EPS:S=0輸入序列自回歸系數(shù)工作數(shù)組工作數(shù)組協(xié)方差， S(0) 即為方差自回歸后的協(xié)方差 滑動(dòng)平均系數(shù) 迭代中用的滑動(dòng)系數(shù) , 中間變量中間變量 自回歸后的序列 a(t) 的方差 迭代的精度DO TAO=0,MDO I=1,N-TAOS(TAO)=S(TAO)+X(I)*X(I+TAO)END DOS(TAO)=S(TAO)/(N-TAO)END DODO I=1,PDO J=1,PA(I,J)=S(ABS(Q+I-J)END DOB(I)=S(Q+I)END DOCALL GAS

17、JDN(A,B,P) FAI(1:P)=B(1:P) FAI(0)=-1A1=0DO I=0,PA1=A1+FAI(I)*FAI(I)END DODO K=0,QA2=0DO I=1,PA3=0DO J=0,P-I A3=A3+FAI(J)*FAI(J+I)END DOA2=A2+A3*(S(K+I)+S(ABS(K-I)END DOSC(K)=A1*S(K)+A2END DOS2B=0THIT=0NN=0DO ! 迭代NN=NN+1WRITE(*,'(" NN=",I3)')NNA1=1DO I=1,QA1=A1+THIT(I)*THIT(I)END DO

18、S2A=SC(0)/A1DO K=1,QTHITA(K)=-SC(K)/S2ADO I=1,Q-KTHITA(K)=THITA(K)+THIT(I)*THIT(K+I) END DO END DO EP1=ABS(S2A-S2B) EP2=MAXVAL(ABS(THIT-THITA) WRITE(*,*)S2A,EP1,EP2 IF(EP1<EPS.AND.EP2<EPS)EXIT THIT=THITA S2B=S2AEND DOEND全選主元高斯約當(dāng)法( Gauss-Jordan )求解 n 階線性代數(shù)方程組 SUBROUTINE GASJDN(A,B,N)REAL(8),DIM

19、ENSION(N,N):AREAL(8),DIMENSION(N):BREAL(8),DIMENSION(N):JAREAL(8):DMAX,DDLL=1DO K=1,NDMAX=0DO I=K,NDO J=K,NIF(ABS(A(I,J)>DMAX)THENDMAX=ABS(A(I,J)JA(K)=JIA=IEND IFEND DOEND DOIF(DMAX+1=1)THENWRITE(*,'(" 主元為 0 ，求解失敗 ")')LL=0RETURNEND IFDO J=K,NDD=A(K,J)A(K,J)=A(IA,J)A(IA,J)=DDEND

20、DODD=B(K)B(K)=B(IA)B(IA)=DDDO I=1,NDD=A(I,K)A(I,K)=A(I,JA(K)A(I,JA(K)=DDEND DODO J=K+1,NA(K,J)=A(K,J)/A(K,K)END DOB(K)=B(K)/A(K,K)DO J=K+1,NDO I=1,NIF(I/=K)A(I,J)=A(I,J)-A(I,K)*A(K,J)END DOEND DODO I=1,NIF(I/=K)THENB(I)=B(I)-A(I,K)*B(K)ENDIFEND DOEND DODO K=N,1,-1DD=B(K)B(K)=B(JA(K) B(JA(K)=DDEND DOEND6.4.6 例以 7.3.6 中資料為例，計(jì)算北京 1951年 1980 年 1月的平均氣溫 2階字回歸和 1 階滑動(dòng) 平均模型的系數(shù)。PROGRAM ARMAPQINTEGER,PARAMETER:N=30INTEGER,PARAMETER:P=2,