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 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 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