C     PROGRAM 14.1  NGSMTH
C
C  ...  NON-GAUSSIAN SMOOTHING  ...
C
C     INPUTS:
C        NOISEV:   TYPE OF SYSTEM NOISE DENSITY (0,1,2,3)
C        NOISEW:   TYPE OF OBSER. NOISE DENSITY (0,1,2,3,4)
C        TAU2:     VARIANCE OR DISPERSION OF SYSTEM NOISE
C        SIG2:     VARIANCE OR DISPERSION OF OBSERVATION NOISE
C        BV:       SHAPE PARAMETER OF SYSTEM NOISE (FOR NOISEV = 2)
C        BW:       SHAPE PARAMETER OF OBSER. NOISE (FOR NOISEW = 2)
C     PARAMETERS:
C        MJ:       ADJUSTABLE DIMENSION OF Y, PS, FS, SS (MJ >= N)
C        K:        NUMBER OF INTERVALS + 1
C     @NFILTER.F:   5/14/85, 6/11/87, 5/19/88, 2/14/91
C     MODIFIED  2/16/93
C
      PARAMETER( MJ=500, K=201 )
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER  TITLE*72
      DIMENSION  Y(MJ), P(K), F(K), S(K), T(K), Q(-K:K)
      DIMENSION  YM(10), ST(7), TREND(MJ,7), LOC(MJ)
      REAL*4    PS(K,MJ), FS(K,MJ), SS(K,MJ)
      COMMON   /COMAAA/  Y
      COMMON   /C91214/  XMIN, XMAX, FIGMIN, FIGMAX, YM
      COMMON   /C91215/  NOISEV, NOISEW, INITD, ITF, ITH
      COMMON   /C91216/  B, OUTMIN, OUTMAX
      COMMON   /C91218/  SIG2, TAU2, BW, BV
      COMMON   /C91219/  N, KK, NS, NFE, NPE
      COMMON   /CMDATA/  TITLE
      COMMON    / DDD /  FF , AIC , SD
C     NAMELIST  /PARAM/  XMIN, XMAX, TAU20, SIG20, BC, IOPT, B,
C    *                   OUTMIN, OUTMAX, NS, NFE, NPE, FIGMIN, FIGMAX,
C    *                   NOISEV, NOISEW, ITF, ITH
C
      EXTERNAL  GAUSS
      EXTERNAL  PEARSN
      EXTERNAL  TWOEXP
      EXTERNAL  USERV
C
C  ...  READ DATA  ...
C
      CALL  READTS( 1,Y,N )
      CALL  MOMENT( Y,N,YM(1),YM(2) )
C
      LOC(1) = 0
      KK = K
      CALL  DEFALT
      READ( 5,* )  NOISEV, TAU2, BV
      READ( 5,* )  NOISEW, SIG2, BW
C
      DX = (XMAX-XMIN)/(K-1)
      CALL  IDIST( F,K,YM(1),YM(2),XMIN,DX )
      CALL  NORMLZ( F,K,DX,SSUM )
      IF( NOISEV.EQ.0 )  CALL  TRANS( USERV ,K,DX,TAU2,BV,Q )
      IF( NOISEV.EQ.1 )  CALL  TRANS( GAUSS ,K,DX,TAU2,BV,Q )
      IF( NOISEV.EQ.2 )  CALL  TRANS( PEARSN,K,DX,TAU2,BV,Q )
      IF( NOISEV.EQ.3 )  CALL  TRANS( TWOEXP,K,DX,TAU2,BV,Q )
      CALL  NGSMTH( Y,P,F,S,T,N,K,DX,XMIN,Q,FF,PS,FS,SS,LOC )
C
      DO 30 I=1,NPE
      DO 10 J=1,K
   10 S(J) = SS(J,I)
      CALL  PINTVL( S,K,XMIN,DX,ST )
      DO 20 J=1,7
   20 TREND(I,J) = ST(J) + DX*LOC(I)
   30 CONTINUE
C
      CALL  PRNGSM( NOISEV,NOISEW,TAU2,SIG2,FF,TREND,MJ,N )
C
C      CALL  PTNGSM( N,NPE,NOISEV,NOISEW,TAU2,SIG2,FF,B,Y,TREND,MJ,
C     *              FIGMIN,FIGMAX,SS,K,LOC )
      STOP
      E N D
      SUBROUTINE  DEFALT
