C PROGRAM 7.3 MARLSQ C C ... MAR model fitting (least squares method) ... C C Inputs: C MT: Input device C LAG: Heighest AR order C Inputs for READMD C TITLE: Tile of the data set C N: Data length C L: Dimension C IFM: Control variable for reading order C FORM: Reading format C Y(N,J): ID-dimensional time series C PARAMETER( MJ=1000,MJ1=200,MJ2=5,MJ3=20,MJ4=MJ2*(MJ3+1),IDEV=1 ) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(MJ,MJ2), X(MJ1,MJ4), D(MJ1) DIMENSION A(MJ2,MJ2,MJ3) , V(MJ2,MJ2) EXTERNAL SETMAR C IPR = 3 READ( 5,* ) LAG C C ... Read multi-variate time series ... C CALL READMD( IDEV,MJ,Y,N,L ) N0 = 0 NMK = N - LAG C C ... Housholder reduction ... C CALL MREDCT( SETMAR,Y,D,NMK,N0,LAG,L,MJ,MJ1,X ) C C ... Least squares method for MAR model ... C CALL MARFIT( X,D,NMK,L,LAG,MJ1,MJ2,MJ3,IPR,A,V,LMAX,AIC ) C C ... Print out fitted model ... C CALL PRMAR2( A,V,L,LMAX,LAG,N,MJ2,MJ3 ) C WRITE(2,*) LMAX C DO 97 I=1,L C 97 WRITE(2,600) (V(I,J),J=1,L) C DO 98 II=1,LMAX C DO 98 I=1,L C 98 WRITE(2,600) (A(I,J,II),J=1,L) C 600 FORMAT( 3D15.7 ) C STOP E N D SUBROUTINE MARFIT( X,D,N,ID,M,MJ1,MJ2,MJ3,IPR,B,E,LMAX,AICS ) C C ... Multivariate AR model fitting (Householder method) ... C C Inputs: C X: Householder reduced form C D: Working area C N: Data length C ID: Dimension of the series C M: Highest AR order C IPR: Print out control C MJ1, MJ2, MJ3: Adjustable dimensions C Outputs: C B: AR-coefficient matrices C E: innovation covariance matrix C LMAX: Order of the MAICE model C AICS: Total AIC of the model C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), D(1), E(MJ2,1), B(MJ2,MJ2,MJ3) DIMENSION A(100), AIC(51), SD(51), EX(10) DIMENSION IND(100), JND(100) DATA PI2/6.28318531D0/ C MD = M*ID MD2 = (M+1)*ID AICSUM = 0.D0 LMAX = 0 DO 20 I=1,MJ2 DO 20 J=1,MJ2 DO 10 K=1,MJ3 10 B(I,J,K) = 0.D0 20 E(I,J) = 0.D0 DO 30 I=1,MD2 30 JND(I) = I C DO 200 II=1,ID C JJ = II - 1 KK = MD + JJ IF( II .GT. 1 ) THEN C C ... Addition of instantaneous variable ... C DO 110 I=1,MD2 J = JND(I) 110 IND(J) = I JND(1) = 1 DO 120 I=1,JJ 120 JND(I) = MD+I DO 130 I=1,MD 130 JND(I+JJ) = I C I1 = JJ + 1 DO 140 I=JJ+1,ID J = MD + I 140 JND(J) = J C CALL HUSHL1( X,D,MJ1,MD2,MD2,1,IND,JND ) END IF C C ... AIC's of the models for II-th variable ... C DO 160 I=1,M+1 K = (I-1)*ID + JJ SUM = 0.0D0 DO 150 IJ=K+1,KK+1 150 SUM = SUM + X(IJ,KK+1)**2 SD(I) = SUM/N 160 AIC(I) = N*DLOG( PI2*SD(I) ) + N + 2*(K+1) C C ... Order determination by AIC ... C CALL MAICE( AIC,SD,M,IPR-2,AICMIN,SDMIN,IMIN ) C K0 = IMIN*ID + JJ CALL SRCOEF( X,K0,KK,N,MJ1,JND,A,EX(II),AICX ) DO 170 I=1,II-1 170 E(II,I) = -A(I) E(II,II) = 1.0D0 DO 180 IJ=1,IMIN DO 180 J=1,ID 180 B(II,J,IJ) = A((IJ-1)*ID+J+II-1) AICSUM = AICSUM + AICMIN IF( IMIN.GT.LMAX ) LMAX = IMIN C 200 CONTINUE C IF(IPR.GE.1) WRITE( 6,3 ) CALL MCOEF( B,E,EX,ID,LMAX,MJ2,MJ3 ) IF(IPR.GE.1) WRITE( 6,611 ) AICSUM AICS = AICSUM C RETURN 3 FORMAT( 1H0,132(1H-) ) 15 FORMAT( 1H0,'***** FINAL ESTIMATE *****' ) 611 FORMAT( 1H0,'AIC =',F15.3 ) END SUBROUTINE MCOEF( B,E,EX,ID,LMAX,MJ2,MJ3 ) C C ... AR coefficient matrices and innovation covariance matrix ... C C Inputs: C B: Regression coefficient matrix C E: B(0) C EX: Residual variances of orthogonarized model C ID: Dimension C LMAX: AR-order C MJ2: Adjustable dimension of B and E C MJ3: djustable dimension of B C Outputs: C B: AR-coefficient matrix C E: Innovation covaraince matrix C IMPLICIT REAL*8 ( A-H,O-Z ) DIMENSION B(MJ2,MJ2,MJ3), E(MJ2,1), EX(1), EE(24,24) MJ5 = 24 C CALL TRIINV( E,ID,MJ2,MJ5,EE ) C DO 30 II=1,LMAX DO 20 I=1,ID DO 20 J=1,ID SUM = 0.D0 DO 10 JJ=1,I 10 SUM = SUM + EE(I,JJ)*B(JJ,J,II) 20 E(I,J) = SUM DO 30 I=1,ID DO 30 J=1,ID 30 B(I,J,II) = E(I,J) DO 50 I=1,ID DO 50 J=1,I SUM = 0.D0 DO 40 II=1,J 40 SUM = SUM + EE(I,II)*EE(J,II)*EX(II) E(I,J) = SUM 50 E(J,I) = SUM DO 99 I=1,ID 99 WRITE(6,*) (E(I,J),J=1,ID) C RETURN E N D SUBROUTINE PRMAR2( A,V,L,M,LAG,N,MJ2,MJ3 ) C C ... Print out AR coefficient matrices and innovation covariance ... C C Inputs: C A: AR coefficient matrix C V: Innovation covariance matrix C L: Dimension C M: AR order C LAG: Highest order C N: Data length C MJ2 and MJ3: Adjustable dimension C IMPLICIT REAL*8 ( A-H,O-Z ) CHARACTER*72 TITLE DIMENSION A(MJ2,MJ2,MJ3), V(MJ2,1) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) N, L, LAG, M WRITE(6,630) DO 10 I=1,L 10 WRITE(6,640) (V(I,J),J=1,L) WRITE(6,650) DO 20 II=1,M WRITE(6,660) II DO 20 I=1,L 20 WRITE(6,670) (A(I,J,II),J=1,L) C RETURN 600 FORMAT( 1H ,'PROGRAM 7.3: MAR MODEL (HOUSEHOLDER METHOD)' ) 610 FORMAT( 1H ,A72 ) 620 FORMAT( 1H ,'N(DATA LENGTH =',I6,5X,'L(DIMENSION) =',I3,5X * 'LAG(HEIGHEST ORDER) =',I3,5X,'M(MAICE ORDER) =',I3 ) 630 FORMAT( 1H ,'INNOVATION COVARIANCE MATRIX' ) 640 FORMAT( 1H ,5D15.7 ) 650 FORMAT( 1H ,'AR COEFFICIENT MATRICES' ) 660 FORMAT( 1H ,'M =',I3 ) 670 FORMAT( 1H ,5F12.6 ) END SUBROUTINE MREDCT( MSETX,Z,D,NMK,N0,LAG,ID,MJ,MJ1,X ) C C ... This subroutine computes Householder reduced form of the C data matrix for fitting MAR model ... C C Inputs: C MSETX: Function name C Z: Multivariate time series C D: Working area C NMK: Actual data length C N0: Origin of new data C LAG: Highest order of the AR model C ID: Dimension of observations C MJ,MJ1: Adjustable dimension C OUTPUT: C X: Householder reduced form (upper triangular form) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(MJ1,1) , D(1), Z(MJ,1) C L = MIN0( NMK,MJ1 ) KD1 = (LAG+1)*ID N1 = L C CALL MSETX( Z,N0,L,LAG,ID,MJ,MJ1,0,X ) CALL HUSHLD( X,D,MJ1,L,KD1 ) C IF( N1 .GE. NMK ) RETURN 100 L = MIN0( NMK-N1,MJ1-KD1 ) LK = L + KD1 N2 = N0 + N1 C CALL MSETX( Z,N2,L,LAG,ID,MJ,MJ1,1,X ) CALL HUSHLD( X,D,MJ1,LK,KD1 ) N1 = N1 + L IF( N1 .LT. NMK ) GO TO 100 C RETURN E N D SUBROUTINE SETMAR( Z,N0,L,LAG,ID,MJ,MJ1,JSW,X ) C C ... Make matrix X for multi-variate AR model ... C C Inputs: C Z: Multivariate time series C N0: Origin of new data C L: Number of new observations C LAG: Highest order of AR model C ID: Dimension of time series C MJ,MJ1: Adjustable dimension C JSW: =0 To make initail L*(LAG+1) data matrix C =1 To augment original (LAG+1)*(LAG+1) matrix X C by L*(LAG+1) data matrix of additional obs. C Outputs C X: L*(LAG+1) matrix if JSW = 0 C (LAG+1+L)*(LAG+1) matrix if JSW = 1 C REAL*8 X(MJ1,1), Z(MJ,1) C KD = LAG*ID KD1 = (LAG+1)*ID I0 = 0 IF( JSW.EQ.1 ) I0 = KD1 C DO 30 II=1,L I1 = N0 + LAG + II I2 = I0 + II DO 10 J=1,ID J2 = KD + J 10 X(I2,J2) = Z(I1,J) DO 30 JJ=1,LAG I1 = I1 - 1 J1 = (JJ-1)*ID DO 20 J=1,ID J2 = J1 + J 20 X(I2,J2) = Z(I1,J) 30 CONTINUE C RETURN E N D SUBROUTINE TRIINV( X,M,MJ,MJ1,Y ) C C ... Lower triangular matrix inversion ... C C Inputs: C X: Triangular matrix with unit diagonal elements C M: Dimension of matrix X C MJ: Adjustable dimension of X C MJ1: Adjustable dimension of Y C Output: C Y: Inverse of X C IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION X(MJ,1) , Y(MJ1,1) DO 10 I=1,M-1 DO 10 J=I,M 10 Y(I,J) = 0.0D0 DO 20 I=1,M 20 Y(I,I) = 1.0D0 DO 40 J=1,M-1 DO 40 I=J+1,M SUM = 0.0D0 DO 30 II=1,I-J JJ = II + J - 1 30 SUM = SUM + X(I,JJ) * Y(JJ,J) 40 Y(I,J) = -SUM RETURN E N D SUBROUTINE MAICE( AIC,SD,K,ISW,AICM,SDM,IMIN ) C C INPUTS: C AIC: VECTOR OF AIC'S C SD: VECTOR OF INNOVATION VARIANCES C K: UPPER LIMIT OF THE ORDER C ISW: =0 OUTPUTS ARE SUPPRESSED C >0 AIC'S ARE DIPLAIED C OUTPUTS: C AICM: MINIMUM AIC C SDM: MAICE INNOVATION VARIANCE C IMIN: MAICE ORDER C IMPLICIT REAL*8(A-H,O-Z) DIMENSION AIC(1), SD(1) C C SEARCH FOR THE MINIMUM OF AIC(I) C IMIN = 0 SDM = SD(1) AICM = AIC(1) DO 20 I=1,K IF( AIC(I+1).LT.AICM ) THEN IMIN = I SDM = SD(I+1) AICM = AIC(I+1) END IF 20 CONTINUE IF( ISW .LE. 0 ) RETURN C C DISPLAY OF AIC'S C WRITE( 6,5 ) DO 30 I=1,K+1 II = I - 1 DIC = AIC(I) - AICM 30 WRITE( 6,6 ) II, SD(I), AIC(I), DIC WRITE( 6,7 ) IMIN, AICM, SDM C RETURN 5 FORMAT( 1H ,4X,'M',9X,'SIG2(M)',13X,'AIC(M)',7X,'AIC(M)-AICMIN' ) 6 FORMAT( 1H ,I5,D20.10,2F16.3 ) 7 FORMAT( 1H0,'MAICE ORDER =',I3,5X,'AIC =',F15.3,5X,'SIG2 =', * D17.10 ) E N D SUBROUTINE READMD( IDEV,NMAX,Y,N,L ) C C ... THIS SUBROUTINE READS IN MULTI-VARIATE TIME SERIES ... C C INPUTS: C IDEV: INPUT FILE SPECIFICATION C NMAX: ADJUSTABLE DIMENSION OF Y(I,J) C OUTPUTS: C Y(I,J): TIME SERIES (I=1,N;J=1,L) C N: DATA LENGTH C L: DIMENSION OF THE TIME SERIES C CHARACTER*72 TITLE REAL*8 Y(NMAX,*) COMMON /CMDATA/ TITLE C C OPEN( IDEV,FILE='temp2.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N, L, IFM IF( IFM.EQ.0 ) READ( IDEV,* ) ((Y(I,J),J=1,L),I=1,N) IF( IFM.GE.1 ) READ( IDEV,* ) ((Y(I,J),I=1,N),J=1,L) C CLOSE( IDEV ) C RETURN 1 FORMAT( A72 ) E N D SUBROUTINE HUSHL1( X,D,MJ,K,L,M,IND,JND ) C C ... Householder trasnformation of the matrix X ... C C Inputs: C X(I,J): (K+1)*(K+1) matrix with the form specified by IND(I) C D(I): Working area C MJ: Adjustable dimension of X C K: Number of columns of X C M: Starting position of the Householder transformation C L: End position of the Householder transformation C IND: Specification of the present form of X C << IND(J) = I means J-th regressor is variable I>> C JND: Specification of the required form of X C << JND(I) = J means I-th regressor is variable J >> C Output: C X(I,J) Transformed matrix with the form specified by JND(I) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(MJ,1), D(40), IND(40), JND(40) TOL = 1.0D-60 C NN = 0 DO 100 II=M,L JJ = JND(II) NN = MAX0( NN,IND(JJ) ) H = 0.0D0 DO 10 I=II,NN D(I) = X(I,JJ) 10 H = H + D(I)**2 IF( H.GT.TOL ) GO TO 20 G = 0.0D0 GO TO 100 20 G = DSQRT( H ) F = X(II,JJ) IF( F.GE.0.0D0 ) G = -G D(II) = F - G H = H - F * G C DO 30 I=II+1,NN 30 X(I,JJ) = 0.D0 IF( II.EQ.K ) GO TO 100 DO 60 J1=II+1,K J = JND(J1) S = 0.0D0 DO 40 I=II,NN 40 S = S + D(I)*X(I,J) S = S / H DO 50 I=II,NN 50 X(I,J) = X(I,J) - D(I)*S 60 CONTINUE 100 X(II,JJ) = G RETURN E N D SUBROUTINE SRCOEF( X,M,K,N,MJ,JND,A,SIG2,AIC ) C C ... Subset regression coefficients, residual variance and AIC ... C C Iuputs: C X: Householder reduced form with index vector JND(I) C M: Number of regressors C K: Heighest order of the models C N: Data length C MJ: Adjustable dimension of X C JND(I): Specification of the I-th regressor C Outputs: C A(I): Regression coefficients C SIG2: Innovation variance C AIC: AIC C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ,1), A(1), JND(1) DATA PI2/6.28318531D0/ C C ... REGRESSION COEFFICIENTS... C L = JND(M) A(M) = X(M,K+1)/X(M,L) DO 10 II=1,M-1 I = M - II SUM = X(I,K+1) DO 20 J=I+1,M L = JND(J) 20 SUM = SUM - A(J)*X(I,L) L = JND(I) 10 A(I) = SUM/X(I,L) C C ... RESIDUAL VARIANCE AND AIC ... C SD = 0.0D0 DO 30 I=M+1,K+1 30 SD = SD + X(I,K+1)**2 SIG2 = SD/N AIC = N*DLOG( PI2*SIG2 ) + N + 2.0D0*(M+1) C RETURN E N D SUBROUTINE HUSHLD( X,D,MJ1,N,K ) C C ... Householder transformation ... C C Inputs: C X(I,J): Original matrix C D(I): Working area C MJ1: Adjustable dimension of X C N: Number of rows of X C K: Number of columns of X C Output: C X(I,J): Householder reduced form (upper triangular form) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), D(MJ1) TOL = 1.0D-60 C DO 100 II=1,K H = 0.0D0 DO 10 I=II,N D(I) = X(I,II) 10 H = H + D(I)**2 IF( H .GT. TOL ) GO TO 20 G = 0.0D0 GO TO 100 20 G = DSQRT( H ) F = X(II,II) IF( F .GE. 0.0D0 ) G = -G D(II) = F - G H = H - F*G C DO 30 I=II+1,N 30 X(I,II) = 0.0D0 DO 60 J=II+1,K S = 0.0D0 DO 40 I=II,N 40 S = S + D(I)*X(I,J) S = S/H DO 50 I=II,N 50 X(I,J) = X(I,J) - D(I)*S 60 CONTINUE 100 X(II,II) = G C RETURN E N D