C     PROGRAM 15.2  NGSIM
C
C  ...  SIMULATION BY NON-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        NOISEW:    TYPE OF OBSERVATIONAL NOISE
C        NOISEV:    TYPE OF SYSTEM NOISE
C        WMIN,WMAX: LOWER AND UPPER BOUND OF OBS. NOISE
C        VMIN,VMAX: LOWER AND UPPER BOUND OF SYSTEM NOISE
C        X(I): INITIAL STATE
C     PARAMETERS:
C        MJ:   MAXIMUM STATE DIMENSION
C        NN:   MAXIMUM DATA LENGTH
C     MODIFIED  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)
      DIMENSION  PARAMV(3), PARAMW(3)
      DATA   TAU1/1.0D0/, TAU2/1.0D0/, TAU3/1.0D0/, SIG2/1.0D0/
      INI = 1992092521
C
      READ( 5,* )  M1, M2, M3, N, NS, INI
C
      READ( 5,* )  NOISEW, WMIN, WMAX
      PARAMW(1) = 0.0D0
      IF( NOISEW.NE.2 )  READ(5,*)  PARAMW(2)
      IF( NOISEW.EQ.2 )  READ(5,*)  PARAMW(2), PARAMW(3)
C
      READ( 5,* )  NOISEV, VMIN, VMAX
      PARAMV(1) = 0.0D0
      IF( NOISEV.NE.2 )  READ(5,*)  PARAMV(2)
      IF( NOISEV.EQ.2 )  READ(5,*)  PARAMV(2), PARAMV(3)
