C PROGRAM 7.2 MARFIT C C ... This program fits multivariate AR model ... C C The following inputs are required in the subroutine READMD. C TITLE: title of the data C N: data length C ID: dimension of the observation C IFM: = 1 ((Y(I,J),J=1,ID),I=1,N) C = 2 ((Y(I,J),I=1,N),J=1,ID) C Y(I,J): multi-variate time series C Parameters: C MJ: adjustable dimension of Y; (MJ.GE.N) C MJ1: adjustable dimension of Y, C, R; (MJ1.GE.ID) C LAG: maximum lag of the cross-covariance function C PARAMETER( MJ0=500,MJ=4,MJ1=10,LAG=MJ1,IDEV=1 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(MJ0,MJ), OUTMIN(10), OUTMAX(10) DIMENSION C(0:MJ1,MJ,MJ), R(0:MJ1,MJ,MJ), YMEAN(10) DIMENSION AMIN(MJ1,MJ,MJ), VMIN(MJ,MJ), AIC(0:MJ1) DATA OUTMIN/10*-1.0D30/, OUTMAX/10*1.0D30/ C C ... read in multivariate time series ... C CALL READMD( IDEV,MJ0,Y,N,L ) C C ... cross-covariance function ... C CALL CRSCOR( Y,N,L,LAG,MJ0,MJ,OUTMIN,OUTMAX,C,R,YMEAN ) C C ... Yule-Walker method for MAR model ... C CALL MYULE( L,LAG,N,C,MJ,MJ1,AMIN,VMIN,MMIN,AIC ) C C ... print out estimated model ... C CALL PRMAR( L,N,LAG,MJ,MJ1,AMIN,VMIN,MMIN,AIC ) STOP E N D SUBROUTINE MYULE( L,LAG,N,C,MJ,MJ1,AMIN,VMIN,MMIN,AIC ) C C ... Yule-Walker method for fitting multivariate AR model ... C C Inputs: C L: Dimension of the time series C N: Data length C LAG: Highest order of fitted AR models C C(K,I,J): Crosscovariance function C MJ, MJ1: Adjustable dimensions C Outputs: C AMIN(K,I,J): AR coefficient of the AIC best model C VMIN(I,J): Innovation covariance matrix of AIC best model C MMIN: MAICE order C AIC(I): AIC's of the AR models C IMPLICIT REAL*8(A-H,O-Z) DIMENSION C(0:MJ1,MJ,MJ) DIMENSION AMIN(MJ1,MJ,MJ), VMIN(MJ,MJ), AIC(0:MJ1) DIMENSION A(20,10,10), A0(20,10,10), B(20,10,10), B0(20,10,10) DIMENSION V(10,10), U(10,10), W(10,10) DATA PI2 /6.28318531D0/ MJ2 = 10 C DO 10 I=1,L DO 10 J=1,L V(I,J) = C(0,I,J) U(I,J) = C(0,I,J) 10 VMIN(I,J) = V(I,J) CALL INVDET( U,UDET,L,MJ2 ) CALL INVDET( V,VDET,L,MJ2 ) AIC(0) = N*L*(DLOG(PI2) + 1) + N*DLOG(VDET) + L*(L+1) AICMIN = AIC(0) MMIN = 0 C DO 300 M=1,LAG C DO 110 I=1,L DO 110 J=1,L SUM = C(M,I,J) DO 100 K=1,M-1 DO 100 IJ=1,L 100 SUM = SUM - A0(K,I,IJ)*C(M-K,IJ,J) 110 W(I,J) = SUM C DO 130 I=1,L DO 130 J=1,L SUM1 = 0.0D0 SUM2 = 0.0D0 DO 120 IJ=1,L SUM1 = SUM1 + W(I,IJ)*U(IJ,J) 120 SUM2 = SUM2 + W(IJ,I)*V(IJ,J) A(M,I,J) = SUM1 130 B(M,I,J) = SUM2 C DO 150 K=1,M-1 DO 150 I=1,L DO 150 J=1,L SUM1 = A0(K,I,J) SUM2 = B0(K,I,J) DO 140 IJ=1,L SUM1 = SUM1 - A(M,I,IJ)*B0(M-K,IJ,J) 140 SUM2 = SUM2 - B(M,I,IJ)*A0(M-K,IJ,J) A(K,I,J) = SUM1 150 B(K,I,J) = SUM2 C DO 170 I=1,L DO 170 J=1,L SUM1 = C(0,I,J) SUM2 = C(0,I,J) DO 160 K=1,M DO 160 IJ=1,L SUM1 = SUM1 - A(K,I,IJ)*C(K,J,IJ) 160 SUM2 = SUM2 - B(K,I,IJ)*C(K,IJ,J) V(I,J) = SUM1 W(I,J) = SUM1 170 U(I,J) = SUM2 C DO 180 K=1,M DO 180 J=1,L DO 180 I=1,L A0(K,I,J) = A(K,I,J) 180 B0(K,I,J) = B(K,I,J) C CALL INVDET( U,UDET,L,MJ2 ) CALL INVDET( V,VDET,L,MJ2 ) AIC(M) = N*L*(DLOG(PI2) + 1) + N*DLOG(VDET) + L*(L+1) + 2*L*L*M C IF( AIC(M).LT.AICMIN ) THEN AICMIN = AIC(M) MMIN = M DO 200 I=1,L DO 200 J=1,L VMIN(I,J) = W(I,J) DO 200 K=1,M 200 AMIN(K,I,J) = A(K,I,J) END IF 300 CONTINUE C RETURN E N D SUBROUTINE PRMAR( L,N,LAG,MJ,MJ1,AMIN,VMIN,MMIN,AIC ) C C ... Print out AR model and AIC's ... C C Inputs: C L: Dimension of the time series C N: Data length C LAG: Heighest order of fitted AR models C AMIN(K,I,J): AR coefficient of the AIC best model C VMIN(I,J): Innovation covariance matrix of AIC best model C MMIN: MAICE order C AIC(I): AIC's of the AR models C IMPLICIT REAL*8(A-H,O-Z) CHARACTER*72 TITLE DIMENSION AMIN(MJ1,MJ,MJ), VMIN(MJ,MJ), AIC(0:LAG) COMMON /CMDATA/ TITLE C WRITE(6,680) TITLE WRITE(6,690) N, LAG WRITE(6,600) WRITE(6,610) DO 10 I=0,LAG 10 WRITE(6,620) I, AIC(I) WRITE(6,625) MMIN C WRITE(6,630) DO 20 I=1,L 20 WRITE(6,640) (VMIN(I,J),J=1,L) C WRITE(6,650) DO 30 K=1,MMIN WRITE(6,660) K DO 30 I=1,L 30 WRITE(6,670) (AMIN(K,I,J),J=1,L) C RETURN 600 FORMAT( 1H ,'** YULE WALKER ESTIMATE **' ) 610 FORMAT( 1H ,'ORDER',10X,'AIC' ) 620 FORMAT( 1H ,I4,F15.4 ) 625 FORMAT( 1H ,'** MAICE ORDER ** MMIN =',I3 ) 630 FORMAT( 1H ,'** INNOVATION COVARIANCE MATRIX **' ) 640 FORMAT( 1H ,5D15.6 ) 650 FORMAT( 1H ,'** AR COEFFICIENT MATRICES **' ) 660 FORMAT( 1H ,'K =',I2 ) 670 FORMAT( 1H ,5F12.6 ) 680 FORMAT( 1H ,A72 ) 690 FORMAT( 1H ,'N(DATA LENGTH) =',I6,5X,'LAG(HEIGHEST ORDER) =',I3 ) E N D SUBROUTINE INVDET( X,DET,M,MJ ) C C ... Inverse and determinant of matrix X ... C C Inputs: C X: M*M square matrix C M: Dimension of X C MJ: Adjustable dimension of X C Outputs: C X: Inverse of X C DET: Determinant of X C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(MJ,MJ), IND(100) C DET = 1.0D0 DO 60 L=1,M XMAX = 0.1D-10 IMAX = 0 DO 10 I=L,M IF( DABS( X(I,L) ).GT.DABS( XMAX ) ) THEN XMAX = X(I,L) IMAX = I END IF 10 CONTINUE IND(L) = IMAX IF( IMAX.NE.L ) THEN IF( IMAX.LE.0 ) THEN DET = 0.0D0 RETURN END IF C C Row interchange DO 20 J=1,M XTEMP = X(IMAX,J) X(IMAX,J) = X(L,J) 20 X(L,J) = XTEMP DET = -DET END IF DET = DET*XMAX X(L,L) = 1.0D0 DO 30 J=1,M 30 X(L,J) = X(L,J)/XMAX DO 50 I=1,M IF(I.NE.L) THEN XTEMP =X(I,L) X(I,L) = 0.0D0 DO 40 J=1,M 40 X(I,J) = X(I,J) - XTEMP*X(L,J) END IF 50 CONTINUE 60 CONTINUE C C Column interchange IF( M.GT.1 ) THEN DO 110 J=1,M-1 JJ = IND(M-J) IF( JJ.NE.M-J ) THEN DO 100 I=1,M XTEMP = X(I,JJ) X(I,JJ) = X(I,M-J) 100 X(I,M-J) = XTEMP END IF 110 CONTINUE END IF RETURN E N D SUBROUTINE CRSCOR( Y,N,ID,LAG,MJ,MJ1,OUTMIN,OUTMAX,C,R,YMEAN ) C C ... cross correlation function computation ... C C Inputs: C Y(I,J): multi-variate time series C N: data length C ID: dimension of the observation C LAG: maximum lag of cross-covariance C MJ: adjustable dimension (MJ.GE.N) C MJ1: adjustable dimension (MJ1.GE.ID) C OUTMIN: bound for outliers in low side C OUTMAX: bound for outliers in high side C Outputs: C C: sample cross-covariance function C R: sample cross-correlation function C YMEAN: sample mean vector C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(MJ,MJ1), OUTMIN(MJ1), OUTMAX(MJ1) DIMENSION C(0:LAG,MJ1,MJ1), R(0:LAG,MJ1,MJ1) DIMENSION YMEAN(MJ1), NSUM(10) C DO 10 J=1,ID 10 CALL MEAN( Y(1,J),N,OUTMIN(J),OUTMAX(J),NSUM(J),YMEAN(J) ) C DO 30 I=1,ID DO 30 J=1,ID WRITE(6,*) I,J DO 30 L=0,LAG SUM = 0.0D0 NNN = 0 DO 20 II=L+1,N IF( Y(II,I).GT.OUTMIN(I).AND.Y(II,I).LT.OUTMAX(I) ) THEN IF( Y(II-L,J).GT.OUTMIN(J).AND.Y(II-L,J).LT.OUTMAX(J) ) THEN SUM = SUM + (Y(II,I)-YMEAN(I))*(Y(II-L,J)-YMEAN(J)) NNN = NNN + 1 END IF END IF 20 CONTINUE 30 C(L,I,J)=SUM/DSQRT( DBLE( NSUM(I)*NSUM(J) ) ) C DO 40 I=1,ID DO 40 J=1,ID DO 40 L=0,LAG 40 R(L,I,J) = C(L,I,J)/DSQRT(C(0,I,I)*C(0,J,J)) RETURN C 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 MEAN( Y,N,OUTMIN,OUTMAX,NSUM,YMEAN ) C C ... THIS SUBROUTINE COMPUTES SAMPLE MEAN ... C C INPUTS: C Y(I): TIME SERIES C N: DATA LENGTH C OUTMIN: BOUND FOR OUTLIERS IN LOW SIDE C OUTMAX: BOUND FOR OUTLIERS IN HIGH SIDE C OUTPUTS: C NSUM: NUMBER OF NON-OUTLIER OBSERVAIONS C YMEAN: SAMPLE MEAN C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N) C NSUM = 0 YSUM = 0.0D0 DO 10 I=1,N IF( Y(I) .GT.OUTMIN .AND. Y(I).LT.OUTMAX ) THEN NSUM = NSUM + 1 YSUM = YSUM + Y(I) END IF 10 CONTINUE YMEAN = YSUM/NSUM C RETURN E N D