C PROGRAM 11.1 POLREG C C ... Polynomial regression model ... C C INPUT: C K: Order of polynomial regression C PARATERS: C MJ: Adjustable dimension of Y (MJ.GE.N) C MJ1: Adjustable dimension of X and D C MJ2: Adjustable dimension of A, SIG2 and AIC C IDEV: Input device C LAG: Number of sine and cosine terms C @TEST.PN51: C PARAMETER( MJ=1000,MJ1=200,MJ2=21,ISW=2,IDEV=1 ) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(MJ), AIC(0:MJ2) DIMENSION X(MJ1,MJ2), D(MJ1), A(MJ2,MJ2), SIG2(0:MJ2) EXTERNAL SETXPL C READ( 5,* ) K CALL READTS( IDEV,Y,N ) C CALL REDUCT( SETXPL,Y,D,N,0,K+1,MJ1,X ) C CALL REGRES( X,K+1,N,MJ1,MJ2,A,SIG2,AIC,IMIN ) C CALL PRPOL( N,K+1,MJ2,A,SIG2,AIC ) C CALL PTPOL( Y,N,A(1,IMIN),IMIN ) C STOP E N D SUBROUTINE SETXPL( Z,N0,L,K,MJ1,JSW,X ) C C ... Data matrix for polynomial regression ... C C Inputs: C Z(I): Data vector C N0: Origin of the current observations C L: Number of current observations C K: Number of regressors C (Order of polynomial = K-1) C MJ1: Adjustable dimension of X C JSW=0: Make initial data matrix C =1: Apend L*(K+1) data matrix below the triangular one C Output: C X(I,J): Data matrix C REAL*8 X(MJ1,1), Z(1) C I0 = 0 IF( JSW .EQ. 1 ) I0 = K+1 DO 10 I=1,L II = I + I0 X(II,K+1) = Z(N0+I) X(II,1) = 1.0D0 SUM = 1.0D0 DO 10 J=1,K-1 SUM = SUM*(N0+I) 10 X(II,J+1) = SUM C RETURN C E N D SUBROUTINE PRPOL( N,K,MJ,A,SIG2,AIC ) C C ... This subroutine prints out fitted polynomial ... C C Inputs: C N: Data length C K: Heighest order C MJ: Adjustable dimension of A, SIG2 and AIC C A(I,M): Regression coefficients of the model with order M C SIG2(M): Residual varance of the nodel with oder M C AIC(M): AIC of the nodel with oder M C IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION A(MJ,MJ), SIG2(0:MJ), AIC(0:MJ) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) N, TITLE IMIN = 0 AICM = AIC(0) DO 10 I=1,K IF( AIC(I).LT.AICM ) THEN IMIN = I AICM = AIC(I) END IF 10 CONTINUE C WRITE( 6,620 ) DO 20 I=0,K DIC = AIC(I) - AICM 20 WRITE(6,630) I, SIG2(I), AIC(I), DIC WRITE(6,640) IMIN, AICM, SIG2(IMIN) C WRITE(6,650 ) DO 30 M=1,K WRITE(6,660) M 30 WRITE(6,670) (A(I,M),I=1,M) C RETURN 600 FORMAT( 1H ,'PROGRAM 5.1: REGRESSION MODEL' ) 610 FORMAT( 1H ,'DATA:',2X,'N =',I5,5X,A72 ) 620 FORMAT( 1H ,'ORDER',5X,'SIG2',10X,'AIC',6X,'AIC-AICMIN' ) 630 FORMAT( 1H ,I3,D14.6,2F12.3 ) 640 FORMAT( 1H0,'IMIN =',I3,5X,'AIC =',F12.3,5X,'SIG2 =',D14.6 ) 650 FORMAT( 1H ,'REGRESSION COEFFICIENTS' ) 660 FORMAT( 1H ,'M =',I3 ) 670 FORMAT( 1H ,10D13.5 ) E N D SUBROUTINE REGRES( X,K,N,MJ1,MJ2,A,SIG2,AIC,IMIN ) C C ... Regression model fitting ... C ... Order of the model is selected by AIC ... C C Inputs: C X: Householder reduced form (upper triangular matrix) C K: Maximum number of regressors C N: Number of data C MJ1: Adjustable dimension of X C MJ2: Adjustable dimension of A, ISG2 and AIC C Outputs: C A(I,M): Regression coefficients of the model with order M C SIG2: Residual variances C AIC: AIC's C IMIN: MAICE order C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), A(MJ2,MJ2), SIG2(0:MJ2), AIC(0:MJ2) C CALL COMAIC( X,N,K,MJ1,SIG2,AIC ) C IMIN = 0 AICM = AIC(0) DO 10 M=1,K IF( AIC(M).LT.AICM ) THEN IMIN = M AICM = AIC(M) END IF CALL RECOEF( X,M,K,MJ1,A(1,M) ) 10 CONTINUE C RETURN E N D SUBROUTINE COMAIC( X,N,K,MJ1,SIG2,AIC ) C C ... This subroutine computes residual variances and AIC's ... C C Inputs: C X(I,J): Householder reduced form C N: Data length C K: Heighest order C MJ1: Adjustable dimension of X C Outputs: C SIG2(I): residual variance of the model with order I C AIC(I): AIC of the model with order I C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ1,1), AIC(0:K), SIG2(0:K) DATA PI2/6.28318531D0/ C PVAR = 0.0D0 DO 10 I=K,0,-1 PVAR = PVAR + X(I+1,K+1)**2 SIG2(I) = PVAR / N 10 AIC(I) = N*DLOG( PI2*SIG2(I) ) + N + 2*(I+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 SUBROUTINE RECOEF( X,M,K,MJ,A ) C C ... Regression coefficients ... C C Inputs: C X(I,J): Householder reduced form C M: Number of actually used regressors C K: Heighest order C MJ: Adjustable dimension of X C Output: C A(I): Vector of regression coefficients C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(MJ,1), A(1) C A(M) = X(M,K+1)/X(M,M) DO 20 I=M-1,1,-1 SUM = X(I,K+1) DO 10 J=I+1,M 10 SUM = SUM - A(J)*X(I,J) 20 A(I) = SUM/X(I,I) C RETURN E N D SUBROUTINE REDUCT( SETX,Z,D,NMK,N0,K,MJ1,X ) C C ... Successive Householder reduction ... C C Inputs: C SETX: Name of the subroutine for making X(I,J) C Z(I): Data vector C D(I): Working area C NMK: Number of actually used observations C N0: Time point of the previous set ofobservations C K: Heighest order of the model C MJ1: Adjustable dimension of X C Output: C X(I,J): data matrix C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION X(MJ1,1) , D(1), Z(1) C L = MIN0( NMK,MJ1 ) K1 = K + 1 N1 = L C CALL SETX( Z,N0,L,K,MJ1,0,X ) CALL HUSHLD( X,D,MJ1,L,K1 ) IF( N1 .GE. NMK ) RETURN C 10 L = MIN0( NMK-N1,MJ1-K1 ) C LK = L + K1 N2 = N0 + N1 CALL SETX( Z,N2,L,K,MJ1,1,X ) CALL HUSHLD( X,D,MJ1,LK,K1 ) N1 = N1 + L IF( N1.LT.NMK ) GO TO 10 C RETURN C E N D SUBROUTINE READTS( IDEV,Y,N ) IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION Y(*) COMMON /CMDATA/ TITLE C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) CLOSE( IDEV ) C RETURN 1 FORMAT( 10A8 ) E N D