*    
*     bayes-gbm.f  
*     (Bayes estimates for a micromovement model with 
*      ^^^^^               
*     Geometric Brownian Motion as the value process)
*     ^         ^        ^
*     Copyright (C) 2003, 2005, 2007
*
*     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     "par12.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     Put all files (bayes-gbm.f, par12.f and smlt9403) in the same directory
c
c     > f77 bayes-gbm.f   (to generate "a.out")
c     > a.out  (to run)
c
c     output files: final posteriors of MU, SG, and RHO: 12MU.LST, 12SG.LST
c     12RO.LST; and midout.flt (for middel outputs). 
c     Warning: it may take quite a long time to finish the program.

C     The model with alpha, beta and gamma of 1/8
C     PROGRAM FOR COMPUTING THE BAYESIAN FILTERING ESTIMATES 
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 'par12.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     ,PMSR(NM,NS,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      write(*,*) "pr", pr
c      pause "pr"
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      WRITE(*,*) "END CALL DONDDIST3"
c      pause "end call main"
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      WRITE(*,*) "END LOOP 900"
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,LD0
      INCLUDE 'par12.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,LD0
      REAL*8 PI1(N),PI2(N)
c     store data in X(I) and DT
      OPEN(UNIT=32,FILE='midout12.flt',STATUS="UNKNOWN")
      OPEN(UNIT=1,FILE='smlt9403.dat',
     &STATUS='OLD')
      DO 6 I=1,NDATA
         READ(1,*) X2,DT
c         write(*,*) X2,dt
c     loop 6: updating the cond prob from time t_1 to t_ndata
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               write(*,*) mu,SG
C     ASSIGN A(K), B(K) ACCORDINGLY------------------------------------
            IF (MU .GE. 0D0) THEN
c               pause "mu ge 0"
               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            write(*,*) "end if of A B",A,B
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
c         print*,lt,nti
         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
c               pause "k=1"
               PI2(K)=PI1(K)+(-(A(K)+B(K))*PI1(K)
     &+B(K+1)*PI1(K))*DT/DBLE(NTI)
            ELSE IF (K .EQ. N) THEN
c               pause "k=n"
               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
            if (pi2(k) .lt. 0d0) then
          write(*,*) "sg",sg,"mu",mu,"RO",RO
               write(*,*) "i,k,l,m,W",i,k,l,m,W,"DT",DT
               write(*,*) "a_k,b_k,a_k-1,b_k+1",a(k),b(k),a(k-1),b(k+1)
               write(*,*) "PI1(K-1),PI1(K),PI1(K+1)",
     &PI1(K-1),PI1(K),PI1(K+1)
               pause "negative PROBABILITY"
            end if
 130            CONTINUE
 55         CONTINUE
            DO 60 K=1,N
               P(K,L,M,W)=PI2(K)
 60         CONTINUE
 133        CONTINUE
c          write (*,*) 'end loop 133'
 110     CONTINUE
 100     CONTINUE
c         if (i .eq. 6) then
c         pause "i=6,end loop 100"
c         end if
c         IF (I .EQ. 1) THEN
c            DO 450 L=1,NM
cc               DO 452 M=1,NS
c                  DO 454 W=1,NR
c                  DO 456 K=1,N
c                    IF (P(K,L,M,W) .GT. 0D0) THEN
c                       WRITE(32,*) K,L,M,W,P(K,L,M,W)
c                     END IF
c                        XA=(DBLE(K)+(XMIN-LL)*8D0-1D0)/8D0 
c                 WRITE(*,*) "XA",XA
c                 PAUSE
c                 IF (X2 .EQ. XA) THEN
c                    WRITE(*,*) "XA",XA,X2,p(k,l,m,W),K,L,M,I
c                    PAUSE
c                 END IF
c 456              CONTINUE
c 454              CONTINUE
c 452           CONTINUE
c 450        CONTINUE
c         END IF
c         pause
c     ----------------------------------------------------------------
C     "IF STRUCTURE": ASSIGN FOR DIFFERENT INTENSITY FUNCTIONS TO OBTAIN
C     P(T_I)
       XB=X2*2D0
       XC=X2*4D0
c       write (*,*) 'Xb',xb,XC
C     CASE 1: X2 IS AN INTEGER OR HALF
         IF (DMOD(X2,1d0) .EQ. 0) THEN
C       write (*,*) "inside if"
            CALL DEND(LD0,P,X2)
C     CASE 2: X2 is an half
         else if (DMOD(XB,1d0) .EQ. 0) THEN
            CALL DEND(LD1,P,X2)
C     CASE 3: X2 IS ODD QUARTERS
         ELSE IF (DMOD(XC,1d0) .EQ. 0) THEN
C            write(*,*) "INSIDE else if"
            CALL DEND(LD2,P,X2)
C     CASE 4: X2 IS ODD EIGHTH
         ELSE
C            write (*,*) "INSIDE else"
            CALL DEND(LD3,P,X2)
         END IF
c         if ((i .eq. 6) .or. (i .eq. 7) .or. ((i .ge. 10) .and. 
c     & (i .le. 34))) then
c         if ( (i .ge. 87) .and. (i .le. 109)) then
         if ((i .eq. 1) .or. (i .eq. 5000) .or. (i .eq. 10000) 
     & .or. (i .eq.15000) .or. (i .eq. 32347) .or. (i .eq. 32500)
     & .or. (i .eq. 20000) .or. (i .eq. 25000) .or. (i .eq. 30000)
     & )  then 
            CALL WOUTPUT(P,I,X2,DT)
         end if
