*    
*     bayes-gbm.f  
*     (Bayes estimates for a micromovement model with 
*      ^^^^^               
*     Geometric Brownian Motion as the value process)
*     ^         ^        ^
*     Copyright (C) 2003, 2005
*
*     Yong Zeng
*     Department of Mathematics and Statistics
*     University of Missouri at Kansas City
*     5100 Rockhill Road
*     Kansas City, MO 64110-2499 USA.
*
*     Permission to use, copy, modify, and distribute this
*     software for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*----------------------------------------------------------------
C     the program computes the joint posterior and the Bayes estimates 
c     for the simple model described
C     in the paper Yong Zeng (2003) "A Partially-Observed Model for 
c     Micro-movement of Asset Prices with Bayes Estimation via Filtering", 
c     Mathematical Finance,
c     (2003), Vol. 13, pp. 411-444.
*----------------------------------------------------------------
c     "par.f" specifies the parameters of this program
c     simulated date set: smlt9403.dat (it has two columes: price, and 
c     duration or waiting time between trades)
c     For a different data set such as msft9403tmdff.dat (MicroSoft,
c     March, 1994), you need to reset the ranges for X, mu, sigma and rho.
c     All files should be in the same directory
c     > f77 bayes-gbm.f   (to generate "a.out")
c     > a.out  (to run)
c     Warning: it may take quite a long time to finish the program.
C----------------------------------------------------------------
C     PX0(N): VECTOR OF INITIAL DIST OF X_0
C     PM(NM): PRIOR OF MU,THEN POSTERIOR
C     PS(NS): PRIOR OF SIGMA,THEN POSTERIOR
C     
      PROGRAM DONDDIST
      include 'par-gbm.f'
      INTEGER I,K,L,M,W
      REAL*8 PX0(N),PM(NM),PS(NS),PR(NR),P(N,NM,NS,NR),MU(NM),SG(NS),
     &RO(NR)
c     ASSIGN INITIAL DISTRIBUTION OF X_0 ---------------
      DO 10 I=1,N
         PX0(I)=0.0d0
 10   CONTINUE
      PX0(I0)=1.0d0
C     --------------------------------------------------
C     ASSIGN PRIOR FOR MU, PM(NM) ----------------------
      DO 20 L=1,NM
         PM(L)=1.0d0/dble(NM)
 20   CONTINUE
C     --------------------------------------------------
C     ASSIGN PRIOR FOR SIGMA, PS(NS) -------------------
      DO 30 M=1,NS
         PS(M)=1.0d0/dble(NS)
 30   CONTINUE
C     --------------------------------------------------
C     ASSIGN PRIOR FOR RHO, PR(NR)
      DO 40 W=1,NR
         PR(W)=1.0D0/DBLE(NR)
 40   CONTINUE
c     COMPUTE INITIAL DIST OF P(X_0,MU,SIGMA | Z_0) ----
      DO 600 K=1,N
         DO 400 L=1,NM
            DO 200 M=1,NS
            DO 190 W=1,NR
                  P(K,L,M,W)=PX0(K)*PM(L)*PS(M)*PR(W)
 190           CONTINUE
 200           CONTINUE
 400           CONTINUE
 600           CONTINUE
C     ------------------------------------------------
      CALL DONDDIST3(P)
C     --------------------------------------------------------------------
C     COMPUTE THE POSTERIOR OF MU-----------------------
      OPEN(UNIT=30,FILE="12MU.LST",STATUS="UNKNOWN")
      DO 900 L=1,NM
        PM(L)=0d0
        MU(L)=MUMIN+(dble(L)-1.000d0)*MUL
        DO 930 K=1,N
           DO 940 W=1,NR
           DO 950 M=1,NS
           PM(L)=PM(L)+P(K,L,M,W)
 950       CONTINUE
 940       CONTINUE
 930    CONTINUE
         WRITE(30,'(g16.6E4,G16.6E4)') MU(L),PM(L)
 900  CONTINUE
      CLOSE(30)
C     -----------------------------------------------------
C     COMPUTE THE POSTERIOR OF SIGMA-----------------------
      OPEN(UNIT=20,FILE="12SG.LST",STATUS="UNKNOWN")
      DO 1000 M=1,NS
         PS(M)=0d0
         SG(M)=SGMIN+(dble(M)-1.00d0)*SGL
         DO 1030 K=1,N
            DO 1040 W=1,NR
            DO 1050 L=1,NM
               PS(M)=PS(M)+P(K,L,M,W)
 1050       CONTINUE
 1040       CONTINUE
 1030    CONTINUE
         WRITE(20,'(G16.6E4,G16.6E4)') SG(M),PS(M)
 1000 CONTINUE
      CLOSE(20)
