时间序列分析

1、第六章时间序列分析6自回归模型（AR）自回归模型中最简单的是一阶自回归模型和二阶自回归模型。为节省篇幅，这里直接给出p阶自回归模型。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子程序SUBROUTINE ARP(X,N,M,R,FAI)IN

2、TEGER:TAO !落后时间REAL(4),DIMENSION(N):XREAL(4),DIMENSION(M,M):FAIREAL(4),DIMENSION(M):RREAL(4),DIMENSION(M):S!协方差REAL(4):S2,A1,A2 !S2:方差, A1,A2:中间变量S=0DO TAO=1,M DO I=1,N-TAO S(TAO)=S(TAO)+X(I)*X(I+TAO) END DO S(TAO)=S(TAO)/(N-TAO)END DOS2=0DO I=1,N S2=S2+X(I)*X(I)END DOS2=S2/NDO TAO=1,M R(TAO)=0 DO I

3、=1,N-TAO R(TAO)=R(TAO)+X(I)*X(I+TAO)/S2 END DO R(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,M A1=0 A2=0 DO K=1,J-1 A1=A1+FAI(K,J-1)*R(J-K) A2=A2+FAI(K,J-1)*R(K) END DO FAI(J,J)=(R(J)-A1)/(1-A2) DO K=1,J-1 FAI(K,J)=FAI(K,J-1)-F

4、AI(J,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

5、(10)XV=0DO I=1,N XV=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,M WRITE(12,'("R(",I2,")=",F8.4," FAI(",I2,")=",F8.4)')I,R(I),I,FAI(I,M)END DOCLOSE(12)END输出结果为：XV= 2

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

7、1 FAI(11)= .1246R(12)= .9508 FAI(12)= .1798从而得到自回归方程为： 注意：以上是距平值，加上平均值即为实际值。63滑动平均模型（MA）6.3.1功能求出q阶滑动平均模型方程的系数，从而得到q阶滑动平均方程。6.3.2方法说明6.3.3子程序语句SUBROUTINE MAQ(X,N,Q,EPS)6.3.4哑元说明X输入参数，实型一维数组，大小为N，存放观测序列值。N输入参数，整型变量，数组的长度。Q输入参数，整型变量，滑动平均的阶数。EPS输入参数，实型变量，存放迭代精度。6.3.5子程序SUBROUTINE MAQ(X,N,Q,EPS)INTEGER:

8、TAO,Q !TAO:落后时间;Q:滑动平均的阶数REAL(8),DIMENSION(N):XREAL(8),DIMENSION(Q):THITA !滑动系数REAL(8),DIMENSION(Q):THIT !迭代中用的滑动系数,中间变量REAL(8),DIMENSION(Q):R !相关系数REAL(8),DIMENSION(Q):S !S协方差REAL(8):S2,A1 !S2:方差, A1:中间变量REAL(8):S2A !S2A:序列a(t)的方差REAL(8):EPS,EP1,EP2 !EPS:迭代的精度S=0DO TAO=1,Q DO I=1,N-TAO S(TAO)=S(TAO

9、)+X(I)*X(I+TAO) END DO S(TAO)=S(TAO)/(N-TAO)END DOS2=0DO I=1,N S2=S2+X(I)*X(I)END DOS2=S2/NDO TAO=1,Q R(TAO)=0 DO I=1,N-TAO R(TAO)=R(TAO)+X(I)*X(I+TAO) END DO R(TAO)=R(TAO)/S2/(N-TAO)END DOTHIT=0S2B=S2NN=0DO NN=NN+1 A1=1 DO I=1,Q A1=A1+THIT(I)*THIT(I) END DO S2A=S2/A1 THITA=-R*S2/S2A DO K=1,Q-1 DO I

10、=1,Q-K THITA(K)=THITA(K)+THIT(I)*THIT(K+I) END DO END DO EP1=ABS(S2A-S2B) EP2=MAXVAL(ABS(THIT-THITA) IF(EP1<EPS.AND.EP2<EPS)EXIT THIT=THITA S2B=S2A PRINT*,'NN=',NNEND DOOPEN(12,FILE='MAQ.DAT')WRITE(12,*)WRITE(12,'("S2=",D12.5)')S2WRITE(12,'("R=",

11、<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月的结果）PROGRAM MAININTEGER,PARAMETER:N=30INTEGER,PARAMETER:Q=2REAL(8),DIMENSION(N):XREAL(8),PARAMETER:EPS=1.

12、0E-5REAL(8):XV !X的平均值OPEN(10,FILE='BEIJING.DAT')READ(10,*)XCLOSE(10)XV=0DO I=1,N XV=XV+X(I)END DOXV=XV/NX=X-XVCALL MAQ(X,N,Q,EPS)END计算结果为：2阶：S2= .11905D+01R= -.82189D-01 .65269D-01S2A= .11782E+01THITA= .77908D-01 -.65949D-01滑动平均模型为：3阶：S2= .11905D+01R= -.82189D-01 .65269D-01 .23275D-01S2A= .1

13、1770E+01THITA= .79343D-01 -.67884D-01 -.23542D-01滑动平均模型为：6.3.7附注64自回归滑动平均模型（ARMA）6.4.1功能求出（p，q）阶自回归滑动平均方程的系数，从而得到（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

14、。 R输出参数，一维实型数组，存放自相关系数。FAI输出参数，一维实型数组，存放自回归系数。THITA输出参数，一维实型数组，存放滑动平均系数。EPS实型常量，存放迭代时要求的精度。6.4.5子程序SUBROUTINE ARMA(X,N,P,Q,M,FAI,THITA,EPS)INTEGER:TAO !落后时间INTEGER:P !自回归阶数INTEGER:Q !滑动平均阶数INTEGER:M !M=P+QREAL(8),DIMENSION(N):X !输入序列REAL(8),DIMENSION(0:P):FAI !自回归系数REAL(8),DIMENSION(P,P):A !工作数组REAL

15、(8),DIMENSION(P):B !工作数组REAL(8),DIMENSION(0:M):S !协方差，S(0)即为方差REAL(8),DIMENSION(0:Q):SC !自回归后的协方差REAL(8),DIMENSION(Q):THITA !滑动平均系数REAL(8),DIMENSION(Q):THIT !迭代中用的滑动系数,中间变量REAL(8):A1,A2,A3 !A1,A2,A3:中间变量REAL(8):S2A !S2A:自回归后的序列a(t)的方差REAL(8):EPS,EP1,EP2 !EPS:迭代的精度S=0DO TAO=0,M DO I=1,N-TAO S(TAO)=S(

16、TAO)+X(I)*X(I+TAO) END DO S(TAO)=S(TAO)/(N-TAO)END DODO I=1,P DO J=1,P A(I,J)=S(ABS(Q+I-J) END DO B(I)=S(Q+I)END DOCALL GASJDN(A,B,P)FAI(1:P)=B(1:P)FAI(0)=-1A1=0DO I=0,P A1=A1+FAI(I)*FAI(I)END DODO K=0,Q A2=0 DO I=1,P A3=0 DO J=0,P-I A3=A3+FAI(J)*FAI(J+I) END DO A2=A2+A3*(S(K+I)+S(ABS(K-I) END DO SC

17、(K)=A1*S(K)+A2END DOS2B=0THIT=0NN=0DO !迭代 NN=NN+1 WRITE(*,'(" NN=",I3)')NN A1=1 DO I=1,Q A1=A1+THIT(I)*THIT(I) END DO S2A=SC(0)/A1 DO K=1,Q THITA(K)=-SC(K)/S2A DO I=1,Q-K THITA(K)=THITA(K)+THIT(I)*THIT(K+I) END DO END DO EP1=ABS(S2A-S2B) EP2=MAXVAL(ABS(THIT-THITA) WRITE(*,*)S2A,EP1

18、,EP2 IF(EP1<EPS.AND.EP2<EPS)EXIT THIT=THITA S2B=S2AEND DOEND!全选主元高斯约当法（Gauss-Jordan）求解n阶线性代数方程组SUBROUTINE GASJDN(A,B,N)REAL(8),DIMENSION(N,N):AREAL(8),DIMENSION(N):BREAL(8),DIMENSION(N):JAREAL(8):DMAX,DDLL=1DO K=1,N DMAX=0 DO I=K,N DO J=K,N IF(ABS(A(I,J)>DMAX)THEN DMAX=ABS(A(I,J)JA(K)=JIA=I

19、 END IFEND DO END DO IF(DMAX+1=1)THEN WRITE(*,'(" 主元为0，求解失败 ")')LL=0RETURN END IF DO J=K,N DD=A(K,J) A(K,J)=A(IA,J) A(IA,J)=DD END DO DD=B(K) B(K)=B(IA) B(IA)=DD DO I=1,N DD=A(I,K) A(I,K)=A(I,JA(K) A(I,JA(K)=DD END DO DO J=K+1,N A(K,J)=A(K,J)/A(K,K) END DO B(K)=B(K)/A(K,K) DO J=K+1,N DO I=1,N IF(I/=K)A(I,J)=A(I,J)-A(I,K)*A(K,J) END DO END DO DO I=1,N IF(I/=K)THEN B(I)=B(I)-A(I,K)*B(K)ENDIF END DOEND DODO K=N,1,-1 DD=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=30 I