# simulate the tick data in three steps

# step 1: simulate the GMB, X_t  in Poisson random time
# dX_t/X_t = a dt + b dW_t
# X0, initial value
# a: drift
# b: sd of BM.
# la: intensity of Poisson random arrival time
# lg: length of X_t
# aph: the probability of moving odd eighths to the nearest
	# quarters but not halves
# bt:	the probability of moving odd eighths to the nearest
	# halves or integers
#options(object.size=5e8)
X0<-83
a<-4.40e-08
b<-1.2e-04
la<-0.06
lg<-32500
aph<-0.225
bt<-0.066
gm<-0.300
ro<-0.8 # that is c in the prg**.f
#sink('s9403.slt')
X<-rep(lg,0)
set.seed(783)
RG<-rgeom(lg,ro)
print(mean(RG)) ## mean is (1-ro)/ro
print(var(RG))
RT<-(rexp(lg,la))
print(mean(RT))
print(var(RT))
RN<-rnorm(lg,0,1)
print(mean(RN))
print(var(RN))
RB<-rbinom(lg,1,0.5)
print(mean(RB))
print(var(RB))
print(sum((-1+2*RB)*RG))
T1<-rep(lg,0)
T1[1]<-RT[1]
for (i in 2:lg)
	{
	T1[i]<-T1[i-1]+RT[i]
	}
#print(RT)
X[1]<-X0*exp((a-0.5*b^2)*RT[1]+b*sqrt(RT[1])*RN[1])
for (i in 2:lg)
	{
	X[i]<-X[i-1]*exp((a-0.5*b^2)*RT[i]+b*sqrt(RT[i])*RN[i])
	}
print(X[1:10])

# step 2': round off to the greatest eighth smaller than its value
y1a<-round(8*X)/8  # floor() or ceiling() sometimes
print(y1a[1:10])
#step 2: pick the y1
RB<-rbinom(lg,1,0.5)
RG<-rgeom(lg,ro)
y1<-((-1+2*RB)*RG+round(8*X+0.00000000001))/8
print(mean(RB))
print(var(RB))
print(mean(RG))
print(var(RG))
print(sum((-1+2*RB)*RG))
print(RB[1:10])
print(RG[1:10])
# step 3: move the odd eighths to quarters and halves
#         with biased probabilities
RU<-runif(lg,0,1)
print(mean(RU))
print(var(RU))
print(RU[1:10])
y2<-rep(lg,0)
for (i in 1:lg)
        {
	if (y1[i]-floor(4*y1[i]+0.00000000001)/4==0)
	{
	 y2[i]<-y1[i]
	}
	else
	{if ((y1[i]-floor(y1[i]+0.00000000001))*8==1)
	   {if (RU[i] < aph)
        	   { y2[i]<-y1[i]+1/8}
                 else
                   {if (RU[i] < aph+bt)
                       {y2[i]<-y1[i]+3/8}
                     else
                       {if (RU[i] < aph+bt+gm)
			{ y2[i]<-y1[i]-1/8}
			else
			{
                         y2[i]<-y1[i]
                       }
                   }
               }
	    }
	 else
	 {if ((y1[i]-floor(y1[i]+0.00000000001))*8==3)
	   {if (RU[i] < aph)
        	   { y2[i]<-y1[i]-1/8}
                 else
                   {if (RU[i] < aph+bt)
                       {y2[i]<-y1[i]+1/8}
                     else
                       {if (RU[i] < aph+bt+gm)
			{ y2[i]<-y1[i]-3/8}
			else
			{
                         y2[i]<-y1[i]
                       }
                   }
               }
	    }
	  else
	  {if ((y1[i]-floor(y1[i]+0.00000000001))*8==5)
	   {if (RU[i] < aph)
        	   { y2[i]<-y1[i]+1/8}
                 else
                   {if (RU[i] < aph+bt)
                       {y2[i]<-y1[i]-1/8}
                     else
                       {if (RU[i] < aph+bt+gm)
			{ y2[i]<-y1[i]+3/8}
			else
			{
                         y2[i]<-y1[i]
                       }
                   }
               }
	    }
           else
	   {if (RU[i] < aph)
        	   { y2[i]<-y1[i]-1/8}
                 else
                   {if (RU[i] < aph+bt)
                       {y2[i]<-y1[i]-3/8}
                     else
                       {if (RU[i] < aph+bt+gm)
			{ y2[i]<-y1[i]+1/8}
			else
			{
                         y2[i]<-y1[i]
                       }
                   }
               }
	    }
	  }
	}
	}}
print(y2[1:10])
#sink()
T2<-T1[1:1000]
X2<-X[1:1000]
y1a2<-y1a[1:1000]
y12<-y1[1:1000]
y22<-y2[1:1000]
postscript("s9403.ps", he=6,wi=10)
par(las=1,mar=c(3.1,3.2,3,1.3),mgp=c(2.1,0.5,0))
#eighths<-y2-floor(y2)
day<-T1/23400
day2<-day[1:1000]
plot(day,X,type="l",ylab="price")
title("GBM")
plot(day,y1a,type="l",ylab="price")
title("Rounding Model: R(X, 1/8)")
plot(day,y1,type="l",ylab="price")
title("Rounding + Non-clustering: R(X,1/8) + V")
plot(day,y2,type="l",ylab="price")
title("Our Model: b[R(X,1/8) + V]")
#plot(day2,X2,type="l")
#plot(day2,y1a2,type="l")
#plot(day2,y12,type="l")
#plot(day2,y22,type="l")
#title(sub="c=0.8")
#hist(eighths) 
graphics.off()


yt0<-cbind(y2,RT)

write.table(yt0,"m9403.dat",row.names=FALSE)
#rm(T2,X2,y1a2,y22,y2,y1a,X0,a,b,la,lg,aph,bt,ro,X,y1)
#rm(RT,T1,RN,RB,RG,yt0,RU)