c         WRITE(*,*) "END IF STRUCTURE"
C         If (I .lt. 4) then
C         write(10,*) "I",I, "p(t_i)", P
C         end if
C         if (I .gt. 995) then
C            write(10,*) I, "p(t_i)", P
C         end if
 6    continue
      close(1)
      CLOSE(32)
         RETURN
         END

      SUBROUTINE DEND(FUNC,P,X)
      include 'par12.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
c         write(*,*) "END LOOP 80"
c         write(10,*) "loop 80,p1",p1
C     SUM UP THE TOTEL TL
c         write(*,*) "func", x, func(x, 78d0)
         TL=0D0
         DO 180 L=1,NM
            DO 182 M=1,NS
               DO 184 W=1,NR
               DO 186 K=1,N
C                 XA=(DBLE(K)-1D0+(XMIN-LL)*8D0)/8D0  (FOR 1/8)
C                 WRITE(*,*) P1(K,L,M),FUNC(X,XA)
C      FOR 1/16
                  XA=(DBLE(K)-1D0+(XMIN-LL)*32D0)/32D0 
                 TL=TL+P1(K,L,M,W)*FUNC(X,XA,W)
c                     IF (P1(K,L,M,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,x,xa,P1(K,L,M,W)
c                      write(32,*) func(x,xa,w),func(x,xa,w)*p1(k,l,m,w)
c                      END IF
 186            CONTINUE
 184            CONTINUE
 182         CONTINUE
 180      CONTINUE
          IF (TL .EQ. 0D0) THEN 
             PAUSE "TL=0"
          END IF
c         write(32,*) "END LOOP 180", "tl",tl
         DO 280 L=1,NM
            DO 282 M=1,NS
               DO 284 W=1,NR
               DO 286 K=1,N
c                 XA=(DBLE(K)-1D0+(XMIN-LL)*8D0)/8D0  (for 1/8)
                  XA=(DBLE(K)-1D0+(XMIN-LL)*32D0)/32D0 
                 P(K,L,M,W)=P1(K,L,M,W)*FUNC(X,XA,W)/TL
c                     IF (P(K,L,M,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,xa,P(K,L,M,W)
c                     END IF
 286            CONTINUE
 284            CONTINUE
 282         CONTINUE
 280      CONTINUE
c         write(*,*) "END LOOP 280"
      RETURN
      END

      
      SUBROUTINE WOUTPUT(P,I,X2,DT)
      INCLUDE 'par12.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
      ERO2=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 'par12.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
C      WRITE(*,*) "INSIDE LD1"
      IF (XA .EQ. 0) THEN
         LD1=(1.0D0-RHO)*(1.0D0+BETA*RHO*(1+RHO**2))
      ELSE IF (IABS(XA) .eq. 1) Then
         LD1=0.50d0*(1.0d0-rho)*(rho+beta*(2+2*RHO**2
     &+RHO**4))
      ELSE IF (IABS(XA) .eq. 2) Then
         LD1=0.50d0*(1.0D0-RHO)*RHO*(rho+beta*(2+rho**2
     &+rho**4))
      ELSE IF (IABS(XA) .eq. 3) Then
         LD1=0.50d0*(1.0D0-RHO)*(RHO**3+beta*(2+rho**2
     &+rho**4+rho**6))
      ELSE
         LD1=0.50d0*(1.0D0-RHO)*(rho**(IABS(XA)-3))
     &*(RHO**3+beta*(1+rho**2+rho**4+rho**6))
      END IF
C      write(*,*) X,Y,W,XA,RHO, ROMIN,ROL,BETA,LD1
C      PAUSE
      RETURN
      END

      FUNCTION LD2(X,Y,W)
      INCLUDE 'par12.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 'par12.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-GAMMA)
      ELSE
        LD3=0.50D0*(1.0D0-ALPHA-BETA-GAMMA)*(1.0D0-RHO)
     &*(RHO**(IABS(XA))) 
      END IF 
      RETURN
      END

      FUNCTION LD0(X,Y,W)
      INCLUDE 'par12.f'
      REAL*8 X,LD0,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
C      WRITE(*,*) "INSIDE LD1"
      IF (XA .EQ. 0) THEN
         LD0=(1.0D0-RHO)*(1.0D0+GAMMA*RHO*(1+RHO**2))
      ELSE IF (IABS(XA) .eq. 1) Then
         LD0=0.50d0*(1.0d0-rho)*(rho+gamma*(2+2*(RHO**2)
     &+RHO**4))
      ELSE IF (IABS(XA) .eq. 2) Then
         LD0=0.50d0*(1.0D0-RHO)*RHO*(rho+gamma*(2+rho**2
     &+rho**4))
      ELSE IF (IABS(XA) .eq. 3) Then
         LD0=0.50d0*(1.0D0-RHO)*(RHO**3+gamma*(2+rho**2
     &+rho**4+rho**6))
      ELSE
         LD0=0.50d0*(1.0D0-RHO)*(rho**(IABS(XA)-3))
     &*(RHO**3+gamma*(1+rho**2+rho**4+rho**6))
      END IF
C      write(*,*) X,Y,W,XA,RHO, ROMIN,ROL,GAMMA,LD1
C      PAUSE
      RETURN
      END