时间序列分析.docx
《时间序列分析.docx》由会员分享,可在线阅读,更多相关《时间序列分析.docx(14页珍藏版)》请在冰豆网上搜索。
时间序列分析
第六章时间序列分析
6.2自回归模型(AR)
自回归模型中最简单的是一阶自回归模型和二阶自回归模型。
为节省篇幅,这里直接给出p阶自回归模型。
6.2.1功能
求出p阶自回归方程的系数,从而得到p阶自回归方程。
6.2.2方法说明
6.2.3子程序语句
SUBROUTINEARP(X,N,M,R,FAI)
6.2.4哑元说明
X——输入参数,一维实型数组,大小为N,存放观测序列值。
N——输入参数,整型变量,为观测序列的长度。
M——输入参数,整型变量,为自回归的阶数。
R——输出参数,一维实型数组,存放自相关系数。
FAI——输出参数,二维实型数组,存放自回归系数。
6.2.5子程序
SUBROUTINEARP(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,A2!
S2:
方差,A1,A2:
中间变量
S=0
DOTAO=1,M
DOI=1,N-TAO
S(TAO)=S(TAO)+X(I)*X(I+TAO)
ENDDO
S(TAO)=S(TAO)/(N-TAO)
ENDDO
S2=0
DOI=1,N
S2=S2+X(I)*X(I)
ENDDO
S2=S2/N
DOTAO=1,M
R(TAO)=0
DOI=1,N-TAO
R(TAO)=R(TAO)+X(I)*X(I+TAO)/S2
ENDDO
R(TAO)=R(TAO)/(N-TAO)
ENDDO
FAI(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)
DOJ=3,M
A1=0
A2=0
DOK=1,J-1
A1=A1+FAI(K,J-1)*R(J-K)
A2=A2+FAI(K,J-1)*R(K)
ENDDO
FAI(J,J)=(R(J)-A1)/(1-A2)
DOK=1,J-1
FAI(K,J)=FAI(K,J-1)-FAI(J,J)*FAI(J-K,J-1)
ENDDO
ENDDO
END
6.2.6例
以某海区的22年的逐月气温为例,计算出自回归系数,并给出自回归方程。
PROGRAMMAIN
INTEGER,PARAMETER:
:
N=264
INTEGER,PARAMETER:
:
M=12
REAL(4),DIMENSION(N):
:
X
REAL(4),DIMENSION(M,M):
:
FAI
REAL(4),DIMENSION(M):
:
R
REAL(4):
:
XV!
X的平均值
OPEN(10,FILE='AA2.DAT')
DOI=1,N
READ(10,'(F8.2)')X(I)
ENDDO
CLOSE(10)
XV=0
DOI=1,N
XV=XV+X(I)
ENDDO
XV=XV/N
X=X-XV
CALLARP(X,N,M,R,FAI)
OPEN(12,FILE='ARP.DAT')
WRITE(12,'(2X,"XV=",F8.4)')XV
DOI=1,M
WRITE(12,'("R(",I2,")=",F8.4,"FAI(",I2,")=",F8.4)')I,R(I),I,FAI(I,M)
ENDDO
CLOSE(12)
END
输出结果为:
XV=22.5718
R
(1)=.8383FAI
(1)=.6094
R
(2)=.4648FAI
(2)=-.1669
R(3)=-.0148FAI(3)=-.0701
R(4)=-.4776FAI(4)=-.0564
R(5)=-.8080FAI(5)=-.1197
R(6)=-.9222FAI(6)=.0477
R(7)=-.8019FAI(7)=-.0471
R(8)=-.4747FAI(8)=-.1702
R(9)=-.0108FAI(9)=.0053
R(10)=.4665FAI(10)=.0977
R(11)=.8211FAI(11)=.1246
R(12)=.9508FAI(12)=.1798
从而得到自回归方程为:
注意:
以上是距平值,加上平均值即为实际值。
6.3滑动平均模型(MA)
6.3.1功能
求出q阶滑动平均模型方程的系数,从而得到q阶滑动平均方程。
6.3.2方法说明
6.3.3子程序语句
SUBROUTINEMAQ(X,N,Q,EPS)
6.3.4哑元说明
X——输入参数,实型一维数组,大小为N,存放观测序列值。
N——输入参数,整型变量,数组的长度。
Q——输入参数,整型变量,滑动平均的阶数。
EPS——输入参数,实型变量,存放迭代精度。
6.3.5子程序
SUBROUTINEMAQ(X,N,Q,EPS)
INTEGER:
:
TAO,Q!
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!
S协方差
REAL(8):
:
S2,A1!
S2:
方差,A1:
中间变量
REAL(8):
:
S2A!
S2A:
序列a(t)的方差
REAL(8):
:
EPS,EP1,EP2!
EPS:
迭代的精度
S=0
DOTAO=1,Q
DOI=1,N-TAO
S(TAO)=S(TAO)+X(I)*X(I+TAO)
ENDDO
S(TAO)=S(TAO)/(N-TAO)
ENDDO
S2=0
DOI=1,N
S2=S2+X(I)*X(I)
ENDDO
S2=S2/N
DOTAO=1,Q
R(TAO)=0
DOI=1,N-TAO
R(TAO)=R(TAO)+X(I)*X(I+TAO)
ENDDO
R(TAO)=R(TAO)/S2/(N-TAO)
ENDDO
THIT=0
S2B=S2
NN=0
DO
NN=NN+1
A1=1
DOI=1,Q
A1=A1+THIT(I)*THIT(I)
ENDDO
S2A=S2/A1
THITA=-R*S2/S2A
DOK=1,Q-1
DOI=1,Q-K
THITA(K)=THITA(K)+THIT(I)*THIT(K+I)
ENDDO
ENDDO
EP1=ABS(S2A-S2B)
EP2=MAXVAL(ABS(THIT-THITA))
IF(EP1THIT=THITA
S2B=S2A
PRINT*,'NN=',NN
ENDDO
OPEN(12,FILE='MAQ.DAT')
WRITE(12,*)
WRITE(12,'("S2=",D12.5)')S2
WRITE(12,'("R=",D12.5)')R
WRITE(12,'("S2A=",E12.5)')S2A
WRITE(12,'("THITA=",D12.5)')THITA
CLOSE(12)
END
6.3.6例
计算北京1951年——1980年1月的平均气温2阶、3阶滑动平均模型的系数(同时也算出了12月、2月的结果)
PROGRAMMAIN
INTEGER,PARAMETER:
:
N=30
INTEGER,PARAMETER:
:
Q=2
REAL(8),DIMENSION(N):
:
X
REAL(8),PARAMETER:
:
EPS=1.0E-5
REAL(8):
:
XV!
X的平均值
OPEN(10,FILE='BEIJING.DAT')
READ(10,*)X
CLOSE(10)
XV=0
DOI=1,N
XV=XV+X(I)
ENDDO
XV=XV/N
X=X-XV
CALLMAQ(X,N,Q,EPS)
END
计算结果为:
2阶:
S2=.11905D+01
R=-.82189D-01.65269D-01
S2A=.11782E+01
THITA=.77908D-01-.65949D-01
滑动平均模型为:
3阶:
S2=.11905D+01
R=-.82189D-01.65269D-01.23275D-01
S2A=.11770E+01
THITA=.79343D-01-.67884D-01-.23542D-01
滑动平均模型为:
6.3.7附注
6.4自回归滑动平均模型(ARMA)
6.4.1功能
求出(p,q)阶自回归—滑动平均方程的系数,从而得到(p,q)阶自回归—滑动平均方程。
6.4.2方法说明
6.4.3子程序语句
SUBROUTINEARMA(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子程序
SUBROUTINEARMA(X,N,P,Q,M,FAI,THITA,EPS)
INTEGER:
:
TAO!
落后时间
INTEGER:
:
P!
自回归阶数
INTEGER:
:
Q!
滑动平均阶数
INTEGER:
:
M!
M=P+Q
REAL(8),DIMENSION(N):
:
X!
输入序列
REAL(8),DIMENSION(0:
P):
:
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!
EPS:
迭代的精度
S=0
DOTAO=0,M
DOI=1,N-TAO
S(TAO)=S(TAO)+X(I)*X(I+TAO)
ENDDO
S(TAO)=S(TAO)/(N-TAO)
ENDDO
DOI=1,P
DOJ=1,P
A(I,J)=S(ABS(Q+I-J))
ENDDO
B(I)=S(Q+I)
ENDDO
CALLGASJDN(A,B,P)
FAI(1:
P)=B(1:
P)
FAI(0)=-1
A1=0
DOI=0,P
A1=A1+FAI(I)*FAI(I)
ENDDO
DOK=0,Q
A2=0
DOI=1,P
A3=0
DOJ=0,P-I
A3=A3+FAI(J)*FAI(J+I)
ENDDO
A2=A2+A3*(S(K+I)+S(ABS(K-I)))
ENDDO
SC(K)=A1*S(K)+A2
ENDDO
S2B=0
THIT=0
NN=0
DO!
迭代
NN=NN+1
WRITE(*,'("NN=",I3)')NN
A1=1
DOI=1,Q
A1=A1+THIT(I)*THIT(I)
ENDDO
S2A=SC(0)/A1
DOK=1,Q
THITA(K)=-SC(K)/S2A
DOI=1,Q-K
THITA(K)=THITA(K)+THIT(I)*THIT(K+I)
ENDDO
ENDDO
EP1=ABS(S2A-S2B)
EP2=MAXVAL(ABS(THIT-THITA))
WRITE(*,*)S2A,EP1,EP2
IF(EP1THIT=THITA
S2B=S2A
ENDDO
END
!
全选主元高斯——约当法(Gauss-Jordan)求解n阶线性代数方程组
SUBROUTINEGASJDN(A,B,N)
REAL(8),DIMENSION(N,N):
:
A
REAL(8),DIMENSION(N):
:
B
REAL(8),DIMENSION(N):
:
JA
REAL(8):
:
DMAX,DD
LL=1
DOK=1,N
DMAX=0
DOI=K,N
DOJ=K,N
IF(ABS(A(I,J))>DMAX)THEN
DMAX=ABS(A(I,J))
JA(K)=J
IA=I
ENDIF
ENDDO
ENDDO
IF(DMAX+1==1)THEN
WRITE(*,'("主元为0,求解失败")')
LL=0
RETURN
ENDIF
DOJ=K,N
DD=A(K,J)
A(K,J)=A(IA,J)
A(IA,J)=DD
ENDDO
DD=B(K)
B(K)=B(IA)
B(IA)=DD
DOI=1,N
DD=A(I,K)
A(I,K)=A(I,JA(K))
A(I,JA(K))=DD
ENDDO
DOJ=K+1,N
A(K,J)=A(K,J)/A(K,K)
ENDDO
B(K)=B(K)/A(K,K)
DOJ=K+1,N
DOI=1,N
IF(I/=K)A(I,J)=A(I,J)-A(I,K)*A(K,J)
ENDDO
ENDDO
DOI=1,N
IF(I/=K)THEN
B(I)=B(I)-A(I,K)*B(K)
ENDIF
ENDDO
ENDDO
DOK=N,1,-1
DD=B(K)
B(K)=B(JA(K))
B(JA(K))=DD
ENDDO
END
6.4.6例
以7.3.6中资料为例,计算北京1951年——1980年1月的平均气温2阶字回归和1阶滑动平均模型的系数。
PROGRAMARMAPQ
INTEGER,PARAMETER:
:
N=30
INTEGER,PARAMETER:
:
P=2,Q=1,M=P+Q!
P自回归阶数,Q滑动平均阶数
REAL(8),DIMENSION(N):
:
X
REAL(8),DIMENSION(0:
P):
:
FAI
REAL(8),DIMENSION(Q):
:
THITA!
滑动系数
REAL(8):
:
XV!
X的平均值
REAL(8):
:
EPS
EPS=1.E-4
OPEN(10,FILE='BEIJING.DAT')
READ(10,*)X
CLOSE(10)
XV=0
DOI=1,N
XV=XV+X(I)
ENDDO
XV=XV/N
X=X-XV
CALLARMA(X,N,P,Q,M,FAI,THITA,EPS)
OPEN(12,FILE='ARMA.DAT')
WRITE(12,'(2X,"XV=",F8.4)')XV
DOI=1,P
WRITE(12,'("FAI(",I2,")=",F8.4)')I,FAI(I)
ENDDO
DOI=1,Q
WRITE(12,'("THITA(",I2,")=",F8.4)')I,THITA(I)
ENDDO
CLOSE(12)
END
计算结果:
XV=-4.5467
FAI
(1)=.4894
FAI
(2)=.1055
THITA
(1)=.5697
因此,自回归滑动平均方程为:
注意:
上式得到的结果加上平均值XV即为实际预报值。