C
C  ...  THIS SUBROUTINE SETS DEFAULT VALUES OF PARAMETERS  ...
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Y(500)
      COMMON   /COMAAA/  Y
      COMMON   /C91214/  XMIN, XMAX, FIGMIN, FIGMAX, YM(10)
      COMMON   /C91215/  NOISEV, NOISEW, INITD, ITF, ITH
      COMMON   /C91216/  B, OUTMIN, OUTMAX
      COMMON   /C91219/  N, KK, NS, NFE, NPE
      NOISEV = 2
      NOISEW = 1
      INITD  = 1
      ITF    = 1
      ITH    = 1
      B = 1.00D0
      OUTMIN = -1.0D30
      OUTMAX =  1.0D30
      CALL  MAXMIN( Y,N,XMIN,XMAX,DY )
      NS = 1
      NFE = N
      NPE = N
      FIGMIN = XMIN
      FIGMAX = XMAX
C
      RETURN
      E N D
      SUBROUTINE  NGSMTH( Y,P,F,S,T,N,K,DX,XMIN,Q,FF,PS,FS,SS,LOC )
C
C  ...  NON-GAUSSIAN SMOOTHER  ...
C
C     INPUTS:
C       Y(I):   TIME SERIES
C       P(I):   INITIAL DENSITY
C       N:      DATA LENGTH
C       K:      NUMBER OF INTERVALS IN STEP FUNCTION APPROXIMATION
C       DX:     WIDTH OF INTERVAL
C       XMIN:   MINIMUM OF THE INTERVAL
C       Q:      SYSTEM NOISE DENSITY
C     OUTPUTS:
C       FF:     LOG-LIKELIHOOD
C       LOC(I): LOCATION OF THE CENTER OF THE INTERVAL AT STEP I
C       SS:     SMOOTHED DENSITY
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  P(K), F(K), S(K), T(K), Y(N), Q(-K:K), LOC(N)
      REAL*4     FS(K,N), PS(K,N), SS(K,N)
      COMMON   /C91215/  NOISEV, NOISEW, INITD, ITF, ITH
      COMMON   /C91216/  B, OUTMIN, OUTMAX
      COMMON   /C91219/  NN, KK, NS, NFE, NPE
C
      FF = 0.0D0
C
      DO 200 II=NS,NPE
C
C  ...  CONVOLUTION (SYSTEM NOISE)  ...
C
      CALL  CONVOL( Q,F,K,P )
      CALL  NORMLZ( P,K,DX,PSUM )
C
C  ...  BAYES FORMULA  ...
C
      IF( Y(II).LE.OUTMIN .OR. Y(II).GE.OUTMAX .OR. II.GT.NFE ) THEN
      DO 110 I=1,K
  110 F(I) = P(I)
      ELSE
      CALL  BAYES( P,K,XMIN,DX,Y(II),F,LOC(II) )
      CALL  NORMLZ( F,K,DX,FINT )
C
C  ...  LIKELIHOOD COMPUTATION  ...
C
      FF = FF + DLOG( FINT )
      IF( MOD(II,10).EQ.0)  WRITE(6,*) II,FF
      END IF
C
C  ...  SAVE FOR SMOOTHING  ...
C
      DO 130 I=1,K
      PS(I,II) = P(I)
  130 FS(I,II) = F(I)
C
C  ...  SHIFT ORIGIN  ...
C
      CALL  SHIFT( F,K,T,II,N,LOC )
C
  200 CONTINUE
C
C  ...  SMOOTHING  ...
C
      DO 190 J=NFE,NPE
      DO 190 I=1,K
  190 SS(I,J) = FS(I,J)
      DO 195 I=1,K
  195 S(I) = FS(I,NFE)
C
      DO 300 II=NFE-1,NS,-1
      IF(MOD(II,10).EQ.0)  WRITE(6,*) II
      DO 210 I=1,K
      T(I) = 0.0D0
      P(I) = 0.0D0
  210 F(I) = FS(I,II)
      DO 220 I=1,K
      J = I - (LOC(II+1)-LOC(II))
      IF( J.GE.1.AND.J.LE.K )  P(I) = PS(J,II+1)
  220 IF( J.GE.1.AND.J.LE.K )  T(I) = S(J)
      DO 230 I=1,K
  230 S(I) = T(I)
C
      CALL  SCONVL( Q,S,P,F,K,T )
      CALL  NORMLZ( T,K,DX,TSUM )
C
      DO 240 I=1,K
      S(I) = T(I)
  240 SS(I,II) = S(I)
  300 CONTINUE
C
      RETURN
      E N D
      SUBROUTINE  PINTVL( P,K,XMIN,DX,Y )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  P(K), Y(7), PROB(7), P1(401)
      DATA  PROB /0.0013D0, 0.0227D0, 0.1587D0, 0.5000D0, 0.8413D0,
     *            0.9773D0,0.9987D0/
