时间序列分析

1、第六章时间序列分析6.2自回归模型(AR)自回归模型中 最简单的是一阶自回归模型和二阶自回归模型。为节省篇幅，这里直接给出 阶自回归模型。6.2.1 功能求出p阶自回归方程的系数，从而得到 p阶自回归方程。6.2.2 方法说明6.2.3 子程序语句SUBROUTINE ARP(X,N,M,R,FAI)6.2.4 哑元说明X 输入参数，一维实型数组，大小为N,存放观测序列值。N 输入参数，整型变量，为观测序列的长度。M输入参数，整型变量，为自回归的阶数。R输出参数，一维实型数组，存放自相关系数。FAI 输出参数，二维实型数组，存放自回归系数。6.2.5 子程序!落后时间!协方差!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 例以某海区的 22 年的逐月气温为例，计算出自回归系数，并给出自回归方程。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输出结果为：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 注意：以上是距平值，加上平均值即为实际值。63 滑动平均模型( MA)6.3.1 功能求出q阶滑动平均模型方程的系数，从而得到q阶滑动平均方程。6.3.2 方法说明6.3.3 子程序语句SUBROUTINE MAQ(X,N,Q,EPS)6.3.4 哑元说明X 输入参数，实型一维数组，大小为N ，存

8、放观测序列值。N 输入参数，整型变量，数组的长度。Q输入参数，整型变量，滑动平均的阶数。EPS输入参数，实型变量，存放迭代精度。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:落后时间

9、；Q:滑动平均的阶数! 滑动系数!迭代中用的滑动系数 ,中间变量! 相关系数!S 协方差!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 例计算北京 1951 年 1980 年 1 月的平均气温 2 阶、3 阶滑动平均模型的系数(同时也算出 了 12 月、2 月的结果)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 计算结果为： 2 阶：S2= .11905D+01R= -.82189D-01 .65269D-01S2A= .11782E

13、+01THITA= .77908D-01 -.65949D-01 滑动平均模型为：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 滑动平均模型为：X t at 0.079343at 1 0.067884at 2 0.023542at 36.3.7 附注64自回归滑动平均模型( ARMA )6.4.1 功能求出(p, q)阶自回归一滑动平均方程的系数，从而得到

14、(p, q)阶自回归一滑动平均方程。6.4.2 方法说明6.4.3 子程序语句SUBROUTINE ARMA(X,N,P,Q,M,R,FAI,THITA,EPS)6.4.4 哑元说明X 输入参数，一维实型数组，大小为N,存放观测序列值。N 输入参数，整型变量，为观测序列的长度。P输入参数，整型变量，为自回归的阶数。Q输入参数，整型变量，为滑动平均的阶数。M 输入参数，整型变量，M=P+Q。R输出参数，一维实型数组，存放自相关系数。FAI 输出参数，一维实型数组，存放自回归系数。THITA 输出参数，一维实型数组，存放滑动平均系数。EPS实型常量，存放迭代时要求的精度。6.4.5 子程序SUBR

15、OUTINE ARMA(X,N,P,Q,M,FAI,THITA,EPS)INTEGER:TAO!落后时间INTEGER:P!自回归阶数INTEGER:Q!滑动平均阶数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输入序列自回归系数工作数组工作数组协方差， S(0) 即为方差自回归后的协方差 滑动平均系数 迭代中用的滑动系数 , 中间变量中间变量 自回归后的序列 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全选主元高斯约当法( Gauss-Jordan )求解 n 阶线性代数方程组 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 中资料为例，计算北京 1951年 1980 年 1月的平均气温 2阶字回归和 1 阶滑动 平均模型的系数。PROGRAM ARMAPQINTEGER,PARAMETER:N=30INTEGER,PARAMETER:P=2,