* * bayes-jsvgbm.f * (Bayes estimates for a micromovement model with * ^^^^^ * Jumping Stochastic Volatility Geometric Brownian Motion as 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 JSV-GBM model described C in the paper Yong Zeng (2004) "Estimating Stochastic Volatility via c Filtering for the Micro-movement of Asset Prices" c IEEE Transactions on Automatic Control, 2004 c Vol. 49, pp. 338-348. *---------------------------------------------------------------- c "par-jsvgbm.f" specifies the parameters of this program c > f77 bayes-jsvgbm.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 (tick size) 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-jsvgbm.f' INTEGER I,K,L,M,W,LD REAL*8 PX(N),PM(NM),PS(NS),PR(NR),PJ(N,NM,NS,NL,NR),MU(NM), &RO(NR),X(N),PL(NL), LBD(NL) c ASSIGN INITIAL DISTRIBUTION OF X_0 --------------- DO 10 I=1,N PX(I)=0.0d0 10 CONTINUE PX(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 Lambda, PL(NL)------------------- DO 50 LD=1,NL PL(LD)=1.0D0/DBLE(NL) 50 CONTINUE c -------------------------------------------------- C ASSIGN PRIOR FOR RHO, PR(NR)---------------------- DO 40 W=1,NR PR(W)=1.0D0/DBLE(NR) 40 CONTINUE c -------------------------------------------------- c COMPUTE INITIAL DIST OF P(X_0,MU,SIGMA | Z_0) ---- DO 600 K=1,N DO 400 L=1,NM DO 300 LD=1,NL DO 200 M=1,NS DO 190 W=1,NR PJ(K,L,M,LD,W)=PX(K)*PM(L)*PS(M)*PL(LD)*PR(W) 190 CONTINUE 200 CONTINUE 300 CONTINUE 400 CONTINUE 600 CONTINUE C ------------------------------------------------ CALL DONDDIST3(PJ) 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 DO 960 LD=1,NL PM(L)=PM(L)+PJ(K,L,M,LD,W) 960 CONTINUE 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="16SG.LST",STATUS="UNKNOWN") DO 1000 M=1,NS PS(M)=0d0 SG=SGMIN+(dble(M)-1.00d0)*SGL DO 1030 K=1,N DO 1040 W=1,NR DO 1050 L=1,NM DO 1060 LD=1,NL PS(M)=PS(M)+PJ(K,L,M,LD,W) 1060 CONTINUE 1050 CONTINUE 1040 CONTINUE 1030 CONTINUE WRITE(20,'(G16.6E4,G16.6E4)') SG,PS(M) 1000 CONTINUE CLOSE(20) C ------------------------------------------------------ C COMPUTE THE POSTERIOR OF LAMBDA----------------------- OPEN(UNIT=25,FILE="16LD.LST",STATUS="UNKNOWN") DO 1500 LD=1,NL PL(LD)=0d0 LBD(LD)=LBDMIN+(dble(LD)-1.00d0)*LBDL DO 1530 K=1,N DO 1540 W=1,NR DO 1550 L=1,NM DO 1560 M=1,NS PL(LD)=PL(LD)+PJ(K,L,M,LD,W) 1560 CONTINUE 1550 CONTINUE 1540 CONTINUE 1530 CONTINUE WRITE(25,'(G16.6E4,G16.6E4)') LBD(LD),PL(LD) 1500 CONTINUE CLOSE(25) c --------------------------------------------------- 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 DO 1160 LD=1,NL PR(W)=PR(W)+PJ(K,L,M,LD,W) 1160 CONTINUE 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 DO 1261 LD=1,NL PX(K)=PX(K)+PJ(K,L,M,LD,W) 1261 CONTINUE 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(PJ) EXTERNAL LD1,LD2,LD3,LD0 INCLUDE 'par-jsvgbm.f' INTEGER I,K,L,M,NTI,J,W,LD REAL*8 PJ(N,NM,NS,NL,NR),A(NM,NS,N),B(NM,NS,N),X2,DT,DDT REAL*8 MU,RO,XB,XC,LD1,LD2,LD3,LD0,LBD REAL*8 PI1(N),PSM(N,NM,NL,NR) 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 SG=SGMIN+(DBLE(K)-1.0000D0)*SGL DO 60 L=1,N X=XMIN-LL+(DBLE(L)-1.0000D0)/64.0D0 IF (MU .GE. 0D0) THEN A(J,K,L)=(SG*SG*X*X*64.0D0*64.0D0)/2.0d0+MU*X*64.0d0 B(J,K,L)=(SG*SG*X*X*64.0d0*64.0d0)/2.0d0 ELSE A(J,K,L)=(SG*SG*X*X*64.0D0*64.0D0)/2.0d0 B(J,K,L)=(SG*SG*X*X*64.0d0*64.0d0)/2.0d0-MU*X*64.0d0 end if 60 CONTINUE 40 CONTINUE 20 CONTINUE OPEN(UNIT=32,FILE='midout.flt',STATUS="UNKNOWN") OPEN(UNIT=25, FILE='SG.dat', STATUS="UNKNOWN") OPEN(UNIT=1, FILE='smlt-jsvgbm.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 NTI=IDNINT(DT/LT)+1 DDT=DT/DBLE(NTI) c ----------------------------------------------------- c refine the partition based on the interarrival time c ----------------------------------------------------- DO 8 J=1,NTI c ------------------------------------------------------ c Define the mean of prob of sigma (PSM) for K,L,LD,W of X,MU, c LD,RO at each step c ------------------------------------------------------ DO 10 K=1,N DO 30 L=1,NM DO 50 LD=1,NL DO 70 W=1,NR PSM(K,L,LD,W)=0D0 DO 90 M=1,NS PSM(K,L,LD,W)=PSM(K,L,LD,W)+PJ(K,L,M,LD,W) 90 CONTINUE PSM(K,L,LD,W)=PSM(K,L,LD,W)/DBLE(NS) 70 CONTINUE 50 CONTINUE 30 CONTINUE 10 CONTINUE c loop 100: approXimate p(K,L,M,W,t_i-) FROM t_{i-1} to each step DO 100 M=1,NS DO 110 L=1,NM DO 120 W=1,NR DO 122 LD=1,NL LBD=LBDMIN+(LD-1.0d0)*LBDL DO 125 K=1,N PI1(K)=PJ(K,L,M,LD,W) 125 CONTINUE DO 130 K=1,N IF (K .EQ. 1) THEN PJ(K,L,M,LD,W)=PI1(K)+(-(A(L,M,K)+B(L,M,K))*PI1(K) &+B(L,M,K+1)*PI1(K+1)+LBD*(PSM(K,L,LD,W)-PI1(K)))*DDT if (PJ(K,L,M,LD,W) .lt. 0) then write(*,*) "l,m,w,j,k,nti,dt", l,m,w,j,k,nti,dt write(*,*) "PJ(K,L,M,LD,W)",PJ(K,L,M,LD,W) write(*,*) "pi1(k)",pi1(k),"pi1(k+1)",pi1(k+1) write(*,*) "A(L,M,K),B(L,M,K),B(L,M,K+1)",A(L,M,K), &B(L,M,K),B(L,M,K+1) pause "PJ(K,L,M,LD,W) .lt. 0" end if ELSE IF (K .EQ. N) THEN PJ(K,L,M,LD,W)=PI1(K)+(A(L,M,K-1)*PI1(K-1)-(A(L,M,K) &+B(L,M,K))*PI1(K)+LBD*(PSM(K,L,LD,W)-PI1(K)))*DDT if (PJ(K,L,M,LD,W) .lt. 0) then write(*,*) "l,m,w,j,k,nti,dt", l,m,w,j,k,nti,dt write(*,*) "PJ(K,L,M,LD,W)",PJ(K,L,M,LD,W),"pi1(k)",pi1(k) write(*,*) "pi1(k-1)", pi1(k-1) write(*,*) "A(L,M,K-1),A(L,M,K),B(L,M,K),", &A(L,M,K-1),A(L,M,K),B(L,M,K) pause "PJ(K,L,M,LD,W) .lt. 0" end if ELSE PJ(K,L,M,LD,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)+LBD*(PSM(K,L,LD,W) &-PI1(K)))*DDT if (PJ(K,L,M,LD,W) .lt. 0) then write(*,*) "l,m,w,j,k,nti,dt", l,m,w,j,k,nti,dt write(*,*) "PJ(K,L,M,LD,W)",PJ(K,L,M,LD,W),"pi1(k)",pi1(k) write(*,*) "pi1(k-1),pi1(k+1)", pi1(k-1),pi1(k+1) write(*,*) "A(L,M,K-1),A(L,M,K),B(L,M,K),B(L,M,K+1)", &A(L,M,K-1),A(L,M,K),B(L,M,K),B(L,M,K+1) pause "PJ(K,L,M,LD,W) .lt. 0" end if END IF 130 CONTINUE 122 CONTINUE 120 CONTINUE 110 CONTINUE 100 CONTINUE 8 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,PJ,X2,I) C CASE 2: X2 IS ODD QUARTERS ELSE IF (DMOD(XC,1d0) .EQ. 0) THEN CALL DEND(LD2,PJ,X2,I) C CASE 3: X2 IS ODD EIGHTH ELSE CALL DEND(LD3,PJ,X2,I) END IF c For the points where the volatility is about to changes, c middle-step output is written. if ((i .eq. 1) .or. (i .eq. 2166) .or. (i .eq. 2676) & .or. (i .eq. 7790) .or. (i .eq. 8113) .or. (i .eq. 8870) & .or. (i .eq. 11793) .or. (i .eq. 11972) .or. (i .eq. 13517) & .or. (i .eq. 15267) .or. (i .eq. 17161) .or. (i .eq. 19235) & .or. (i .eq. 20012) .or. (i .eq. 20299) .or. (i .eq. 22355) & .or. (i .eq. 24964) .or. (i .eq. 29856) .or. (i .eq. 31732) & .or. (i .eq. 31830) .or. (i .eq. 32053) .or. (i .eq. 34738) & .or. (i .eq. 35231) .or. (i .eq. 35767) .or. (i .eq. 38157) & .or. (i .eq. 39158) .or. (i .eq. 39607) .or. (i .eq. 41360) & .or. (i .eq. 43226) .or. (i .eq. 44478) .or. (i .eq. 51559) & .or. (i .eq. 53894) .or. (i .eq. 57307) .or. (i .eq. 58145) & .or. (i .eq. 61812) .or. (i .eq. 62334) .or. (i .eq. 62627) & .or. (i .eq. 64298) .or. (i .eq. 64900) .or. (i .eq. 71690) & .or. (i .eq. 73699) .or. (i .eq. 74971) .or. (i .eq. 75661) & .or. (i .eq. 76336) .or. (i .eq. 81416) .or. (i .eq. 84814) & .or. (i .eq. 90000)) then CALL WOUTPUT(PJ,I,X2,DT) end if CALL WOUTPUTsg(PJ) 6 continue close(1) close(25) CLOSE(32) RETURN END SUBROUTINE DEND(FUNC,PJ,X,I) include 'par-jsvgbm.f' INTEGER K,L,M,W,LD REAL*8 FUNC,PJ(N,NM,NS,NL,NR),X,XA,TL C obtain a measure after a trade 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)*64.0D0)/64.0D0 DO 188 LD=1,NL PJ(K,L,M,LD,W)=PJ(K,L,M,LD,W)*FUNC(X,XA,W) TL=TL+PJ(K,L,M,LD,W) 188 CONTINUE 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 XA=(DBLE(K)-1D0+(XMIN-LL)*64.0D0)/64.0D0 DO 288 LD=1,NL PJ(K,L,M,LD,W)=PJ(K,L,M,LD,W)/TL 288 CONTINUE 286 CONTINUE 284 CONTINUE 282 CONTINUE 280 CONTINUE RETURN END SUBROUTINE WOUTPUT(PJ,I,X2,DT) INCLUDE 'par-jsvgbm.f' INTEGER I,K,L,M,W REAL*8 PM,PS,PR,PL,PJ(N,NM,NS,NL,NR),MU,SG, &RO,X2,DT,EMU,EMU2,SMU,ESG,ESG2,SSG,ERO,ERO2,SRO,EX,EX2, &SX,PX,X,LBD,ELBD, ELBD2, SLBD EMU=0D0 EMU2=0D0 ESG=0D0 ESG2=0D0 ELBD=0D0 ELBD2=0D0 ERO=0D0 ERO2=0D0 EX=0D0 EX2=0D0 WRITE(32,*) "I,X2,DT",I,X2,DT WRITE(32,*) "posterior of MU" DO 901 L=1,NM PM=0d0 MU=MUMIN+(dble(L)-1.000d0)*MUL DO 931 K=1,N DO 941 W=1,NR DO 951 M=1,NS DO 961 LD=1,NL PM=PM+PJ(K,L,M,LD,W) 961 CONTINUE 951 CONTINUE 941 CONTINUE 931 CONTINUE WRITE(32,'(g16.6E4,G16.6E4)') MU,PM EMU=EMU+MU*PM EMU2=EMU2+MU*MU*PM 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=0d0 SG=SGMIN+(dble(M)-1.00d0)*SGL DO 1031 K=1,N DO 1041 W=1,NR DO 1051 L=1,NM DO 1061 LD=1,NL PS=PS+PJ(K,L,M,LD,W) 1061 CONTINUE 1051 CONTINUE 1041 CONTINUE 1031 CONTINUE if (PS .gt. .05d0) then WRITE(32,'(G16.6E4,G16.6E4)') SG,PS end if ESG=ESG+SG*PS ESG2=ESG2+SG*SG*PS 1001 CONTINUE SSG=DSQRT(ESG2-ESG*ESG) WRITE(32,*) "BAYESIAN ESTIMATE OF SG AND ITS SD",ESG,SSG WRITE(32,*) "POSTERIOR OF LBD" DO 1401 LD=1,NL PL=0d0 LBD=LBDMIN+(dble(LD)-1.00d0)*LBDL DO 1431 K=1,N DO 1441 W=1,NR DO 1451 L=1,NM DO 1461 M=1,NS PL=PL+PJ(K,L,M,LD,W) 1461 CONTINUE 1451 CONTINUE 1441 CONTINUE 1431 CONTINUE WRITE(32,'(G16.6E4,G16.6E4)') LBD,PL ELBD=ELBD+LBD*PL ELBD2=ELBD2+LBD*LBD*PL 1401 CONTINUE SLBD=DSQRT(ELBD2-ELBD*ELBD+1.0d-12) WRITE(32,*) "BAYESIAN ESTIMATE OF LBD AND ITS SD",ELBD,SLBD WRITE(32,*) "POSTEROR OF RHO" DO 1101 W=1,NR PR=0d0 RO=ROMIN+(dble(W)-1.00d0)*ROL DO 1131 K=1,N DO 1141 M=1,NS DO 1151 L=1,NM DO 1161 LD=1,NL PR=PR+PJ(K,L,M,LD,W) 1161 CONTINUE 1151 CONTINUE 1141 CONTINUE 1131 CONTINUE WRITE(32,'(G16.6E4,G16.6E4)') RO,PR ERO=ERO+RO*PR ERO2=ERO2+RO*RO*PR 1101 CONTINUE SRO=DSQRT(ERO2-ERO*ERO+1d-12) WRITE(32,*) "BAYESIAN ESTIMATE OF RO AND ITS SD",ERO,SRO WRITE(32,*) "POSTERIOR OF X" DO 1201 K=1,N PX=0d0 X=XMIN-LL+(dble(K)-1.00d0)/64.0d0 DO 1231 W=1,NR DO 1241 M=1,NS DO 1251 L=1,NM DO 1261 LD=1,NL PX=PX+PJ(K,L,M,LD,W) 1261 CONTINUE 1251 CONTINUE 1241 CONTINUE 1231 CONTINUE EX=EX+X*PX EX2=EX2+X*X*PX IF ( PX .GT. 0.05d0) THEN WRITE(32,'(G16.6E4,G16.6E4)') X, PX END IF 1201 CONTINUE SX=DSQRT(EX2-EX*EX) WRITE(32,*) "BAYESIAN ESTIMATE OF X AND ITS SD",EX,SX RETURN END SUBROUTINE WOUTPUTsg(PJ) INCLUDE 'par-jsvgbm.f' INTEGER I,K,L,M,W REAL*8 PS,PJ(N,NM,NS,NL,NR),SG,ESG,ESG2,SSG ESG=0D0 ESG2=0D0 DO 10010 M=1,NS PS=0d0 SG=SGMIN+(dble(M)-1.00d0)*SGL DO 10310 K=1,N DO 10410 W=1,NR DO 10510 L=1,NM DO 10610 LD=1,NL PS=PS+PJ(K,L,M,LD,W) 10610 CONTINUE 10510 CONTINUE 10410 CONTINUE 10310 CONTINUE ESG=ESG+SG*PS ESG2=ESG2+SG*SG*PS 10010 CONTINUE SSG=DSQRT(ESG2-ESG*ESG) write(25,*) ESG, SSG RETURN END FUNCTION LD1(X,Y,W) INCLUDE 'par-jsvgbm.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-jsvgbm.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-jsvgbm.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