C
      P1(1) = 0.0
      DO 10 I=2,K
   10 P1(I) = P1(I-1) + (P(I-1) + P(I))*DX/2.0
C
      DO 30 J=1,7
      PP = PROB(J)
      DO 20 I=2,K
      IF(P1(I-1).LE.PP .AND. P1(I).GT.PP)  GO TO 30
   20 CONTINUE
   30 Y(J) = XMIN + (I-2)*DX + DX*(PP - P1(I-1))/(P1(I) - P1(I-1))
C
      RETURN
      E N D
      SUBROUTINE  CONVOL( Q,S,K,P )
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  S(K), P(K), Q(-K:K)
C
      DO 20 I=1,K
      J1 = 1-I
      J2 = K-I
      SUM = 0.0
      DO 10 J=J1,J2
   10 SUM = SUM + S(I+J)*Q(J)
   20 P(I) = SUM
C
      RETURN
      E N D
      SUBROUTINE  SCONVL( Q,P,R,S,K,T )
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  S(K), P(K), R(K), T(K), Q(-K:K)
C
      DO 20 I=1,K
      J1 = 1-I
      J2 = K-I
      SUM = 0.0D0
      DO 10 J=J1,J2
   10 IF(P(I+J).GT.0.0D0)  SUM = SUM + P(I+J)/R(I+J)*Q(J)
   20 T(I) = S(I)*SUM
      RETURN
      E N D
      SUBROUTINE  TRANS( FUNCT,K,DX,TAU2,BV,Q )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Q(-K:K), PARAM(3)
      EXTERNAL   FUNCT
C
      PARAM(1) = 0.0D0
      PARAM(2) = TAU2
      PARAM(3) = BV
C
      DO 20 I=1-K,K-1
      X0 = -DX*I - DX/2
      SUM = (FUNCT(X0,PARAM) + FUNCT(X0+DX,PARAM))/2
      DO 10 J=1,49
      X = X0 + (DX*J)/50
   10 SUM = SUM + FUNCT( X,PARAM )
   20 Q(I) = SUM*DX/50
C
      RETURN
      E N D
      SUBROUTINE  NORMLZ( P,K,DX,SUM )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  P(K)
C
      SUM = 0.0D0
      DO 10 I=1,K
   10 SUM = SUM + P(I)
      SUM = SUM*DX
      DO 30 I=1,K
   30 P(I) = P(I)/SUM
C
      RETURN
      E N D
      SUBROUTINE  IDIST( P,K,P1,P2,XMIN,DX )
      IMPLICIT REAL*8( A-H,O-Z )
      DIMENSION  P(K), PARAM(3)
      COMMON   /C91215/  NOISEV, NOISEW, INITD, ITF, ITH
C
      PARAM(1) = P1
      PARAM(2) = P2
C
      DO 10 I=1,K
      X = XMIN + DX*(I-1)
      IF( INITD.EQ.0 )  P(I) = USERI( X,PARAM )
      IF( INITD.EQ.1 )  P(I) = GAUSS( X,PARAM )
      IF( INITD.EQ.2 )  P(I) = UNIF ( X,PARAM )
   10 CONTINUE
      RETURN
      END
      SUBROUTINE  BAYES( P,K,XMIN,DX,Y,F,LSHIFT )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  P(K), F(K), PARAM(3)
      COMMON   /C91215/  NOISEV, NOISEW, INITD, ITF, ITH
      COMMON   /C91218/  SIG2, TAU2, BW, BV
      EXTERNAL  USERW
      EXTERNAL  GAUSS
      EXTERNAL  PEARSN
      EXTERNAL  TWOEXP
      EXTERNAL  DBLEXP
C
      PARAM(2) = SIG2
      PARAM(3) = BW
C
      DO 10 I=1,K
      PARAM(1) = XMIN + DX*(I-1+LSHIFT)
      IF( NOISEW.EQ.0 )  F(I) = P(I)*USERW ( Y,PARAM )
      IF( NOISEW.EQ.1 )  F(I) = P(I)*GAUSS ( Y,PARAM )
      IF( NOISEW.EQ.2 )  F(I) = P(I)*PEARSN( Y,PARAM )
      IF( NOISEW.EQ.3 )  F(I) = P(I)*TWOEXP( Y,PARAM )
      IF( NOISEW.EQ.4 )  F(I) = P(I)*DBLEXP( Y,PARAM )
   10 CONTINUE
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  USERW( Y,PARAM )
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PARAM(3)
      DATA  C1  /2.506628275D0/