C     COMPUTE THE POSTERIOR OF RHO-----------------------
      OPEN(UNIT=15,FILE="12RO.LST",STATUS="UNKNOWN")
      DO 1100 W=1,NR
         PR(W)=0d0
         RO(W)=ROMIN+(dble(W)-1.00d0)*ROL
         DO 1130 K=1,N
            DO 1140 M=1,NS
            DO 1150 L=1,NM
               PR(W)=PR(W)+P(K,L,M,W)
 1150       CONTINUE
 1140       CONTINUE
 1130    CONTINUE
         WRITE(15,'(G16.6E4,G16.6E4)') RO(W),PR(W)
 1100 CONTINUE
      CLOSE(15)
C     -----------------------------------------------------
      END

      SUBROUTINE DONDDIST3(P)
      EXTERNAL LD1,LD2,LD3
      INCLUDE 'par-gbm.f'
      INTEGER I,K,L,M,NTI,J,W
      REAL*8 P(N,NM,NS,NR),A(N),B(N),X2,DT
      REAL*8 MU,SG,RO,XB,XC,LD1,LD2,LD3
      REAL*8 PI1(N),PI2(N)
c     'midout.flt' is for reporting middle results
      OPEN(UNIT=32,FILE='midout.flt',STATUS="UNKNOWN")
c     read in and store data in X2 and DT one by one
      OPEN(UNIT=1,FILE='smlt9403.dat', STATUS='OLD')
c     loop 6: updating the cond prob from time t_1 to t_ndata
      DO 6 I=1,NDATA
         READ(1,*) X2,DT
c     loop 100: approXimate p(K,L,M,W,t_i-)
         DO 100 M=1,NS
            SG=SGMIN+(dble(M)-1.00d0)*SGL
            DO 110 L=1,NM
               MU=MUMIN+(dble(L)-1.00000d0)*MUL
