* * bayes-lbm.f * (Bayes estimates for a micromovement model with * ^^^^^ * Linear Brownian Motion as Value Process * ^ ^ ^ * Copyright (C) 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 JSV-GBM model described C in the paper Yong Zeng (2005) "Bayes Estimation via Filtering c for a Simple Micro-movement Model of Asset Price with Discrete Noises" c Nonlinear Analysis Series A: Theory and Methods c (in press) *---------------------------------------------------------------- c "par-lbm.f" specifies the parameters of this program c date set: MSFT940102.dat (it has two columes: price, and c duration or waiting time between trades) c All files should be in the same directory c > f77 bayes-lbm.f (to generate "a.out") c > a.out (to run) c Warning: it may take a long time to finish the program. C---------------------------------------------------------------- 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 'par-lbm.f' INTEGER I,K,L,M,W REAL*8 PX0(N),PM(NM),PS(NS),PR(NR),P(N,NM,NS,NR),MU(NM), &RO(NR),X(N),PX(N) 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="16MU.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="16SG.LST",STATUS="UNKNOWN") DO 1000 M=1,NS PS(M)=0d0 c 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="16RO.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 COMPUTE THE POSTERIOR OF X----------------------------- OPEN(UNIT=10,FILE="16X.LST",STATUS="UNKNOWN") 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 IF ( PX(K) .GT. 0.01d0) THEN WRITE(10,'(G16.7E2,G16.12E1)') X(K),PX(K) END IF 1201 CONTINUE CLOSE(10) END SUBROUTINE DONDDIST3(P) EXTERNAL LD1,LD2,LD3,LD0 INCLUDE 'par-lbm.f' INTEGER I,K,L,M,NTI,J,W REAL*8 P(N,NM,NS,NR),A(NM,NS,N),B(NM,NS,N),X2,DT REAL*8 MU,RO,XB,XC,LD1,LD2,LD3,LD0 REAL*8 PI1(N),PI2(N) c asign A(NM,NS,N),B(NM,NS,N) accordingly------------------ DO 20 J=1,NM MU=MUMIN+(DBLE(J)-1.0000D0)*MUL DO 40 K=1,NS c SG=SGMIN+(DBLE(K)-1.0000D0)*SGL DO 60 L=1,N c X=XMIN-LL+(DBLE(L)-1.0000D0)/32.0D0 IF (MU .GE. 0D0) THEN A(J,K,L)=(SG(k)*SG(k)*32.0D0*32.0D0)/2.0d0+MU*32.0d0 B(J,K,L)=(SG(k)*SG(k)*32.0d0*32.0d0)/2.0d0 ELSE A(J,K,L)=(SG(k)*SG(k)*32.0D0*32.0D0)/2.0d0 B(J,K,L)=(SG(k)*SG(k)*32.0d0*32.0d0)/2.0d0-MU*32.0d0 end if 60 CONTINUE 40 CONTINUE 20 CONTINUE OPEN(UNIT=15, FILE='MUb.dat', STATUS="UNKNOWN") OPEN(UNIT=25, FILE='SGb.dat', STATUS="UNKNOWN") OPEN(UNIT=35, FILE='ROb.dat', STATUS="UNKNOWN") OPEN(UNIT=45, FILE='Xb.dat', STATUS="UNKNOWN") OPEN(UNIT=1,FILE=' MSFT940102.dat', STATUS='OLD') DO 6 I=1,NDATA READ(1,*) 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-) FROM t_{i-1} DO 100 M=1,NS DO 110 L=1,NM DO 120 W=1,NR c ----------------------------------------------------- c refine the partition based on the interarrival time c ----------------------------------------------------- NTI=IDNINT(DT/LT)+1 DO 125 J=1,NTI DO 50 K=1,N PI1(K)=P(K,L,M,W) 50 CONTINUE DO 130 K=1,N IF (K .EQ. 1) THEN P(K,L,M,W)=PI1(K)+(-(A(L,M,K)+B(L,M,K))*PI1(K) &+B(L,M,K+1)*PI1(K+1))*DT/DBLE(NTI) ELSE IF (K .EQ. N) THEN P(K,L,M,W)=PI1(K)+(A(L,M,K-1)*PI1(K-1)-(A(L,M,K) &+B(L,M,K))*PI1(K))*DT/DBLE(NTI) ELSE P(K,L,M,W)=PI1(K)+(A(L,M,K-1)*PI1(K-1)-(A(L,M,K) &+B(L,M,K))*PI1(K)+B(L,M,K+1)*PI1(K+1))*DT/DBLE(NTI) END IF 130 CONTINUE 125 CONTINUE 120 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 4: X2 IS ODD EIGHTH ELSE CALL DEND(LD3,P,X2) END IF CALL WOUTPUT(P,I,X2,DT) 6 continue close(1) CLOSE(15) CLOSE(25) CLOSE(35) CLOSE(45) RETURN END SUBROUTINE DEND(FUNC,P,X) include 'par-lbm.f' INTEGER K,L,M,W REAL*8 FUNC,P(N,NM,NS,NR),X,XA,TL C obtain a measure after jump and SUM UP THE TOTEL TL TL=0.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 P(K,L,M,W)=P(K,L,M,W)*FUNC(X,XA,W) TL=TL+P(K,L,M,W) 186 CONTINUE 184 CONTINUE 182 CONTINUE 180 CONTINUE IF (TL .EQ. 0D0) THEN PAUSE "TL=0" END IF c normalize the measure to get probability 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)=P(K,L,M,W)/TL 286 CONTINUE 284 CONTINUE 282 CONTINUE 280 CONTINUE RETURN END SUBROUTINE WOUTPUT(P,I,X2,DT) INCLUDE 'par-lbm.f' INTEGER I,K,L,M,W REAL*8 PM(NM),PS(NS),PR(NR),P(N,NM,NS,NR),MU(NM), &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 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 EMU=EMU+MU(L)*PM(L) EMU2=EMU2+MU(L)*MU(L)*PM(L) 901 CONTINUE SMU=DSQRT(EMU2-EMU*EMU) WRITE(15,'(G16.6E4,g16.6E4)') EMU,SMU 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 ESG=ESG+SG(M)*PS(M) ESG2=ESG2+SG(M)*SG(M)*PS(M) 1001 CONTINUE SSG=DSQRT(ESG2-ESG*ESG+1.0E-20) WRITE(25,'(G16.6E4,g16.6E4)') ESG,SSG 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 ERO=ERO+RO(W)*PR(W) ERO2=ERO2+RO(W)*RO(W)*PR(W) 1101 CONTINUE SRO=DSQRT(ERO2-ERO*ERO) WRITE(35,'(G16.6E4,g16.6E4)') ERO,SRO 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) 1201 CONTINUE SX=DSQRT(EX2-EX*EX) WRITE(45,'(G16.6E4,g16.6E4)') EX,SX RETURN END FUNCTION LD1(X,Y,W) INCLUDE 'par-lbm.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-lbm.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-lbm.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