C
      YMEAN = 0.0D0
      VAR   = PARAM(1)
      USERW = DEXP( -(Y-YMEAN)**2/(2*VAR) )/(C1*DSQRT( VAR ))
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  USERV( X,PARAM )
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PARAM(3)
      DATA  C1  /2.506628275D0/
C
      USERV = DEXP( -(X-PARAM(1))**2/(2*PARAM(2)) )
     *                /(C1*DSQRT( PARAM(2) ))
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  TWOEXP( X,PARAM )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  PARAM(2)
C
      TWOEXP = DEXP( -DABS(X-PARAM(1))*PARAM(2) )*PARAM(2)/2.0D0
      RETURN
C
      E N D
      SUBROUTINE  SHIFT( F,K,T,II,N,LOC )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  F(K), T(K), LOC(*)
C
C  ...  FIND THE POSTERIOR MODE AND SHIFT ORIGIN  ...
C
      PMAX = 0.0D0
      DO 10  I=1,K
      IF( F(I).LE.PMAX )  GO TO 10
         PMAX = F(I)
         IMAX = I
   10 CONTINUE
      IF(II.LT.N)  LOC(II+1) = LOC(II) + IMAX - (K+1)/2
      DO 20 I=1,K
      J = I + IMAX - (K+1)/2
      T(I) = 0.0D0
   20 IF(J.GE.1.AND.J.LE.K )  T(I) = F(J)
      DO 30 I=1,K
   30 F(I) = T(I)
C
      RETURN
      E N D
      SUBROUTINE  PRNGSM( NOISEV,NOISEW,TAU2,SIG2,FF,TREND,MJ,N )
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER  TITLE*72
      DIMENSION  TREND(MJ,7)
      COMMON   /CMDATA/  TITLE
C
      WRITE(6,600)
      WRITE(6,610) TITLE
      WRITE(6,620) NOISEV, NOISEW
      WRITE(6,630) TAU2, SIG2, FF
      WRITE(6,640)
      WRITE(6,650)  (TREND(I,4),I=1,N)
      RETURN
  600 FORMAT( 1H1,'PROGRAM 14.1:  NON-GAUSSIAN SMOOTHING' )
  610 FORMAT( 1H ,A72 )
  620 FORMAT( 1H ,'NOISEV =',I2,3X,'NOISEW =',I2 )
  630 FORMAT( 1H ,'TAU1 =',D12.5,3X,'SIG2 =',D12.4,3X,'LOG-LIKELIHOOD ='
     *           ,F15.5 )
  640 FORMAT( 1H ,'***  ESTIMATED  TREND  ***' )
  650 FORMAT( 1H ,5D15.6 )
      E N D
      SUBROUTINE  READTS( IDEV,Y,N )
      REAL*8        Y(N)
      CHARACTER*72  TITLE
      COMMON  /CMDATA/  TITLE
C
      OPEN( IDEV,FILE='temp.dat' )
      READ( IDEV,1 )  TITLE
      READ( IDEV,* )  N
      READ( IDEV,* )  (Y(I),I=1,N)
      CLOSE( IDEV )
      RETURN
    1 FORMAT( A72 )
      E N D
      SUBROUTINE  MOMENT( Y,N,YMEAN,VAR )
C
C  ...  Mean and variance of the data  ...
C
C     Inputs:
C       Y(I):   data
C       N:      data length
C     Outputs:
C       YMEAN:  mean
C       VAR:    variance
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  Y(N)
C
      SUM = 0.0D0
      DO 10 I=1,N
   10 SUM = SUM + Y(I)
      YMEAN = SUM/N
      SUM = 0.0D0
      DO 20 I=1,N
   20 SUM = SUM + (Y(I)-YMEAN)**2
      VAR = SUM/N
C
      RETURN
      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  USERI( X,PARAM )
C
C  ...  User supplied density function  ...
C       (The following is an example of two-sided exponential dist.)
C
C     Inputs:
C        X:
C        PARAM(1):  mean
C        PARAM(2):  lambda
C     Output:
C        USERI:     density at X
C
      IMPLICIT  REAL*8(A-H,O-Z)
      DIMENSION  PARAM(2)
C
      SIGMA = DSQRT( PARAM(2) )
         USERI = SIGMA*DEXP( -SIGMA*DABS(X-PARAM(1)) )/2
      RETURN
      E N D
      DOUBLE PRECISION FUNCTION  UNIF( X,PARAM )
C
C  ...  uniform distribution  f(x) = 1  ...
C
C     Output:
C        UNIF:    density at X (=1)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      UNIF = 1.0D0
      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