C PROGRAM 8.1 LSAR1 C C ... Decomposition of time interval to stationary subintevals ... C C Inputs: C LAG: Highest order od AR model C NS: Basic local span C The following inputs are required in the subroutine READTS. C TITLE: Caption of the data set C N: Data length C Y(I): Time series, (i=1,...,N) C Parameters: C IDEV: Input device for time series C NMAX: Adjustable dimension of Y (NMAX.GE.N) C KMAX,MJ2: Adjustable dimensions C NF: Number of frequencies for computing spectrum C PARAMETER( NMAX=300,KMAX=20,MJ2=KMAX+1,MJ1=100,NF=100,IDEV=1) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), X(MJ1,MJ2), D(MJ1), U(MJ2,MJ2) DIMENSION SP(0:NF), AS(KMAX), A(KMAX), B(1) EXTERNAL SETXAR C C ISW = 0 READ( 5,* ) LAG, NS C C ... Read time series ... C CALL READTS( IDEV,Y,N ) C IF = 0 CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) C NBLOCK = N/NS DO 100 II=1,NBLOCK C L = NS*(II-1) WRITE(6,600) L IF( II.EQ.NBLOCK ) NS = N - NS*(II-1) - LAG LK = L + LAG C C ... Locally stationary time series ... C CALL LOCAL( SETXAR,Y,X,U,D,LAG,L,NS,LAG,IF,MJ1,MJ2,A,MF,SIG2) C IF( IF .EQ. 2 ) LK0 = LK + 1 IF( IF.EQ.2 .AND. II.GT.1 ) THEN CALL ARMASP( AS,MFS,B,0,SIG2S,NF,SP ) CALL PTLSP( Y,N,SP,NF,LK0S,LK2S ) END IF MFS = MF SIG2S = SIG2 LK0S = LK0 LK2S = LK + NS DO 20 I=1,MF 20 AS(I) = A(I) C 100 CONTINUE CALL ARMASP( AS,MFS,B,0,SIG2S,NF,SP ) CALL PTLSP( Y,N,SP,NF,LK0S,LK2S ) CALL PLOTE C call plot( 0.0,0.0,999 ) C STOP 600 FORMAT( 1H ,'L =',I5 ) E N D SUBROUTINE LOCAL( SETX,Z,X,U,D,LAG,N0,NS,K,IF,MJ1,MJ2, * A,MF,SDF ) C C ... Locally stationary AR model ... C C Inputs: C SETX: Name of the subroutine for making X(I,J) C Z(I): Data vector C D,U: Working area C LAG: Highest order of the model C N0: Time point of the previous set ofobservations C NS: Number of new observations C MJ1,MJ2: Adjustable dimension C Output: C A(I): AR coefficients of the current model C MF: Order of the current model C SDF: Innovation variance of the current model C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Z(1) DIMENSION X(MJ1,1), U(MJ2,1), D(1), A(1) DIMENSION B(20), AA(20,20), AIC(0:20), SIG2(0:20) EXTERNAL SETX C K1 = K + 1 K2 = K1*2 NN0 = N0 + LAG + 1 NN1 = N0 + LAG + NS WRITE(6,600) NN0, NN1 C CALL REDUCT( SETX,Z,D,NS,N0,K,MJ1,X ) CALL REGRES( X,K,NS,MJ1,20,AA,SIG2,AIC,MS ) C SDS = SIG2(MS) DO 10 I=1,MS 10 B(I) = AA(I,MS) IF( IF.EQ.0 ) THEN CALL COPY( X,K1,0,0,MJ1,MJ2,U ) AICS = AIC(MS) AIC0 = AIC(MS) WRITE(6,610) NS, MS, SDS, AICS ELSE C AICS = AIC(MS) + AICF AIC0 = AIC(MS) WRITE(6,620) NF, NS, MS, SDS, AICS CALL COPY( X,K1,0,K2,MJ1,MJ1,X ) CALL COPY( U,K1,0,K1,MJ2,MJ1,X ) CALL HUSHLD( X,D,MJ1,K2,K1 ) NP = NF + NS CALL REGRES( X,K,NP,MJ1,20,AA,SIG2,AIC,MP ) C AICP = AIC(MP) SDP = SIG2(MP) DO 20 I=1,MP 20 A(I) = AA(I,MP) WRITE(6,630) NP, MP, SDP, AICP IF( AICS.GE.AICP ) GO TO 40 WRITE(6,640) CALL COPY( X,K1,K2,0,MJ1,MJ2,U ) END IF C IF = 2 NF = NS MF = MS AICF = AIC0 DO 30 I=1,MF 30 A(I) = B(I) SDF = SDS C GO TO 50 40 IF = 1 CALL COPY( X,K1,0,0,MJ1,MJ2,U ) WRITE(6,650) SDF = SDP MF = MP AICF = AICP NF = NF + NS 50 CONTINUE C RETURN 600 FORMAT( 1H0,'<<< NEW DATA (N =',I4,' ---',I4,') >>>' ) 610 FORMAT( 1H ,' INITIAL MODEL: NS =',I4,/,10X,'MS =',I2,3X, 1 'SDS =',D13.6,3X,'AICS =',F12.3 ) 620 FORMAT( 1H ,' SWITCHED MODEL: (NF =',I4,', NS =',I4,1H), 2 /,10X,'MS =',I2,3X,'SDS =',D13.6,3X,'AICS =',F12.3 ) 630 FORMAT( 1H ,' POOLED MODEL: (NP =',I4,1H), 3 /,10X,'MP =',I2,3X,'SDP =',D13.6,3X,'AICP =',F12.3 ) 640 FORMAT( 1H ,30X,'*** SWITCHED MODEL ACCEPTED ***' ) 650 FORMAT( 1H ,30X,'*** POOLED MODEL ACCEPTED ***' ) E N D SUBROUTINE PTLSP( Y,N,SP,M,N1,N2 ) C C ... Plot original data and estimated local spectrum ... C C Inputs: C Y(I): Original time series C N: Data length C SP: Power spectrum of the current model C N1: Initial point of the current block C N2: End point of the current block C IMPLICIT REAL*8(A-H,O-Z) CHARACTER VNAME*8 DIMENSION Y(N), SP(0:M) DIMENSION VNAME(5), YM(5) DATA ICOUNT/0/ C WX = 24.0 WY = 6.0 WX1 = 3.0 ICOUNT = ICOUNT + 1 IF( ICOUNT.EQ.1 ) THEN CALL HEADER( 'LOCALLY STATIONARY AR MODEL',27,0,VNAME,YM ) CALL DRAWY( 'ORIGINAL DATA',13,2.5D0,10.5D0,Y,N,WX,WY,1,0) CALL PLOT( 20.0,-1.0,-3 ) END IF IF( ICOUNT.EQ.13 ) THEN C CALL PLOTI C call plot( 0.0,0.0,777 ) CALL PLOT( 2.5,17.0,-3 ) END IF IF( MOD(ICOUNT,6).EQ.1 ) CALL PLOT( -20.0,-4.5,-3 ) IF( MOD(ICOUNT,6).NE.1 ) CALL PLOT( 4.0, 0.0,-3 ) CALL NEWPEN( 1 ) CALL NUMBER( 0.0,3.2,0.25,SNGL(N1),0.0,-1 ) CALL PLOT( 1.2,3.3,3 ) CALL PLOT( 1.5,3.3,2 ) CALL NUMBER( 1.7,3.2,0.25,SNGL(N2),0.0,-1 ) CALL MAXMIN( SP(0),M,SPMIN,SPMAX,DSP ) CALL AXISXY( 0.0D0,0.0D0,WX1,WX1,0.0D0,0.5D0,SPMIN,SPMAX,0.5D0, * DSP,0.25D0,1,5,1 ) CALL NEWPEN( 1 ) CALL PLOTY ( SP(0),M+1,SPMIN,SPMAX,WX1,WX1,IPOS,1 ) RETURN E N D SUBROUTINE AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY) C C ... This subroutine draws X and Y axes. C C Inputs: C X,Y: location of the left bottom of the figure C WX,WY: width and height of the figure C X0,X1: lower and upper bounds of the X axis C Y0,Y1: lower and upper bounds of the Y axis C DX,DY: step width in X and Y axes C DWC: size of characters C IBOX: = 0 draw X and Y axes only C = 1 draw window doundary C IX,IY: number of additional click in each step C Modified: 8/31/90, 11/21/90 C REAL*8 X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC C FX = X FY = Y FWX = WX FWY = WY FWC = DWC CALL NEWPEN( 2 ) CALL PLOT( FX,FY,-3 ) CALL PLOT( 0.0,FWY,3 ) CALL PLOT( 0.0,0.0,2 ) CALL PLOT( FWX,0.0,2 ) IF ( IBOX.EQ.1 ) THEN CALL PLOT( FWX,FWY,2 ) CALL PLOT( 0.0,FWY,2 ) END IF C C ... draw X axis ... C WC = DWC NX = (X1-X0)/DX + 0.001 NY = (Y1-Y0)/DY + 0.001 MX = -DLOG10( DX ) + 0.001 IF(DX.LT.1.0D0) MX = -DLOG10( DX ) + 0.999 IF( MX.LE.0 ) MX = -1 D = WX*DX/(X1-X0) DO 10 I=1,NX+1 XN = X0 + DX*(I-1) XX = D*(I-1) CALL PLOT( XX,-FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF LX = 0 IF(XN.GT.0.0D0) LX = ALOG10( XN ) + 0.001 IF( MX.GT.0 ) XXX = WC*(LX+MX)*0.5 + WC*4/7.0 IF( MX.LT.0 ) XXX = WC*LX*0.5 + WC*4/7.0 CALL NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX ) 10 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 20 I=1,NX+1 DO 20 J=1,IX XX = D*(I-1) + (D*J)/IX IF( XX.LE.WX ) THEN CALL PLOT( XX,FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF END IF 20 CONTINUE CALL NEWPEN( 2 ) END IF C C ... draw Y axis ... C D = WY*DY/(Y1-Y0) MY = -DLOG10( DY ) + 0.99999 IF( MY.LE.0 ) MY = -1 DO 30 I=1,NY+1 YN = Y0 + DY*(I-1) YY = D*(I-1) CALL PLOT( -FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF LY = 1 IF(YN.LT. 0.0) LY = 2 IF(YN.GE. 1.0) LY = ALOG10( YN )+1.00001 IF(YN.LE.-1.0) LY = ALOG10(-YN )+2.00001 IF( MY.GT.0 ) YYY = WC*(LY+MY+1) + 0.20 IF( MY.LT.0 ) YYY = WC*LY + 0.20 CALL NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY ) 30 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 40 I=1,NY+1 DO 40 J=1,IY YY = D*(I-1) + (D*J)/IY IF( YY.LE.WY ) THEN CALL PLOT( FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF END IF 40 CONTINUE CALL NEWPEN( 2 ) END IF C RETURN E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... This subroutine draw title of data, program name, names C and values of parameters ... C C Inputs: C PNAME: program name C NP: number of characters of program name C M: number of parameters C VNAME(I): name of the I-th parameter C VALUE(I): value of the I-th parameter C Input via common: C TITLE: title of the data set C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 3.5 ELSE X1 = X1 + 5.0 X2 = X2 + 5.0 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE DELX( XMIN,XMAX,DX ) C C ... This subroutine determines step width in drawing axis ... C C Inputs: C XMIN: the minimum value in the axis C XMAX: the maximum value in the axis C Output: C DX: step size C IMPLICIT REAL*8(A-H,O-Z) DDX = XMAX - XMIN DX = INT( DLOG10( DDX ) ) IF( DDX.LT.1.0D0 ) DX = DX - 1 DX = 10.0D0**DX IF( DDX/DX.LT.3.0D0 ) DX = DX/2.0 IF( DDX/DX.LT.3.0D0 ) DX = DX/2.0 RETURN E N D SUBROUTINE DRAWY( FNAME,NC,XO,YO,Y,N,WX,WY,IPOS,ISW ) C C ... This subroutine draws time series data ... C C Inputs: C FNAME: title of the data set C NC: number of charcters of the name C XO,YO: origin of the figure C Y(I): time series C N: data length C WX,WY: width and height of the figure C IPOS: starting position of time series C ISW: = 1 connect data by line C = 2 connect data by dashed line C IMPLICIT REAL*8(A-H,O-Z) CHARACTER FNAME*80 DIMENSION Y(N) C CALL DELX( 0.0D0,DFLOAT(N),DX ) CALL MAXMIN( Y,N,YMIN,YMAX,DY ) CALL AXISXY( XO,YO,WX,WY,0.0D0,DBLE(N-1+IPOS),YMIN,YMAX, * DX,DY,0.25D0,1,10,1 ) IF( ISW.LE.1 ) CALL PLOTY (Y,N,YMIN,YMAX,WX,WY,IPOS,1 ) IF( ISW.EQ.1 ) CALL PLOTY2(Y,N,YMIN,YMAX,WX,WY,IPOS,ISW,0.D0) IF( ISW.EQ.2 ) CALL PLOTY2(Y,N,YMIN,YMAX,WX,WY,IPOS,ISW,YMIN) CALL NEWPEN( 1 ) CALL SYMBOL( 0.0,SNGL(WY)+0.2,0.3,FNAME,0.0,NC ) RETURN END SUBROUTINE MAXMIN( X,N,XMIN0,XMAX0,DXL ) C C ... This subroutine determines the minimum, the maximum and C the step width ... C C Inputs: C X(I): data C N: data length C Outputs: C XMIN0: the minimum bound for the figure C XMAX0: the maximum bound for the figure C DXL: step width C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C XMIN = 1.0D30 XMAX =-1.0D30 C DO 10 I=1,N IF( X(I) .LT. XMIN ) XMIN = X(I) 10 IF( X(I) .GT. XMAX ) XMAX = X(I) DX = XMAX-XMIN IF( DLOG10(DX) .GE. 0.0D0 ) DXL = INT( DLOG10(DX) ) IF( DLOG10(DX) .LT. 0.0D0 ) DXL = INT( DLOG10(DX) )-1.0 DXL = 10.0**DXL IF( DX/DXL.GT.6.0D0 ) DXL = DXL*2.0 DIF = INT( DX/DXL ) XMIN0 = INT( XMIN/DXL )*DXL XMAX0 = XMIN0 + DIF*DXL IF( XMIN0 .GT. XMIN ) XMIN0 = XMIN0 - DXL 30 IF( XMAX0 .GE. XMAX ) GO TO 40 XMAX0 = XMAX0 + DXL GO TO 30 40 CONTINUE C RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... This subroutine draws time vs data plot of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 connect data by straight line C = 2 connect data by dashed line C = 3 express data by character (character code=ISW) C Modified: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN E N D SUBROUTINE PLOTY2( Y,N,YMIN,YMAX,WX,WY,IOFF,ISW,YLEVEL ) C C ... This subroutine draws box car graph of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 draw box with thick line C = 2 draw box with thin line C YLEVEL: ground level for the box car C REAL*8 Y(N), YMIN, YMAX, WX, WY, YLEVEL C DX = WX/(N-1+IOFF) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 3 ) IF( N.GT.75 ) CALL NEWPEN( 2 ) IF( N.GT.150) CALL NEWPEN( 1 ) IF( ISW.EQ.2 ) CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IOFF) YY = (Y(I) - YMIN)*DY Y0 = (YLEVEL-YMIN)*DY CALL PLOT( XX,Y0,3 ) CALL PLOT( XX,YY,2 ) 100 CONTINUE CALL NEWPEN( 2 ) C RETURN E N D SUBROUTINE READTS( IDEV,Y,N ) REAL*8 Y(N) CHARACTER*72 TITLE COMMON /CMDATA/ TITLE C C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) C CLOSE( IDEV ) RETURN 1 FORMAT( A72 ) 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 SETXAR( Z,N0,L,K,MJ1,JSW,X ) C C ... Data matrix for AR model ... C C Inputs: C Z(I): Time series C N0: Origin of the current observations C L: Number of current observations C K: Number of regressors 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 JJ = N0 + K + I X(II,K+1) = Z(JJ) DO 10 J=1,K JJ = JJ - 1 10 X(II,J) = Z(JJ) C RETURN E N D SUBROUTINE COPY( X,K,II,JJ,MJ1,MJ2,Y ) C C ... Make a copy of X on Y C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION X(MJ1,1) , Y(MJ2,1) C DO 10 I=1,K I1 = I + II I2 = I + JJ DO 10 J=1,K 10 Y(I2,J) = X(I1,J) C RETURN E N D SUBROUTINE ARMASP( A,M,B,L,SIG2,NF,SP ) C C ... Logarithm of the power spectrum of the ARMA model ... C C Inputs: C M: AR order C L: MA order C A(I): AR coefficient C B(I): MA coefficient C SIG2: Innovation variance C NF: Number of frequencies C Output: C SP(I): Power spectrum (in log scale) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M), B(L) DIMENSION SP(0:NF), H(0:500), FR(0:500), FI(0:500) C H(0) = 1.0D0 DO 10 I=1,M 10 H(I) = -A(I) C CALL FOURIE( H,M+1,NF+1,FR,FI ) C DO 20 I=0,NF 20 SP(I) = SIG2/( FR(I)**2 + FI(I)**2 ) C H(0) = 1.0D0 DO 30 I=1,L 30 H(I) = -B(I) CALL FOURIE( H,L+1,NF+1,FR,FI ) DO 40 I=0,NF 40 SP(I) = SP(I)*( FR(I)**2 + FI(I)**2 ) C DO 50 I=0,NF 50 SP(I) = DLOG10( SP(I) ) C RETURN E N D SUBROUTINE FOURIE( X,N,M,FC,FS ) C C ... Discrete Fourier transformation by Goertzel method ... C C Inputs: C X(I): data (I=1,N) C N: data length C M: number of Fourier components C FC(J): Fourier cosine transform (J=1,M) C FS(J): Fourier sine transform (J=1,M) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N), FC(M), FS(M) DATA PI/3.14159265358979D0/ C W = PI/(M-1) DO 20 I=1,M CI = DCOS(W*(I-1)) SI = DSIN(W*(I-1)) T1 = 0.0 T2 = 0.0 DO 10 J=N,2,-1 T0 = 2*CI*T1 - T2 + X(J) T2 = T1 10 T1 = T0 FC(I) = CI*T1 - T2 + X(1) 20 FS(I) = SI*T1 C RETURN E N D