C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE ZREGR(IERROR,N,M,T,TMIN,TMAX,F,WK,A,AINV,S,C,RESID,COV) C The modelled regession time series T(I) (I=1,2,...,N) is given by C (see vector notation on the RIGHT): C M C T(I) = SUM C(J)*F(I,J) <== T=FC (1) C J=1 C where C(J)=regression coefficient for J-th input surrogate time C series F(I,J) C C Method: Least squares analysis applied to (1) yields the following C linear matrix equation: C AC = S <==> (F'F)C=F'T (2) C C where C = model coefficient vector to be solved and S = data source C vector. It follows that C N C S(J) = SUM T(I)*F(I,J) <== F'T (3) C I=1 C and N C A(J,K) = SUM F(I,J)*F(I,K) <== F'F (4) C I=1 C NOTE(s): The solution is (in vector notation): C -1 -1 C C = (F'F) F'T = A S (5) C C The covariance matrix of the C coefficients is: C 2 -1 C COV = sigma (F'F) (6) C where C 2 1 N 2 C sigma = --- SUM Residual(I) (7) C N-M I=1 C with Residual(I) being the original time series T(I) MINUS the time C series T(I) modelled by (1) above. Under the null hypothesis C (all cross correlations between the coefficients are zero) that C Var(T|F)=sigma**2 I, where I is the N X N unit matrix. C C Input Data: C IERROR = Error flag. If set to 1, standard error analysis is C performed for the coefficient vector C. Setting IERROR C to any other integer SKIPS this C N = Time series length C M = Number of surrogate series F applied in (1) C (Note: This routine effectively solves N equations for M C constants in which M must be less than N) C T(N) = measured time series being modelled C TMIN,TMAX = Values of input series T that are less than TMIN C or greater than TMAX will NOT be included in the C regression (these are the flagging values) C F(N,M) = J-th (J=1,2,..,M) surrogate series, each of length N C Output Data: C C(M) = solved coefficient vector C RESID(N) = Residual error time series (Original series T(I) C MINUS the modelled T(I) from (1) above) C COV(M,M) = Covariance matrix for coefficients C (SKIPPED if C IERROR is not equal to 1) C EXAMPLE: C PARAMETER(IERROR=1,TMIN=-1.1,TMAX=1.1) C PARAMETER(N=20,M=7) C REAL T(N),F(N,M),WK(N,2*N),A(M,M),AINV(M,M),S(M),C(M), C & RESID(N),COV(M,M) C C DO I=1,N C T(I)=COS(2.*I) C F(I,1)=1.0 C F(I,2)=COS(2.05*I) C F(I,3)=SIN(0.4*I) C F(I,4)=COS(0.6*I) C F(I,5)=SIN(0.8*I) C F(I,6)=COS(0.3*I) C F(I,7)=SIN(0.9*I) C ENDDO C CALL ZREGR(IERROR,N,M,T,TMIN,TMAX,F,WK,A,AINV,S,C,RESID,COV) C WRITE(*,*)'I,T(I),T_model(I)' C DO I=1,N C WRITE(*,*) I,T(I),T(I)-RESID(I) C ENDDO C C STOP C END C C Address questions/comments to: C Dr. Jerry Ziemke C NASA/Goddard Space Flight Center PH: (301) 614-6034 C Code 614 FAX: (301) 614-5903 C Greenbelt, MD 20771 C (Affil: GESTAR, Morgan St. Univ., Baltimore, MD) C e-mail: jerald.r.ziemke@nasa.gov REAL T(N),F(N,M),WK(N,2*N),A(M,M),AINV(M,M),S(M),C(M), & RESID(N),COV(M,M) NGOOD=0 DO I=1,N IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) NGOOD=NGOOD+1 ENDDO DO J=1,M S(J)=0. DO I=1,N IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) THEN S(J)=S(J)+T(I)*F(I,J) ! <== S=F'T ENDIF ENDDO ENDDO !Construct matrix A: DO K=1,M DO J=1,M A(J,K)=0. DO I=1,N IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) THEN A(J,K)=A(J,K)+F(I,J)*F(I,K) ! <== A=F'F ENDIF ENDDO ENDDO ENDDO !Solve for coef vector C: IPRINT=0 ! ==> If 1, print out solution each CALL ICHECK=0 ! ==> If 1, print out error analysis each CALL DO I=1,M DO J=1,M WK(I,J)=A(I,J) ENDDO ENDDO DO I=1,M DO J=M+1,2*M IF (M+I.EQ.J) THEN WK(I,J)=1.0 ELSE WK(I,J)=0.0 ENDIF ENDDO ENDDO DO K=1,M DO I=1,M IF (I.NE.K) THEN IF (WK(K,K).EQ.0.0) THEN L=1 30 IF (WK(L,K).EQ.0.0) THEN L=L+1 GOTO 30 ENDIF DO J=K,2*M WK(K,J)=WK(K,J)+WK(L,J) ENDDO ENDIF U=-WK(I,K)/WK(K,K) DO J=K+1,2*M WK(I,J)=WK(I,J)+U*WK(K,J) ENDDO ENDIF ENDDO ENDDO DO J=1,M DO I=1,M AINV(I,J)=WK(I,J+M)/WK(I,I) ENDDO ENDDO IF (IPRINT.EQ.1) THEN DO J=1,M WRITE(*,'(1X,A7,I3,A19)')'Column ',J,' of inverse matrix:' DO I=1,M WRITE(*,*) A(I,J) ENDDO ENDDO ENDIF IF (ICHECK.EQ.1) THEN WRITE(*,*)'-----------------------------------------------' WRITE(*,*)'Simple error check by multiplying A by AINV:' WRITE(*,*)'(this should ideally yield the identity matrix)' DO J=1,M WRITE(*,'(1X,A7,I3,A27)') & 'Column ',J,' of product matrix A*AINV:' DO I=1,M SUM=0.0 DO K=1,M SUM=SUM+A(I,K)*AINV(K,J) ENDDO WRITE(*,*) SUM ENDDO ENDDO ENDIF !Coefficient vector solution: DO I=1,M C(I)=0. DO J=1,M C(I)=C(I)+AINV(I,J)*S(J) ENDDO ENDDO IF (IERROR.EQ.1) THEN DO I=1,N IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) THEN SUM=0 DO J=1,M SUM=SUM+C(J)*F(I,J) ENDDO RESID(I)=T(I)-SUM ENDIF ENDDO SUM=0 DO I=1,N IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) SUM=SUM+RESID(I)**2 ENDDO SIGSQR=SUM/(NGOOD-M) DO K=1,M DO J=1,M COV(J,K)=SIGSQR*AINV(J,K) ENDDO ENDDO ENDIF RETURN END