# R --vsize=50M --nsize=1M
#it is modified from filter/lsigma/lsigma.s 
#This one is corrected from filter/lsigma/lsigma.s 
#	by replacing STT[S1T[i]] by STT[ST1[i]+1]
# simulate the tick data in three steps

# step 1: simulate the GMB, X_t  in Poisson random time
# dX_t/X_t = mu dt + sigma dW_t with sigma change over time according to
# ds_t = (S_{N(t)} - s_t) dN(t)
# X0, initial value
# mu: drift
# sgv: vector of sd of BM, for different periods
# ma: mean of inter-trade occurrance time
# ms: mean of sigma change time
# lg: length of X_t
# lm: length of mu
# ls: length of sigma
# 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<-80
mu<-4.5e-08
ma<-0.9
lg<-90000
ms<-3.75e-04
sgmin<-0.00004
sgmax<-0.0004
ls<-50
aph<-0.4
bt<-0.2
ro<-0.8 # that is c in the prg**.f
sink("sigma.lis")
set.seed(283)
#generate standard normal r.v. for diffusion
RN<-rnorm(lg,0,1)
print(mean(RN))
print(var(RN))
#generate the duration for  sigma
RST<-rexp(ls,ms)
cat("mean(RST)")
print(mean(RST))
cat("var(RST)")
print(var(RST))
#generate the duration of tradeing time
RT<-(rexp(lg,ma))
cat("mean(RT)")
print(mean(RT))
cat("var(RT)")
print(var(RT))
#generate Bin(1,0.5) for V
RB<-rbinom(lg,1,0.5)
print(mean(RB))
print(var(RB))
#generate geometric r.v. for V
RG<-rgeom(lg,ro)
print(mean(RG))
print(var(RG))
print(sum((-1+2*RB)*RG))
#generate uniform for Biasing
RU<-runif(lg,0,1)
print(mean(RU))
print(var(RU))
#compute the changing time of sigma
ST<-rep(ls,0)
ST[1]<-RST[1]
if (ls >1)
  {
for (i in 2:ls)
  {
    ST[i]<-ST[i-1]+RST[i]
  }
}
cat("ST:")
print(ST)
#compute the tradeing time
T1<-rep(0,lg)
T1[1]<-RT[1]
for (i in 2:lg)
	{
	T1[i]<-T1[i-1]+RT[i]
	}
#print(RT)
cat("T1[lg]")
print(T1[lg])
#adjust the change time of sigma
for (i in 1:ls)
  {
    if (ST[ls-i+1] > T1[lg])
      {
        ST<-ST[-ls+i-1]
      }
  }
ls<-length(ST)
#generate sigma vector
sg<-runif(ls+1, sgmin,sgmax)
cat("sg")
print(sg)
# determine the change point of sigma, ST1
ST<-append(1,ST)
ST1<-rep(0,ls+1)
ls1<-ls+1
for (j in 2:ls1)
  {
    k<-ST1[j-1]+1
      for (i in k:lg)
        {
          if (T1[i] <= ST[j] && ST[j]< T1[i+1])
            {
              ST1[j]<-i
             }
         }
   }
if (ST1[ls] < lg)
  {
    ST1<-append(ST1,lg)
  }
DST<-ST1[-1] - ST1[-(ls+2)]
cat("ST1")
print(ST1)
cat("DST")
print(DST)
#generate vector for sigma at trading times
#old approach
ST1<-ST1[-1]
STT<-rep(sg[1],ST1[1])
for (i in 2:ls1)
  {
    b<-rep(sg[i], ST1[i]-ST1[i-1])
    STT<-append(STT,b)
  }
#new approach
STTn<-rep(sg,DST)
#check both approach get the same answer
print(sum(STT-STTn))
#adjust STT at the change points
ST<-ST[-1]
for (i in 1:ls)
  {
    print(STT[ST1[i]+1])
    STT[ST1[i]+1]<-sqrt((sg[i]^2*(ST[i]-T1[ST1[i]])+sg[i+1]^2*(T1[ST1[i]+1]-ST[i]))/(T1[ST1[i]+1]-T1[ST1[i]]))
# corrected
    print(STT[ST1[i]+1])
  }

#generate the diffusion process with mu and sigma changeing over time
X<-rep(0,lg)
#X[1]<-X0+mu*RT[1]+STT[1]*sqrt(RT[1])*RN[1]
X[1]<-X0*exp((mu-0.5*STT[1]^2)*RT[1]+STT[1]*sqrt(RT[1])*RN[1])
for (i in 2:lg)
     {
#      X[i]<-X[i-1]+mu*RT[i]+STT[i]*sqrt(RT[i])*RN[i]
      X[i]<-X[i-1]*exp((mu-0.5*STT[i]^2)*RT[i]+STT[i]*sqrt(RT[i])*RN[i])
      }

print(X[1:10])
cat("X[lg]")
print(X[lg])
# 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
y1<-((-1+2*RB)*RG+round(8*X+0.00000000001))/8
print(RB[1:10])
print(RG[1:10])
# step 3: move the odd eighths to quarters and halves
#         with biased probabilities
print(RU[1:10])
y2<-rep(0,lg)
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]+0.00000000001))*8==1 || (y1[i]-floor(y1[i]+0.00000000001))*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]+0.00000000001))*8==1 || (y1[i]-floor(y1[i]+0.00000000001))*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)
#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()

sink()
#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,"gsgj.dat",row.names=FALSE)
write(X, "gsgj-x.dat",ncolumns=1)