C
      MM = 0
      IF( M1.GT.0 )  THEN
         READ( 5,* )  (X(I),I=1,M1)
         MM = M1
      END IF
      IF( M2.GT.0 )  THEN
         READ( 5,* )  PERIOD
         READ( 5,* )  (X(MM+I),I=1,PERIOD-1)
         MM = MM + PERIOD-1
      END IF
      IF( M3.GT.0 )  THEN
         READ( 5,* )  (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  NGSIM( NOISEV,NOISEW,PARAMV,PARAMW,VMIN,VMAX,WMIN,WMAX,
     *             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 NGSIM( NOISEV,NOISEW,PV,PW,X0,X1,Y0,Y1,F,G,H,Q,R,
     *                  X,N,M,L,K,IX,MJ,NN,Y )
C
C  ...  SIMULATION BY NON-GAUSSIAN STATE SPACE MODEL  ...
C
C     INPUTS:
C        NOISEV:  TYPE OF THE SYSTEM NOISE
C        NOISEW:  TYPE OF THE OBSERVATIONAL NOISE
C        PV:      PARAMETER OF THE SYSTEM NOISE DENSITY
C        PW:      PARAMETER OF THE OBSER. NOISE DENSITY
C        X0,X1:   DOMAIN OF SYSTEM NOISE DENSITY
C        Y0,Y1:   DOMAIN OF OBSER. NOISE DENSITY
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)
      DIMENSION FV(0:400), FW(0:400), XV(0:400), XW(0:400)
      DIMENSION PV(3), PW(3)
      EXTERNAL  GAUSS
      EXTERNAL  PEARSN
      EXTERNAL  DBLEXP
      EXTERNAL  USERV
      EXTERNAL  USERW
C
C  ...  DISTRIBUTION FUNCTION OF SYSTEM NOISE AND OBS. NOISE  ...
C
      IF(NOISEV.EQ.0)  CALL  DISTRI( USERV ,PV,X0,X1,FV,XV,DXV )
      IF(NOISEV.EQ.1)  CALL  DISTRI( GAUSS ,PV,X0,X1,FV,XV,DXV )
      IF(NOISEV.EQ.2)  CALL  DISTRI( PEARSN,PV,X0,X1,FV,XV,DXV )
      IF(NOISEV.EQ.3)  CALL  DISTRI( DBLEXP,PV,X0,X1,FV,XV,DXV )
C
      IF(NOISEW.EQ.0)  CALL  DISTRI( USERW ,PW,Y0,Y1,FW,XW,DXW )
      IF(NOISEW.EQ.1)  CALL  DISTRI( GAUSS ,PW,Y0,Y1,FW,XW,DXW )
      IF(NOISEW.EQ.2)  CALL  DISTRI( PEARSN,PW,Y0,Y1,FW,XW,DXW )
      IF(NOISEW.EQ.3)  CALL  DISTRI( DBLEXP,PW,Y0,Y1,FW,XW,DXW )
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  NGNOIS( NOISEV,FV,XV,DXV,IX,SQ,K,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  NGNOIS( NOISEW,FW,XW,DXW,IX,SR,L,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  NGNOIS( NOISE,F,X,DX,IX,Q,K,V )
C
C  ...  K-DIMENSIONAL NON-GAUSSIAN RANDOM NUMBER  ...
C
C     INPUTS:
C        NOISE:NOISE TYPE
C        F:    DISTRIBUTION FUNCTION OF NON-GAUSSIAN DISTRIBUTION
C        IX:   INITIAL SEED FOR RANDOM NUMBER GENERATOR
C        Q:    COVARIANCE
C        K:    DIMENSION OF THE RANDOM NUMBER
C     OUTPUT:
C        V:    K-DIMENSIONAL RANDOM NUMBER WITH COVARIANCE Q
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Q(10,10), W(10), V(K), F(0:400), X(0:400)
C
      DO 10 I=1,K
      IF( NOISE.LT.0 )  W(I) = RNG2( IX,NOISE )
   10 IF( NOISE.GE.0 )  W(I) = RNG( IX,F,X,DX )
      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
      SUBROUTINE  DISTRI( FUNCT,PARAM,XMIN,XMAX,F,X,DX )
C
C  ...  DISTRIBUTION FUNCTION  ...
C
C     INPUTS:
C        FUNCT:  FUNCTION NAME
C        PARAM:  PARAMETER OF THE DENSITY
C        XMIN,XMAX:  DOMAIN OF THE  DENSITY
C     OUTPUTS:
C        F:      DISTRIBUTION FUNCTION
C        X:      LOCATION OF NODE
C        DX:     SIZE OF EACH BIN
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  F(0:400), P(0:400), X(0:400), PARAM(3)
      K = 400
      DX = (XMAX-XMIN)/K
C
      DO 10 I=0,K
      X(I) = XMIN + I*DX
   10 P(I) = FUNCT( X(I),PARAM )
      DO 20 I=0,K
   20 F(I) = 0.0D0
      DO 30 I=1,K
   30 F(I) = F(I-1) + ( P(I-1) + P(I) )*DX/2.0D0
      DO 40 I=1,K
   40 F(I) = F(I)/F(K)
C
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  RNG( IX,F,X,DX )
C
C  ...  NON-GAUSSIAN RANDON NUMBER GENERATOR  ...
C
C     INPUTS:
C        IX:     SEED
C        F:      DISTRIBUTION FUNCTION
C        X:      LOCATION
C        DX:     SIZE OF EACH BIN
      IMPLICIT REAL*8( A-H,O-Z)
      DIMENSION  F(0:400), X(0:400)
C
      IF( IX.EQ.0 )  IX = 1990103011
      U = RUNI(IX)
C
      I = 0
   10 I = I+1
      IF( U.GT.F(I) )  GO TO 10
      IF( U.EQ.F(I) )  THEN
        V = X(I)
      ELSE
        V =((U-F(I-1))/( F(I)-F(I-1)) )*DX+X(I-1)
      END IF
      RNG = V
C
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  RNG2( IX,NOISE )
C
C  ...  NON-GAUSSIAN RANDON NUMBER GENERATOR  ...
C
C     NOISE = -1:  CAUCHY RANDOM NUMBER
C           = -2:  EXPONENTIAL DISTRIBUTION
C           = -3:  DOUBLE EXPONENTIAL DISTRIBUTION
C
      IMPLICIT REAL*8( A-H,O-Z)
      DATA  PI /3.1415926535D0/
C
      IF( IX.EQ.0 )  IX = 1990103011
   10 U = RUNI(IX)
C
C  ...  CAUCHY DISTRIBUTION  ...
C
      IF( NOISE.EQ.-1 )  THEN
         IF( U.EQ.0.5D0 )  GO TO 10
         RNG2 = DTAN( PI*U )
      END IF
C
C  ...  EXPONENTIAL DISTRIBUTION  ...
C
      IF( NOISE.EQ.-2 )  THEN
         IF( U.EQ.0.0D0 )  GO TO 10
         RNG2 = -DLOG( U )
      END IF
C
C  ...  DOUBLE EXPONENTIAL DISTRIBUTION  ...
C
      IF( NOISE.EQ.-3 )  THEN
         IF( U.EQ.0.0D0 )  GO TO 10
         RNG2 = DEXP( -DEXP( U ) )
      END IF
C
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  USERV( X,PARAM )
C
C  ...  User supplied density function  ...
C       (The double exp. dist. with mean 0 )
C     Inputs:
C        X:
C        PARAM(1):  lambda
C     Output:
C        USERV:     density at X
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PARAM(2)
      DATA  C/0.5771D0/
C
      USERV = DEXP( X-C - DEXP(X-C) )
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  USERW( X,PARAM )
C
C  ...  User supplied density function  ...
C       (The double exp. dist. with mean 0 )
C     Inputs:
C        X:
C        PARAM(1):  lambda
C     Output:
C        USERW:     density at X
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PARAM(2)
      DATA  C/0.5771D0/
C
      USERW = DEXP( X-C - DEXP(X-C) )
      RETURN
      E N D
      FUNCTION ID( M )
      ID = 0
      IF( M.GT.0 )  ID = 1
      RETURN
      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  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  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
      DOUBLE PRECISION FUNCTION  GAUSS( X,PARAM )
C
C  ...  Gaussian (normal) distribution  ...
C
C     Inputs:
C        X:
C        PARAM(1):  mean
C        PARAM(2):  variance
C     Output:
C        GAUSS:     density at X
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PARAM(2)
      DATA  C1  /2.506628275D0/
C
      GAUSS = DEXP( -(X-PARAM(1))**2/(2*PARAM(2)) )/(C1*DSQRT(PARAM(2)))
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  PEARSN( X,PARAM )
C
C  ...  Pearson family of  distributions  ...
C
C     Inputs:
C        X:
C        PARAM(1):  location parameter, mu
C        PARAM(2):  dispersion parameter, tau square
C        PARAM(3):  shape parameter
C     Output:
C        PEARSN:    density at X
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  PARAM(3)
      DATA  PI/3.1415926535D0/
C
      PEARSN = DGAMMA(PARAM(3))/DGAMMA(PARAM(3)-0.5D0)
     *                  /DSQRT(PI)*PARAM(2)**(PARAM(3)-0.5D0)
     *                  /((X-PARAM(1))**2 + PARAM(2))**PARAM(3)
      RETURN
C
      END
      DOUBLE PRECISION  FUNCTION  DGAMMA( X )
C
C  ...  Gamma function  ...
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  A(0:10)
      DATA  A /0.999999999999269D0, 0.42278433696202D0,
     *  0.41184025179616D0, 0.08157821878492D0, 0.0742379070629D0,
     * -0.00021090746731D0, 0.01097369584174D0,-0.00246674798054D0,
     *  0.00153976810472D0,-0.00034423420456D0, 0.00006771057117D0/
C
      DGAM = 1.0D0
      Y = X
      IF( X.GT.3.0D0 )  THEN
   10   Y = Y-1
        DGAM = DGAM*Y
        IF( Y.GT.3.0D0 )  GO TO 10
      END IF
      IF( X.LE.2.0D0 )  THEN
   20   DGAM = DGAM/Y
        Y = Y+1
        IF( Y.LE.2.0D0 )  GO TO 20
      END IF
C
      Z   = 1.0D0
      SUM = 0.0D0
      DO 30 I=0,10
      SUM = SUM + A(I)*Z
   30 Z = Z*(Y-2)
      DGAMMA = DGAM*SUM
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  DBLEXP( X,PARAM )
C
C  ...  double exponential distribution  f(x) = exp(x - exp(x))  ...
C
C     Inputs:
C        X:
C     Output:
C        DBLEXP:    density at X
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      DBLEXP = DEXP( X-DEXP(X) )
      RETURN
C
      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