rm(list=ls(all=TRUE)) 
cat("\014") 

#########################
EMHMMExp<-function(data, k) {
T=length(data)
stat=matrix(rep(data,k),nrow=k,byrow=TRUE) #k*T
#hidden state
a=matrix(rexp(T*k,1),nrow=k) #k*T
a=a/matrix(rep(colSums(a),k),nrow=k,byrow=TRUE)
#Parameters
c=rowMeans(a) #k*1
nu=rowSums(a*stat)/rowSums(a) #k*1
P=matrix(rexp(k*k,1),nrow=k) #k*k, transition matrix,
P=P/matrix(rep(rowSums(P),k),nrow=k) ##Pij=Prob(q_{t+1}=j|q_t=i)
#Set up
threshold=10^(-5)
continue=TRUE
while(continue){
density=exp(-(1/nu)%*%t(data))/matrix(rep(nu,T),nrow=k) #k*T

#E step alpha-beta algorithm
scaledAlpha=as.matrix(1/k*density[,1]) #alpha(q0) k*1 
scale=sum(scaledAlpha)
scaledAlpha=scaledAlpha/scale
for (t in 2:T) {
scaledAlpha=cbind(scaledAlpha,t(density[,t]*(t(scaledAlpha[,t-1])%*%P)))
scale=c(scale,sum(scaledAlpha[,t]))
scaledAlpha[,t]=scaledAlpha[,t]/scale[t]
}
#scaledAlpha[,t]=1/scale[t]*alpha[,t]
scaledBeta=cbind(P%*%density[,T-1]/scale[T-1],density[,T]/scale[T]) 
#scaledBeta[,t]=1/scale[t]beta[,t]
for (t in (T-2):1) {
scaledBeta=cbind(P%*%(density[,t+1]*scaledBeta[,1])/scale[t],
               scaledBeta)
}
aNew=scaledAlpha*scaledBeta
aNew=aNew/matrix(rep(colSums(aNew),k),nrow=k,byrow=TRUE)

#M step
cNew=rowMeans(aNew)
nuNew=rowSums(aNew*stat)/rowSums(aNew)
PNew=mapply(function(x) rowSums(matrix(rep(aNew[x,1:(T-1)],k),nrow=k,byrow=TRUE)*aNew[,2:T]),
        1:k)
PNew=t(PNew)
PNew=PNew/matrix(rep(rowSums(PNew),k),nrow=k)
#Stop criterion
if(max(c(abs(cNew-c),abs(nuNew-nu),abs(c(P-PNew))))<threshold) {
continue=FALSE
}
a=aNew
c=cNew
nu=nuNew
P=PNew
}
logL=sum(log(c%*%density))
return (list(c=c,lambda=1/nu,P=P,logL=logL,stateProb=a))
}

#############
fullData=read.csv("D:/Dropbox/Dropbox/Limit Order Book/market making/A.csv", nrow=1000,
                  header=TRUE)
fullData=read.csv("/Users/jingguo/Dropbox/Limit Order Book/market making/A.csv", nrow=1000,
                  header=TRUE)
end=nrow(fullData)
#time=fullData$time
#date=fullData$date
time=as.POSIXct(paste(fullData$date,fullData$time))
timeDiff=as.numeric(time[2:end]-time[1:(end-1)])
#######################


if (FALSE) {
k=10
HMMExp=EMHMMExp(timeDiff,k)
####
logL=-10^5
rep=10
for (i in 1:rep) {
  temp=EMHMMExp(timeDiff,k)
  print(temp)
  if (logL<temp$logL) {
    logL=temp$logL
    HMMExp=temp
  }
}

####

qqplot(unlist(mapply(rexp,round(end*HMMExp$c),HMMExp$lambda)),
       timeDiff,xlab=paste(k,"-Mixture Exponential"),
       main="Marginal Distribution Fitting by HMM")
#qqline(timeDiff,col='blue')
lines(0:max(timeDiff),0:max(timeDiff),col='blue')
}

###################
upperLowerMedian<-function(alpha,lambda,n) {
  rVec=alpha%*%t(mapply(function(x) rexp(n,x), lambda))
  return(c(quantile(rVec,c(.05,.5,.95))))
}
upperLowerMedianProcess<-function(alpha,lambda,n) {
  #alpha k*T
  #k=length(alpha[,1])
  T=length(alpha[1,])
  return(mapply(function(t) upperLowerMedian(alpha[,t],lambda,n), 1:T))
}

forwardFilter<-function(data,pi0,lambda,P,n) {
  T=length(data)
  k=length(pi0)
  density=exp(-(lambda)%*%t(data))*matrix(rep(lambda,T),nrow=k) #k*T
  
  #E step alpha-beta algorithm
  scaledAlpha=as.matrix(t(density[,1]*(t(pi0)%*%P))) #alpha(q0) k*1 
  scale=sum(scaledAlpha)
  scaledAlpha=scaledAlpha/scale
  for (t in 2:(T-1)) {
    scaledAlpha=cbind(scaledAlpha,t(density[,t]*(t(scaledAlpha[,t-1])%*%P)))
    scale=c(scale,sum(scaledAlpha[,t]))
    scaledAlpha[,t]=scaledAlpha[,t]/scale[t]
  }
  return(upperLowerMedianProcess(scaledAlpha,lambda,n))
}


k=10
fitProportion=.95
fitData=timeDiff[1:(end*fitProportion)]
testData=timeDiff[(end*fitProportion+1):(end-1)]
fit=EMHMMExp(fitData,k)
filter=forwardFilter(testData,fit$stateProb[,(end*fitProportion)],fit$lambda,
                     fit$P,1000)
#plot(time[(end*fitProportion+1):(end-1)],testData,xlab='time')
#lines(time[(end*fitProportion+2):(end-1)],filter[1,],col='blue')
#lines(time[(end*fitProportion+2):(end-1)],filter[2,],col='red')
#lines(time[(end*fitProportion+2):(end-1)],filter[3,],col='blue')

plot((floor(end*fitProportion)+2):(end-1),
     timeDiff[(floor(end*fitProportion)+2):(end-1)], 
     xlab=paste("Out-of-Sample Arrival Time, Index",
                (floor(end*fitProportion)+1),":",
                (end-1)),
     ylab = paste("Arrival Time",(floor(end*fitProportion)+1),":",
                  (end-1)),
     main=paste("Real Versus Predicted Arrival Time,",
                "In-Sample 1:",floor(end*fitProportion)))
legend('topleft',legend = c("Real Data","95% Upper Bound",
                            "Median","5% Lower Bound"),
       col = c('black','blue','red','blue'),
       lty = c(1,2,1,2) )
lines((floor(end*fitProportion)+2):(end-1), filter[1,],col='blue',lty=2)
lines((floor(end*fitProportion)+2):(end-1), filter[2,],col='red')
lines((floor(end*fitProportion)+2):(end-1), filter[3,],col='blue',lty=2)