时间序列分析(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)

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

DO 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 DO

S2=0

DO I=1,N

S2=S2+X(I)*X(I)

END DO

S2=S2/N

DO TAO=1,M

R(TAO)=0

DO I=1,N-TAO

R(TAO)=R(TAO)+X(I)*X(I+TAO)/S2

END DO

R(TAO)=R(TAO)/(N-TAO)

END DO 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)

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)-FAI(J,J)*FAI(J-K,J-1)

END DO

END DO

END

6.2.6例

以某海区的22年的逐月气温为例，计算出自回归系数，并给出自回归方程。

PROGRAM MAIN

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')

DO I=1,N

READ(10,'(F8.2)')X(I)

END DO

CLOSE(10)

XV=0

DO I=1,N

XV=XV+X(I)

END DO

XV=XV/N

X=X-XV

CALL ARP(X,N,M,R,FAI)

OPEN(12,FILE='ARP.DAT')

WRITE(12,'(2X,"XV=",F8.4)')XV

DO I=1,M

WRITE(12,'("R(",I2,")=",F8.4," FAI(",I2,")=",F8.4)')I,R(I),I,FAI(I,M)

END DO

CLOSE(12)

END 输出结果为：

XV= 22.5718

R( 1)= .8383 FAI( 1)= .6094

R( 2)= .4648 FAI( 2)= -.1669

R( 3)= -.0148 FAI( 3)= -.0701

R( 4)= -.4776 FAI( 4)= -.0564

R( 5)= -.8080 FAI( 5)= -.1197

R( 6)= -.9222 FAI( 6)= .0477

R( 7)= -.8019 FAI( 7)= -.0471

R( 8)= -.4747 FAI( 8)= -.1702

R( 9)= -.0108 FAI( 9)= .0053

R(10)= .4665 FAI(10)= .0977

R(11)= .8211 FAI(11)= .1246

R(12)= .9508 FAI(12)= .1798

从而得到自回归方程为：

12t11t10t9t8t7t6t5t4t3t2t1ttx1798.x1246.x0977.x0053.x1702.x0471.x0477.x1197.x0564.x0701.x1669.x6094.x

注意：以上是距平值，加上平均值即为实际值。

6．3滑动平均模型（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)::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

DO TAO=1,Q DO I=1,N-TAO

S(TAO)=S(TAO)+X(I)*X(I+TAO)

END DO

S(TAO)=S(TAO)/(N-TAO)

END DO

S2=0

DO I=1,N

S2=S2+X(I)*X(I)

END DO

S2=S2/N

DO 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 DO

THIT=0

S2B=S2

NN=0

DO

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=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

THIT=THITA

S2B=S2A

PRINT*,'NN=',NN

END DO

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月的结果）

PROGRAM MAIN

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

DO I=1,N

XV=XV+X(I)

END DO

XV=XV/N

X=X-XV

CALL MAQ(X,N,Q,EPS)

END

计算结果为：

2阶：

S2= .11905D+01

R= -.82189D-01 .65269D-01

S2A= .11782E+01

THITA= .77908D-01 -.65949D-01

滑动平均模型为：

2t1ttta065949.0a0779083.0aX

3阶：

S2= .11905D+01

R= -.82189D-01 .65269D-01 .23275D-01

S2A= .11770E+01

THITA= .79343D-01 -.67884D-01 -.23542D-01

滑动平均模型为：

3t2t1ttta023542.0a067884.0a079343.0aX

6.3.7附注

6．4自回归滑动平均模型（ARMA）