#############READ-IN DATA##################
#DATA PATH
#setwd("C:/Users/markus/Documents/Forschung/co-jumps/volajumps/Anwendung/DATA")

#LOBSTER DATA READ-IN; SELECT ASSET AND DATE
tab<-read.csv("GE_2009-03-18_34200000_57600000_message_1.csv",header=FALSE)
#tab<-read.csv("MMM_2009-03-18_34200000_57600000_message_1.csv",header=FALSE)

#PREPARE DATA TRADED PRICES
t<-(tab[,1]-34200)/1000
t<-t/max(t)
trade<-tab[,5][sort(c(which(tab[,2]==4),which(tab[,2]==5)))]/10000
tt<-t[sort(c(which(tab[,2]==4),which(tab[,2]==5)))]
trade<-log(trade)

#SELECT NO OF BLOCKS
K<-50

########SET DGP#########################
Y<-trade
time<-tt
Deltay<-diff(Y)
n<-length(Y)-1

#SELECT SPECTRAL CUTOFF
J<-40

w<-array(0,dim=c(J,K));w2<-array(0,dim=c(J,K));Spec<-array(0,dim=c(J,K));Spec2<-array(0,dim=c(J,K));Spec3<-array(0,dim=c(J,K))

#############SPECTRAL BASIS####################
phi<-function(j,k,x){
(sqrt(2/K)/(pi*j))*(sin(j*pi*(x-k/K)*K))*((k-1)/K<=x)*(x<=k/K)}
a<-1/2
phi2<-function(j,k,x){
(sqrt(2/K)/(pi*j))*(sin(j*pi*(x-(k-1-a)/K)*K))*((k-1-a)/K<=x)*(x<=(k-a)/K)}
################NOISE LEVEL ESTIMATION#############################
nlevel<-(2*length(diff(Y)[which(diff(Y)!=0)]))^(-1)*sum(diff(Y)^2)
FN<-function(x){
length(time[which(time<=x)])/length(time)}
level<-numeric(K)
for(k in 1:K){
level[k]<-nlevel/(K*(FN(k/K)-FN((k-1)/K)))}

#ODD FREQUENCIES
F[1]<-1
for(i in 1:J){
F[i+1]<-F[i]+2}


##########PILOT SPOT ESTIMATORS####################################
L<-5 #smoothing window (no of blocks) for pilot estimation
JPILOT<-15 #Spectral cutoff for pilot estimator
Volpil<-numeric(K)   #Vector of pilot estimates of spot volatility, local averages of PIVblock	
PIVblock<-numeric(K) #Vector of block-wise pilot spectral estimates with equal weights

#SPECTRAL STATISTICS
for(k in 1:K){
for(j in 1:J){
       Spec[j,k]<-sum(phi(j,k,0.5*(time[2:(n+1)]+time[1:((n+1)-1)]))*Deltay)
       Spec2[j,k]<-sum(phi(F[j],k,0.5*(time[2:(length(time))]+time[1:((length(time))-1)]))*Deltay)
       Spec3[j,k]<-sum(phi2(F[j],k,0.5*(time[2:(length(time))]+time[1:((length(time))-1)]))*Deltay)}
PIVblock[k]<-sum((1/JPILOT)*K^2*(1:JPILOT)^2*pi^2*((Spec[1:JPILOT,k])^2-level[k]/n))
}
#truncation step to eliminate price-jumps
TRP<-mean(PIVblock)*2*log(K) #PILOT CONSTANT THRESHOLD
if(max(PIVblock)>TRP){
PIVblock[which(PIVblock>TRP)]<-TRP
}
for(k in (floor(L/2)+1):(K-floor(L/2))){
Volpil[k]<-mean(PIVblock[(k-floor(L/2)):(k+floor(L/2))])}
for(k in 1:(floor(L/2))){
Volpil[k]<-mean(PIVblock[1:(k+floor(L/2))])}
for(k in ((K-floor(L/2)+1):K)){
Volpil[k]<-mean(PIVblock[(k-floor(L/2)):K])}

#SELECT LOCALLY ADAPTIVE TRUNCATION
trsh<-Volpil*2
#################SPECTRAL ESTIMATORS###########################
IVblock<-numeric(K);IVblocktrsh<-numeric(K);IVblock2<-numeric(K);IVblock3<-numeric(K);IVblocktest<-numeric(K)

