C PROGRAM 15.1 SIMSSM C C ... SIMULATION BY GAUSSIAN STATE SPACE MODEL ... C C INPUTS: C M1: TREND ORDER (0,1,2) C M2: SEASONAL ORDER (0 OR 1) C M3: AR ORDER (0,...,20) C N: SIMULATION INTERVAL C NS: Y(NS+1),...,Y(N) ARE ACTUALLY USED C INI: INITIAL VALUE FOR RANDOM NUMBER GENERATOR C PARAMETERS: C MJ: MAXIMUM STATE DIMENSION C NN: MAXIMUM DATA LENGTH C MODIFED 2/16/93 C PARAMETER( MJ=20,NN=1000 ) IMPLICIT REAL*8(A-H,O-Z) INTEGER PERIOD DIMENSION X(MJ), F(MJ,MJ), G(MJ,MJ), H(MJ) DIMENSION Q(MJ,MJ), Y(NN), AR(20), R(1,1) INI = 1992092521 MM = 0 C READ( 5,* ) M1, M2, M3, N, NS, INI READ( 5,* ) SIG2 IF( M1.GT.0 ) THEN READ( 5,* ) TAU1 READ( 5,* ) (X(I),I=1,M1) MM = M1 END IF IF( M2.GT.0 ) THEN READ( 5,* ) TAU2, PERIOD READ( 5,* ) (X(MM+I),I=1,PERIOD-1) MM = MM + PERIOD-1 END IF IF( M3.GT.0 ) THEN READ( 5,* ) TAU3, (AR(I),I=1,M3) READ( 5,* ) (X(MM+I),I=1,M3) END IF INI0 = INI C C ... SET STATE SPACE MODEL ... C CALL SETSEA( M1,M2,M3,PERIOD,AR,TAU1,TAU2,TAU3,SIG2,F,G,H, * Q,R,M,K,MJ ) C C ... SIMULATION ... C CALL SIMSSM( F,G,H,Q,R,X,N,M,1,K,INI,MJ,NN,Y ) C C ... PLOT AND PRINT OUT SIMULATED TIME SERIES ... C CALL PTSIM( Y,N,NS,K,M1,M2,M3,INI0 ) CALL PRSIM( Y,N,NS,M1,M2,M3,INI0 ) C STOP E N D SUBROUTINE SETSEA( M1,M2,M3,IPER,AR,TAU1,TAU2,TAU3,SIG2, * F,G,H,Q,R,M,K,MJ ) C C ... SET STATE SPACE MODEL FOR SEASONAL ADJUSTMENT ... C C INPUTS: C M1,M2,M3: MODEL ORDERS C IPER: NUMBER OF SEASONS IN ONE PERIOD C AR(I): AR COEFFICIENTS C TAU1-TAU3: SYSTEM NOISE VARIANCES C MJ: ADJUSTABLE DIMENSION C OUTPUTS: C F,G,H: MATRICES OF STATE SPACE MODEL C Q,R: SYSTEM NOISE AND OBS NOISE VARIANCES C M: STATE DIMENSION C K: DIMENSION OF SYSTEM NOISE C MODIFIED 2/16/93 C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(MJ,MJ), G(MJ,MJ), H(MJ), Q(MJ,MJ), AR(*), R(1,1) C M = M1 + M2*(IPER-1) + M3 K = ID(M1) + ID(M2) + ID(M3) DO 10 I=1,M H(I) = 0.0D0 DO 10 J=1,M 10 F(I,J) = 0.0D0 DO 20 I=1,M DO 20 J=1,K 20 G(I,J) = 0.0D0 DO 30 I=1,K 30 Q(I,I) = 0.0D0 C IF( M1.GT.0 ) THEN IF( M1.EQ.1 ) F(1,1) = 1.0D0 IF( M1.EQ.2 ) THEN F(1,1) = 2.0D0 F(1,2) =-1.0D0 F(2,1) = 1.0D0 END IF G(1,1) = 1.0D0 H(1) = 1.0D0 Q(1,1) = TAU1 END IF C C ... SEASONAL COMPONENT ... C IF( M2.GT.0 ) THEN L1 = ID(M1) + 1 DO 40 I=1,IPER-1 40 F(M1+1,M1+I) = -1.0D0 DO 50 I=2,IPER-1 50 F(M1+I,M1+I-1) = 1.0D0 G(M1+1,L1) = 1.0D0 H(M1+1) = 1.0D0 Q(L1,L1) = TAU2 END IF C C ... AR COMPONENT ... C IF( M3.GT.0 ) THEN M12 = M1 + M2*(IPER-1) L12 = ID(M1) + ID(M2) + 1 DO 60 I=1,M3 60 F(M12+1,M12+I) = AR(I) DO 70 I=2,M3 70 F(M12+I,M12+I-1) = 1.0D0 G(M12+1,L12) = 1.0D0 H(M12+1) = 1.0D0 Q(L12,L12) = TAU3 END IF C R(1,1) = SIG2 C RETURN E N D SUBROUTINE SIMSSM( F,G,H,Q,R,X,N,M,L,K,IX,MJ,NN,Y ) C C ... SIMULATION BY STATE SPACE MODEL ... C C INPUTS: C F,G,H: COEFFICIENT MATRICES OF STATE SPACE MODEL C Q,R: SYSTEM NOISE AND OBS. NOISE VARIANCES C X: INITIAL STATE VECTOR C M,L,K: DIMENSIONS OF THE STATE, SYSTEM NOISE AND OBS. C IX: INITAL SEED FOR RANDOM NUMBER GENERATOR C MJ: ADJUSTABLE DIMENSION C NN: NUMBER OF DATA C OUTPUT: C Y(I): SIMULATED TIME SERIES C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ), F(MJ,MJ), G(MJ,MJ), H(L,MJ) DIMENSION Q(MJ,MJ), R(L,L), SQ(10,10), SR(10,10) DIMENSION XT(40), T(20), V(20), W(10), Y(NN,L) C C ... CHOLESKY DECOMPOSITION OF Q AND R ... C CALL CHOLES( Q,MJ,K,SQ,10 ) CALL CHOLES( R,L,L,SR,10 ) C DO 100 II= 1,N C C ... X = F*X + V ... C CALL WHITE( K,SQ,IX,V ) C DO 10 I=1,M 10 XT(I) = 0.0D0 DO 20 J=1,M DO 20 I=1,M 20 XT(I) = XT(I) + F(I,J)*X(J) DO 30 J=1,K DO 30 I=1,M 30 XT(I) = XT(I) + G(I,J)*V(J) DO 40 I=1,M 40 X(I) = XT(I) C C ... T = H*X + W ... C CALL WHITE( L,SR,IX,W ) DO 50 I=1,L 50 T(I) = 0.0D0 DO 60 I=1,L DO 60 J=1,M 60 T(I) = T(I) + H(I,J)*X(J) DO 80 I=1,L 80 Y(II,I) = T(I) + W(I) 100 CONTINUE C RETURN E N D SUBROUTINE CHOLES ( X,MJ,N,Y,NJ ) C C ... CHOLESKY DECOMPOSITION ... C C INPUTS: C X(I,J): SYMMETRIC MATRIX C N: DIMENSION OF Z C MJ: ADJUSTABLE DIMENSION OF X C NJ: ADJUSTABLE DIMENSION OF Y C OUTPUT: C Y(I,J): LOWER TRIANGULAR MATRIX; X = Y*Y' C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(MJ,MJ), Y(NJ,NJ) C DO 10 I=1,N DO 10 J=1,N 10 Y(I,J) = 0.0D0 C DO 100 J = 1,N SUM1 = X(J,J) DO 20 K=1,J-1 20 SUM1 = SUM1 - Y(J,K)**2 IF( SUM1.GT.0.0D0 ) Y(J,J) = DSQRT( SUM1 ) IF( SUM1.EQ.0.0D0 ) Y(J,J) = 0.0D0 DO 40 I=J+1,N SUM2 = 0.0D0 DO 30 K=1,J-1 30 SUM2 = SUM2 + Y(I,K)*Y(J,K) 40 Y(I,J) = ( X(I,J) - SUM2 ) / Y(J,J) 100 CONTINUE C RETURN E N D SUBROUTINE WHITE( K,Q,IX,V ) C C ... GENERATE K-DIMENSIONAL GAUSSIAN RANDOM NUMBER ... C C INPUTS: C K: DIMENSION OF THE RANDOM NUMBER C Q: COVARIANCE C IX: INITIAL SEED FOR RANDOM NUMBER GENERATOR C OUTPUT: C V: K-DIMENSIONAL RANDOM NUMBER WITH COVARIANCE Q C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Q(K,K),W(100),V(K) DO 10 I=1,K 10 W(I) = RGAUSS( IX ) DO 30 I=1,K SUM = 0.0D0 DO 20 J=1,K 20 SUM = SUM + Q(J,I)*W(J) 30 V(I) = SUM RETURN E N D DOUBLE PRECISION FUNCTION RUNI( IX ) C C ... UNIFORM RANDOM NUMBER GENERATOR C REAL*8 A, EM DATA EM /2147483647D0/, A/314159269D0/ C IX = MOD( IX*A,EM ) RUNI = IX/EM RETURN E N D DOUBLE PRECISION FUNCTION RUNIF( IX ) C C ... UNIFORM RANDOM NUMBER GENERATOR (CONGRUENTIAL METHOD) ... C REAL*8 A, C, EM DATA EM /1664501D0/, A/1229D0/, C/351750D0/ C IX = MOD( IX*A+C,EM ) RUNIF = IX/EM RETURN E N D DOUBLE PRECISION FUNCTION RGAUSS( IX ) C C ... GAUSSIAN RANDON NUMBER GENERATOR (MARSAGLIA'S ALGORITHM) ... C IMPLICIT REAL*8( A-H,O-Z) DATA ICOUNT /0/ C IF( IX.EQ.0 ) IX = 1990103011 IF( ICOUNT.EQ.0 ) THEN 10 U1 = RUNI(IX) U2 = RUNI(IX) V1 = 2*U1-1 V2 = 2*U2-1 S = V1**2 + V2**2 IF(S.GE.1) GO TO 10 RGAUSS = V1*SQRT( -2*DLOG(S)/S ) ICOUNT = 1 C ELSE RGAUSS = V2*SQRT( -2*DLOG(S)/S) ICOUNT = 0 END IF RETURN E N D SUBROUTINE PTSIM( Y,N,NS,K,M1,M2,M3,INI ) C C ... PLOT SIMULATED TIME SERIES ... C IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72, VNAME*8 DIMENSION Y(N), VNAME(5), VALUE(5) COMMON /CMDATA/ TITLE C VNAME(1) = 'M1 =' VNAME(2) = 'M2 =' VNAME(3) = 'M3 =' VNAME(4) = 'K =' VNAME(5) = 'INI =' VALUE(1) = M1 VALUE(2) = M2 VALUE(3) = M3 VALUE(4) = K VALUE(5) = INI TITLE = ' ' DO 10 I=1,17 10 TITLE = TITLE//' ' CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 15.X: SIMULATION BY S.S.M. ',38,5, * VNAME,VALUE ) WY = 5.0 WX = 16.0 DX = 20.0 FN = N-NS IY = 10 CALL NEWPEN( 2 ) C CALL MAXMIN( Y(NS+1),N-NS,YMIN1,YMAX1,DY1 ) CALL PLOT( 3.0,16.0-SNGL(WY),-3 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY,0.0D0,FN,YMIN1,YMAX1,DX,DY1, * 0.2D0,1,IY,2 ) C CALL NEWPEN( 1 ) CALL PLOTY( Y(NS+1),N-NS,YMIN1,YMAX1,WX,WY,1,1 ) C CALL PLOTE C call plot( 0.0,0.0,999 ) RETURN E N D SUBROUTINE PRSIM( Y,N,NS,M1,M2,M3,INI ) C C ... PRINT OUT SIMULATED TIME SERIES ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N) C NN = N - NS WRITE(6,600) M1, M2, M3, INI, NS WRITE(6,610) NN WRITE(6,620) (Y(I),I=NS+1,N) C RETURN 600 FORMAT( 1H ,'SIMULATED TIME SERIES: M1 =',I2,3X,'M2 =',I2, * 3X,'M3 =',I2,3X,'INI =',I10,3X,'NS =',I4 ) 610 FORMAT( 1H ,I5 ) 620 FORMAT( 1H ,5D13.5 ) E N D FUNCTION ID(M) ID = 0 IF( M.GT.0 ) ID = 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 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