C     ASSIGN A(K), B(K) ACCORDINGLY------------------------------------
            IF (MU .GE. 0D0) THEN
               DO 120 K=1,N
         A(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & + MU*(DBLE(K)-1.0D0+(XMIN-LL)*32.0D0)
         B(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
 120           CONTINUE
            ELSE 
               DO 125 K=1,N
         A(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
         B(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & - MU*(DBLE(K)-1.0D0+(XMIN-LL)*32.0D0)
 125           CONTINUE
            END IF
c     -----------------------------------------------------
c     refine the partition based on the interarrival time 
c     -----------------------------------------------------
      DO 133 W=1,NR
C         RO=ROMIN+(DBLE(W)-1.00D0)*ROL
         DO 50 K=1,N
            PI2(K)=P(K,L,M,W)
 50      CONTINUE
         NTI=IDNINT(DT/LT)+1
         DO 55 J=1,NTI
            DO 58 K=1,N
               PI1(K)=PI2(K)
 58         CONTINUE
                DO 130 K=1,N
            IF (K .EQ. 1) THEN
               PI2(K)=PI1(K)+(-(A(K)+B(K))*PI1(K)
     &+B(K+1)*PI1(K))*DT/DBLE(NTI)
            ELSE IF (K .EQ. N) THEN
               PI2(K)=PI1(K)+(A(K-1)*PI1(K-1)-(A(K)+B(K))
     &*PI1(K))*DT/DBLE(NTI)
            ELSE
               PI2(K)=PI1(K)+(A(K-1)*PI1(K-1)-(A(K)+B(K))
     &*PI1(K)+B(K+1)*PI1(K+1))*DT/DBLE(NTI)
            END IF
 130            CONTINUE
 55         CONTINUE
            DO 60 K=1,N
               P(K,L,M,W)=PI2(K)
 60         CONTINUE
 133        CONTINUE
 110     CONTINUE
 100     CONTINUE
c     ----------------------------------------------------------------
C     "IF STRUCTURE": ASSIGN FOR DIFFERENT INTENSITY FUNCTIONS TO OBTAIN
C     P(T_I)
       XB=X2*2D0
       XC=X2*4D0
C     CASE 1: X2 IS AN INTEGER OR HALF
         IF (DMOD(XB,1d0) .EQ. 0) THEN
            CALL DEND(LD1,P,X2)
C     CASE 2: X2 IS ODD QUARTERS
         ELSE IF (DMOD(XC,1d0) .EQ. 0) THEN
            CALL DEND(LD2,P,X2)
C     CASE 3: X2 IS ODD EIGHTH
         ELSE
            CALL DEND(LD3,P,X2)
         END IF
         if ((i .eq. 1) .or. (i .eq. 5000) .or. (i .eq. 10000) 
     & .or. (i .eq.15000) .or. (i .eq. 11213) 
     & .or. (i .eq. 20000) .or. (i .eq. 25000) .or. (i .eq. 30000)
     & )  then 
            CALL WOUTPUT(P,I,X2,DT)
         end if
 6    continue
      close(1)
      CLOSE(32)
         RETURN
         END

      SUBROUTINE DEND(FUNC,P,X)
      include 'par-gbm.f'
      INTEGER K,L,M,W
      REAL*8 FUNC,P(N,NM,NS,NR),P1(N,NM,NS,NR),X,XA,TL
C     PUT P TO P1
         DO 80 K=1,N
            DO 82 L=1,NM
               DO 86 M=1,NS
                  DO 88 W=1,NR
                 P1(K,L,M,W)=P(K,L,M,W)
 88              CONTINUE
 86            CONTINUE
 82         CONTINUE
 80      CONTINUE
         TL=0D0
         DO 180 L=1,NM
            DO 182 M=1,NS
               DO 184 W=1,NR
               DO 186 K=1,N
                  XA=(DBLE(K)-1D0+(XMIN-LL)*32D0)/32D0 
                 TL=TL+P1(K,L,M,W)*FUNC(X,XA,W)
 186            CONTINUE
 184            CONTINUE
 182         CONTINUE
 180      CONTINUE
          IF (TL .EQ. 0D0) THEN 
             PAUSE "TL=0"
          END IF
         DO 280 L=1,NM
            DO 282 M=1,NS
               DO 284 W=1,NR
               DO 286 K=1,N
                  XA=(DBLE(K)-1D0+(XMIN-LL)*32D0)/32D0 
                 P(K,L,M,W)=P1(K,L,M,W)*FUNC(X,XA,W)/TL
 286            CONTINUE
 284            CONTINUE
 282         CONTINUE
 280      CONTINUE
      RETURN
      END

      
      SUBROUTINE WOUTPUT(P,I,X2,DT)
      INCLUDE 'par-gbm.f'
      INTEGER I,K,L,M,W
      REAL*8 PM(NM),PS(NS),PR(NR),P(N,NM,NS,NR),MU(NM),SG(NS),
     &RO(NR),X2,DT,EMU,EMU2,SMU,ESG,ESG2,SSG,ERO,ERO2,SRO,EX,EX2,
     &SX,PX(N),X(N)
      EMU=0D0
      EMU2=0D0
      SMU=0D0
      ESG=0D0
      ESG2=0D0
      SSG=0D0
      ERO=0D0
      ER02=0D0
      SRO=0d0
      EX=0D0
      EX2=0D0
      SX=0D0
      WRITE(32,*) "I,X2,DT",I,X2,DT
      WRITE(32,*) "posterior of MU"
      DO 901 L=1,NM
        PM(L)=0d0
        MU(L)=MUMIN+(dble(L)-1.000d0)*MUL
        DO 931 K=1,N
           DO 941 W=1,NR
           DO 951 M=1,NS
           PM(L)=PM(L)+P(K,L,M,W)
 951    CONTINUE
 941  CONTINUE
 931  CONTINUE
         WRITE(32,'(g16.6E4,G16.6E4)') MU(L),PM(L)
         EMU=EMU+MU(L)*PM(L)
         EMU2=EMU2+MU(L)*MU(L)*PM(L)
 901  CONTINUE
        SMU=DSQRT(EMU2-EMU*EMU)
      WRITE(32,*) "BAYESIAN ESTIMATE OF MU AND ITS SD",EMU,SMU
      WRITE(32,*) "POSTERIOR OF SG"
      DO 1001 M=1,NS
         PS(M)=0d0
         SG(M)=SGMIN+(dble(M)-1.00d0)*SGL
         DO 1031 K=1,N
            DO 1041 W=1,NR
            DO 1051 L=1,NM
               PS(M)=PS(M)+P(K,L,M,W)
 1051       CONTINUE
 1041    CONTINUE
 1031 CONTINUE
         WRITE(32,'(G16.6E4,G16.6E4)') SG(M),PS(M)
        ESG=ESG+SG(M)*PS(M)
         ESG2=ESG2+SG(M)*SG(M)*PS(M)
 1001 CONTINUE
        SSG=DSQRT(ESG2-ESG*ESG)
      WRITE(32,*) "BAYESIAN ESTIMATE OF SG AND ITS SD",ESG,SSG
         WRITE(32,*) "POSTEROR OF RHO"
      DO 1101 W=1,NR
         PR(W)=0d0
         RO(W)=ROMIN+(dble(W)-1.00d0)*ROL
         DO 1131 K=1,N
            DO 1141 M=1,NS
            DO 1151 L=1,NM
               PR(W)=PR(W)+P(K,L,M,W)
 1151       CONTINUE
 1141    CONTINUE
 1131 CONTINUE
         WRITE(32,'(G16.6E4,G16.6E4)') RO(W),PR(W)
        ERO=ERO+RO(W)*PR(W)
         ERO2=ERO2+RO(W)*RO(W)*PR(W)
 1101 CONTINUE
        SRO=DSQRT(ERO2-ERO*ERO)
      WRITE(32,*) "BAYESIAN ESTIMATE OF RO AND ITS SD",ERO,SRO
         WRITE(32,*) "POSTERIOR OF X"
        EX=0D0
         EX2=0D0
      DO 1201 K=1,N
         PX(K)=0d0
         X(K)=XMIN-LL+(dble(K)-1.00d0)*0.031250d0
         DO 1231 W=1,NR
            DO 1241 M=1,NS
            DO 1251 L=1,NM
               PX(K)=PX(K)+P(K,L,M,W)
 1251       CONTINUE
 1241    CONTINUE
 1231 CONTINUE
        EX=EX+X(K)*PX(K)
         EX2=EX2+X(K)*X(K)*PX(K)
         IF ( PX(K) .GT. 0.01d0) THEN
           WRITE(32,'(G16.6E4,G16.6E4)') X(K), PX(K)
         END IF
 1201 CONTINUE
        SX=DSQRT(EX2-EX*EX)
      WRITE(32,*) "BAYESIAN ESTIMATE OF X AND ITS SD",EX,SX
      RETURN
      END

      FUNCTION LD1(X,Y,W)
      INCLUDE 'par-gbm.f'
      REAL*8 X,LD1,Y,RHO
      INTEGER XA,W
c     add 0.000001 for protection of roundoff error
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=ROMIN+(DBLE(W)-1.0D0)*ROL
      IF (XA .EQ. 0) THEN
         LD1=(1.0D0-RHO)*(1.0D0+BETA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD1=0.50d0*(1.0d0-rho)*rho*(1.0d0+rho*beta)
     &+(1.0d0-rho)*beta
      ELSE
         LD1=0.50d0*(1.0D0-RHO)*(RHO**(IABS(XA)))*
     & (1.0D0+BETA*(RHO+1.0d0/RHO))
      END IF
      RETURN
      END

      FUNCTION LD2(X,Y,W)
      INCLUDE 'par-gbm.f'
      REAL*8 X,LD2,Y,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=ROMIN+(DBLE(W)-1.0D0)*ROL
      IF (XA .EQ. 0) THEN
         LD2=(1.0D0-RHO)*(1.0D0+ALPHA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD2=0.50d0*(1.0d0-rho)*rho*(1.0d0+rho*alpha)
     &+(1.0d0-rho)*alpha
      ELSE
         LD2=0.50D0*(1.0D0-RHO)*(RHO**(IABS(XA)))*
     & (1.0D0+ALPHA*(RHO+1/RHO))
      END IF
      RETURN
      END
      
      FUNCTION LD3(X,Y,W)
      INCLUDE 'par-gbm.f'
      REAL*8 X,Y,LD3,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=ROMIN+(DBLE(W)-1.0D0)*ROL
      IF (XA .EQ. 0) THEN
        LD3=(1.0D0-RHO)*(1.0D0-ALPHA-BETA)
      ELSE
        LD3=0.50D0*(1.0D0-ALPHA-BETA)*(1.0D0-RHO)*(RHO**(IABS(XA))) 
      END IF 
      RETURN
      END