时间序列分析doc.docx

资源描述

时间序列分析doc.docx

《时间序列分析doc.docx》由会员分享，可在线阅读，更多相关《时间序列分析doc.docx（14页珍藏版）》请在冰豆网上搜索。

时间序列分析doc.docx

时间序列分析doc

第六章时间序列分析

6．２自回归模型（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）:

REAL（4）,DIMENSION（M,M）:

FAI

REAL（4）,DIMENSION（M）:

协方差

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

END

6.2.6例

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

PROGRAMMAIN

INTEGER,PARAMETER:

N=264

INTEGER,PARAMETER:

M=12

REAL（4）,DIMENSION（N）:

REAL（4）,DIMENSION（M,M）:

FAI

REAL（4）,DIMENSION（M）:

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

（1）=.8383FAI

（1）=.6094

（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）:

REAL（8）,DIMENSION（Q）:

THITA!

滑动系数

REAL（8）,DIMENSION（Q）:

THIT!

迭代中用的滑动系数,中间变量

REAL（8）,DIMENSION（Q）:

相关系数

REAL（8）,DIMENSION（Q）:

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

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

EP1=ABS（S2A-S2B）

EP2=MAXVAL（ABS（THIT-THITA））

IF（EP1

THIT=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）:

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:

自回归阶数

INTEGER:

滑动平均阶数

INTEGER:

M=P+Q

REAL（8）,DIMENSION（N）:

输入序列

REAL（8）,DIMENSION（0:

P）:

FAI!

自回归系数

REAL（8）,DIMENSION（P,P）:

工作数组

REAL（8）,DIMENSION（P）:

工作数组

REAL（8）,DIMENSION（0:

M）:

协方差，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

EP1=ABS（S2A-S2B）

EP2=MAXVAL（ABS（THIT-THITA））

WRITE（*,*）S2A,EP1,EP2

IF（EP1

THIT=THITA

S2B=S2A

ENDDO

END

全选主元高斯——约当法（Gauss-Jordan）求解n阶线性代数方程组

SUBROUTINEGASJDN（A,B,N）

REAL（8）,DIMENSION（N,N）:

REAL（8）,DIMENSION（N）:

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

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

DOI=1,N

IF（I/=K）THEN

B（I）=B（I）-A（I,K）*B（K）

ENDIF

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）:

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即为实际预报值。

展开阅读全文