######################ADAPTIVE EFFICIENT WEIGHTS##################
for(k in 1:K){for(j in 1:J){
w[j,k]<-(Volpil[k]+(level[k]/n)*pi^2*j^2*K^2)^{-2}}}
for(k in 1:K){
w[,k]<-w[,k]/sum(w[,k])}
for(k in 1:K){for(j in 1:J){
w2[j,k]<-(Volpil[k]+(level[k]/n)*pi^2*F[j]^2*K^2)^{-2}}
w2[,k]<-w2[,k]/sum(w2[,k])}
############LOCAL ESTIMATES##################
for(k in 1:K){
IVblock[k]<-max(sum(w[1:J,k]*K^2*(1:J)^2*pi^2*((Spec[1:J,k])^2-level[k]/n)),0)
IVblock2[k]<-sum(w2[1:J,k]*K^2*(F[1:J])^2*pi^2*Spec2[1:J,k]*Spec2[1:J,k])
IVblock3[k]<-sum(w2[1:J,k]*K^2*(F[1:J])^2*pi^2*Spec3[1:J,k]*Spec3[1:J,k])
IVblocktest[k]<-(IVblock[k])*(abs(IVblock[k])>trsh[k])+(IVblock[k])*(abs(IVblock[k])<=trsh[k])*sample(c(-0.1,0.1),1,prob=c(0.5,0.5))}
#TRUNCATION TO EXCLUDE JUMPS
IVblocktrsh[which(IVblock<=trsh)]<-IVblock[which(IVblock<=trsh)]
if(min(IVblocktrsh)==0){
IVblocktrsh[which(IVblocktrsh==0)]<-mean(IVblock)}

#############################TRUNCATION STEP####################


###ESTIMATORS#############
SPEV<-mean(IVblock) #original spectral estimator
trSPEV<-mean(IVblock[which(IVblock<=trsh)]) #truncated integrated covolatility estimator
SPEVadj<-0.5*mean(IVblock2)+0.5*mean(IVblock3)
SPECJ<-mean(IVblock[which(abs(IVblock)>trsh)])*1/K
#########################SPOT VOLA ESTIMATES######################
#SELECT WINDOW LENGTH SPOT ESTIMATORS
R<-5
volaleft<-numeric(K)
volaright<-numeric(K)
for(k in 1:(K-R)){
volaright[k]<-mean(IVblocktrsh[(k+1):(k+R)])}
for(k in (K-R+1):K){
volaright[k]<-mean(IVblocktrsh[(k+1):K])}
volaright[K]<-IVblocktrsh[K]
for(k in (R+1):K){
volaleft[k]<-mean(IVblocktrsh[(k-R):(k-1)])}
volaleft[1]<-IVblocktrsh[1]
for(k in 2:R){
volaleft[k]<-mean(IVblocktrsh[1:(k-1)])}
Vol<-0.5*(volaleft+volaright)


#estimate jumps block-wise
JUMPS<-numeric(K)
for(k in 1:K){
if(IVblock[k]>trsh[k]){
JUMPS[k]<-IVblock[k]}}
JUMPS<-sqrt(JUMPS)*sign(Spec[1,])
#estimate volatility motion
VOLAMOTION<-volaright-volaleft
#TEST FOR JUMPS
test<-(n/2)^(1/4)*sum(IVblocktest)*1/K
#STANDARDIZED
testsd<-(0.01*mean(IVblock^2))^(-1/2)*test
#TEST FOR PRICE VOLA CO-JUMPS
test2<-(R)*(nlevel)^(-1/2)*sum(2*sqrt((volaright[which(IVblock[(R+1):(K-R)]>trsh[(R+1):(K-R)])+R]+volaleft[which(IVblock[(R+1):(K-R)]>trsh[(R+1):(K-R)])+R])/2)-sqrt(volaright[which(IVblock[(R+1):(K-R)]>trsh[(R+1):(K-R)])+R])-sqrt(volaleft[which(IVblock[(R+1):(K-R)]>trsh[(R+1):(K-R)])+R]))
#P-VALUES
print("p-value test for jump test");(1-pnorm(testsd))
print("p-value test for price-vola cojumps");(1-pchisq(test2,df=length(which(IVblock[(R+1):(K-R)]>trsh[(R+1):(K-R)]))))


print("estimated variation");SPEV
print("estimated integrated volatility");trSPEV
print("estimated jump variation");SPECJ
print("block-wise jump estimates");JUMPS
print("spot volatility estimate");plot(Vol)