# 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<-24
a<-0.00000005
b<-0.00004
la<-0.0003
lg<-300
aph<-0.4
bt<-0.2
ro<-0.9 # that is c in the prg**.f

X<-rep(lg,0)
set.seed(211)
RT<-(rexp(lg,la))
print(mean(RT))
print(var(RT))
RN<-rnorm(lg,0,1)
print(mean(RN))
print(var(RN))
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))/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 (RU[i] < aph || RU[i]==aph)
	   {
	      if ((y1[i]-floor(y1[i]))*8==1 || (y1[i]-floor(y1[i]))*8==5)
	      {
	      y2[i]<-y1[i]+1/8
	      }
	      else
	      {
	      y2[i]<-y1[i]-1/8
	      }
	    }
	    else
	    {
              if (RU[i] < aph+bt || RU[i]==aph+bt)
	      {
	        if ((y1[i]-floor(y1[i]))*8==1 || (y1[i]-floor(y1[i]))*8==5)
	        {
	        y2[i]<-y1[i]-1/8
	        }
	        else
	        {
	        y2[i]<-y1[i]+1/8
	        }
	      }
	      else
	      {
	      y2[i]<-y1[i]
	      }
	    }
	  }
	}
print(y2[1:10])
T2<-T1[1:1000]
X2<-X[1:1000]
y1a2<-y1a[1:1000]
y22<-y2[1:1000]
postscript("tem1.ps", he=6,wi=9.5)
#eighths<-y2-floor(y2)
plot(T1,X,type="l")
plot(T1,y1a,type="l")
plot(T1,y2,type="l")
plot(T2,X2,type="l")
plot(T2,y1a2,type="l")
plot(T2,y22,type="l")
#title(sub="c=0.8")
#hist(eighths) 
graphics.off()

postscript("tem1a.ps", he=11,wi=8)
par(mfrow=c(3,2))
hist(RT)
hist(RN)
hist(RB)
hist(RG)
hist(RU)
graphics.off()
yt0<-cbind(y2,RT)
write.table(yt0,"temaa.dat")
#rm(list=ls())