1、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,
2、N,M,R,FAI)INTEGER: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,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 DO2FAI(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=0 A2=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)XV=0DO I=1,NXV=XV+X
5、(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)END3输出结果为：XV= 22.5718R( 1)= .8383 FAI( 1)= .6094R( 2)= .4648 FAI( 2)= -.1669R( 3)= -.0148 FAI( 3)= -.0701R( 4)= -.4776 FA
6、I( 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)= .8211 FAI(11)= .1246R(12)= .9508 FAI(12)= .1798从而得到自回归方程为：12t1t10t9t8t7t 654321t x798.26.x7.05.x.04. 04569x 注意
7、：以上是距平值，加上平均值即为实际值。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:TAO,Q !TAO:落后时间;Q:滑动平均的阶数REAL(8),DIMENSION(N):
8、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,Q4DO I=1,N-TAOS(TAO)=S(TAO)+X(I)*X(I+TAO)END DOS(TAO)=S(TAO)/(N-TAO)END
9、 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=S2NN=0DO NN=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 DO END DOEP1=ABS
10、(S2A-S2B)EP2=MAXVAL(ABS(THIT-THITA)IF(EP1D12.5)R5WRITE(12,(“S2A=“,E12.5)S2AWRITE(12,(“THITA=“,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.0E-5REAL(8):
11、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+01THITA= .77908D-01 -.65949D-01滑动平均模型为： 2t1ttt a06594.a07983.aX3 阶：S2= .11905D+01R= -.82189D-01 .65269D-01 .23275D-01S2A
12、= .11770E+01THITA= .79343D-01 -.67884D-01 -.23542D-01滑动平均模型为： 3t2t1ttt a054.a06784.a07934.aX6.3.7 附注64 自回归滑动平均模型（ARMA ）66.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输入参数，整型变量，
13、为自回归的阶数。Q输入参数，整型变量，为滑动平均的阶数。M输入参数，整型变量，M=P+Q。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):
14、FAI !自回归系数REAL(8),DIMENSION(P,P):A !工作数组REAL(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 !
15、EPS:迭代的精度S=0DO 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 DO7B(I)=S(Q+I)END DOCALL GASJDN(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-IA3=A3+FAI(J)*FAI(J+I)END DOA2=A2+A
16、3*(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 DOS2A=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 DOEP1=ABS(S2A-S2B)EP2=MAXVAL(ABS(THIT-THITA)WRITE(*,*)S2A,EP1,EP2IF(E
17、P1DMAX)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 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
18、=1,N9IF(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=30 INTEGER,PARAMETER:P=2,Q=1,M=P+Q !P
19、自回归阶数，Q 滑动平均阶数REAL(8),DIMENSION(N):XREAL(8),DIMENSION(0:P):FAIREAL(8),DIMENSION(Q):THITA !滑动系数REAL(8):XV !X 的平均值REAL(8):EPSEPS=1.E-4OPEN(10,FILE=BEIJING.DAT)READ(10,*)XCLOSE(10)XV=0DO I=1,NXV=XV+X(I)END DOXV=XV/NX=X-XVCALL ARMA(X,N,P,Q,M,FAI,THITA,EPS)OPEN(12,FILE=ARMA.DAT)WRITE(12,(2X,“XV=“,F8.4)XVDO I=1,PWRITE(12,(“ FAI(“,I2,“)=“,F8.4)I,FAI(I)END DODO I=1,QWRITE(12,(“ THITA(“,I2,“)=“,F8.4)I,THITA(I)10END DOCLOSE(12)END计算结果：XV= -4.5467FAI( 1)= .4894FAI( 2)= .1055THITA( 1)= .5697因此，自回归滑动平均方程为：1tt2t1tt a5697.0x05.489.0x 注意：上式得到的结果加上平均值 XV 即为实际预报值。
