setwd()
###################
#US bond data: 13:00 (46800) to 15:30 (55800) -- press release is at 13:00
#Thomson Reuters data: 30 seconds, 2.5hours (1h before event 1.5h-2h after event), 301 observations per bond and event
# Jump test conducted at observation 120 (50370 sec after midnight which is 12:59:30pm) to 131 (50700 sec after midnight is 1:05pm)

library(yuima)

n<-300						#(60*2.5)*2
t<-seq(1,n+1)

c<-1/15;	
M_n<-ceiling(c*sqrt(n))		#pre-average window length jump return
#120, 122 (30sec before and after tau)
Jstart_time=120			#119 1min before 120 30 sec before 8:30am==30600 sec after midnight\\ FOMC 2:00pm== 50400 sec after midnight 
Jwindow=122				#123 1 min #125 2min cut out #131 # 5 min cut out 

alpha<-0.01					#test level

event<-"12bonds_CPIdates" #"12bonds_NFPdates"#"12bonds_FOMCdates"##		

days<-list.files(path=paste(getwd(),event,sep="/")) #list of all available bonds

#Output
Y1J<-array(NA,dim=c(length(days),2));Y2J<-array(NA,dim=c(length(days),2));SJ<-array(NA,dim=c(length(days),2));
IUTchi2<-array(NA,dim=c(length(days),2)); LM<-array(NA,dim=c(length(days),2)); IUTBon<-array(NA,dim=c(length(days),2));
obs<-array(NA,dim=c(length(days),6)); rho<-array(NA,dim=c(length(days),2)) 

for(d in 1:length(days)){
	#if(d==28) next		#CPI, FOMC

	bond<-list.files(path=paste(getwd(),event,days[d],sep="/"))	
	#load Bonds order on folder: U.S. bonds 1=1Y, 2=2Y, 3=5Y, 4=10Y, 5=15Y, 6=20Y 
		N5<-read.csv(paste(getwd(),event,days[d],bond[6],sep="/"), header = FALSE, sep = ",", dec = ".") # 2 # 6
		N10<-read.csv(paste(getwd(),event,days[d],bond[4],sep="/"), header = FALSE, sep = ",", dec = ".") # 6 #13
	
	#load TIPS	order on folder: 8=1Y, 9=2Y, 10=5Y, 11=10Y, 12=15Y, 13=20Y 
		R5<-read.csv(paste(getwd(),event,days[d],bond[2],sep="/"), header = FALSE, sep = ",", dec = ".") # 9 # 2
		R10<-read.csv(paste(getwd(),event,days[d],bond[2],sep="/"), header = FALSE, sep = ",", dec = ".") # 13 # 9

	#calculate spread
		S<-N5[,4]-R5[,4]; S<-cbind(S, N10[,4]-R10[,4]) 

	#calculate covola
		Y1<-cbind(N5[,4],R5[,4]); 		rownames(Y1)<-t
		Y2<-cbind(N10[,4],R10[,4]); 	rownames(Y2)<-t

		obs[d,1:2]<-colSums(!diff(Y1)==0); obs[d,3:4]<-colSums(!diff(Y2)==0); obs[d,5:6]<-colSums(!diff(S)==0)

		u1<-setData(Y1); u2<-setData(Y2); 				
		Sig1<-cce(u1,method="PTHY")				 	 #pre-average MRC of Christensen et al.(2010)
		Sig2<-cce(u2,method="PTHY")

		eta_hat21<-1/(2*n)*t(diff(Y1))%*%diff(Y1) 			#noise estimation
		eta_hat22<-1/(2*n)*t(diff(Y2))%*%diff(Y2) 	

	#Event return \tau
		Y1J[d,]<-colMeans(Y1[(Jwindow:(Jwindow+M_n)),])-colMeans(Y1[((Jstart_time-M_n):Jstart_time),])
		Y2J[d,]<-colMeans(Y2[(Jwindow:(Jwindow+M_n)),])-colMeans(Y2[((Jstart_time-M_n):Jstart_time),])
		SJ[d,]<-colMeans(S[(Jwindow:(Jwindow+M_n)),])-colMeans(S[((Jstart_time-M_n):Jstart_time),])

	#Asymptotic variance of pre-average returns (test statistic)	
		Gamma1<-2/3*Sig1$covmat*c+2/c*eta_hat21	
		Gamma2<-2/3*Sig2$covmat*c+2/c*eta_hat22

		rho[d,1]<-Gamma1[1,2]/sqrt(Gamma1[1,1]*Gamma1[2,2])
		rho[d,2]<-Gamma2[1,2]/sqrt(Gamma2[1,1]*Gamma2[2,2])

	#IUT for spread jumps
		
			#Step 1
				chi2test1<-(1-pchisq(sqrt(max(obs[d,c(1,2)]))*Y1J[d,]%*%solve(Gamma1)%*%Y1J[d,], df=2))
				chi2test2<-(1-pchisq(sqrt(max(obs[d,c(3,4)]))*Y2J[d,]%*%solve(Gamma2)%*%Y2J[d,], df=2))

				pLM1<-1- pnorm(obs[d,c(1,2)]^(1/4)*abs(Y1J[d,])/sqrt(diag(Gamma1)) )#p-value of Lee-Mykland (2012) Lemma1 test
				pLM2<-1- pnorm(obs[d,c(3,4)]^(1/4)*abs(Y2J[d,])/sqrt(diag(Gamma2)) )	

				#chi2test1<-(1-pchisq(sqrt(n)*Y1J[d,]%*%solve(Gamma1)%*%Y1J[d,], df=2))
				#chi2test2<-(1-pchisq(sqrt(n)*Y2J[d,]%*%solve(Gamma2)%*%Y2J[d,], df=2))

				#pLM1<-1- pnorm(n^(1/4)*abs(Y1J[d,])/sqrt(diag(Gamma1)) )#p-value of Lee-Mykland (2012) Lemma1 test
				#pLM2<-1- pnorm(n^(1/4)*abs(Y2J[d,])/sqrt(diag(Gamma2)) )

			#Step 2
				pLMS1<-1- pnorm(obs[d,5]^(1/4)*abs(SJ[d,1])/sqrt(Gamma1[1,1]+Gamma1[2,2]-2*Gamma1[1,2])) #Lee-Mykland spread
				pLMS2<-1- pnorm(obs[d,6]^(1/4)*abs(SJ[d,2])/sqrt(Gamma2[1,1]+Gamma2[2,2]-2*Gamma2[1,2])) 

				#pLMS1<-1- pnorm(n^(1/4)*abs(SJ[d,1])/sqrt(Gamma1[1,1]+Gamma1[2,2]-2*Gamma1[1,2])) #Lee-Mykland spread
				#pLMS2<-1- pnorm(n^(1/4)*abs(SJ[d,2])/sqrt(Gamma2[1,1]+Gamma2[2,2]-2*Gamma2[1,2])) 

				if( (pLMS1<(alpha/2)) & chi2test1<alpha){IUTchi2[d,1]<-1}else{IUTchi2[d,1]<-0}   #BEI test is two sided: alpha/2
				if( (pLMS2<(alpha/2)) & chi2test2<alpha){IUTchi2[d,2]<-1}else{IUTchi2[d,2]<-0}   

				if( (pLMS1<(alpha/2)) & (pLM1[1]<((alpha/2)/2) | pLM1[2]<((alpha/2)/2) ) ){IUTBon[d,1]<-1}else{IUTBon[d,1]<-0}
				if( (pLMS2<(alpha/2)) & (pLM2[1]<((alpha/2)/2) | pLM2[2]<((alpha/2)/2) ) ){IUTBon[d,2]<-1}else{IUTBon[d,2]<-0}

				if( (pLMS1<(alpha/2))){LM[d,1]<-1}else{LM[d,1]<-0}   #BEI test is two sided: alpha/2
				if( (pLMS2<(alpha/2))){LM[d,2]<-1}else{LM[d,2]<-0} 

}

colSums(LM,na.rm=TRUE)/length(days); 
colSums(IUTBon,na.rm=TRUE)/length(days); 
colSums(IUTchi2,na.rm=TRUE)/length(days)

colSums(LM,na.rm=TRUE); 
colSums(IUTBon,na.rm=TRUE); 
colSums(IUTchi2,na.rm=TRUE)

#conf<-LM-IUTchi2
conf<-LM-IUTBon

Gamma1[1,2]/sqrt(Gamma1[1,1]*Gamma1[2,2])
Gamma2[1,2]/sqrt(Gamma2[1,1]*Gamma2[2,2])

mean(abs(SJ[which(IUTchi2[,1]==1),1]))*100;
mean(abs(SJ[which(IUTchi2[,2]==1),2]))*100

plot(Y2J[which(conf[,2]==1),1],Y2J[which(conf[,2]==1),2] ,pch=19,col="orange", xlim=c(-0.06,0.06), ylim=c(-0.06,0.06) )
points(Y2J)					#, 
points(Y2J[which(IUTBon[,2]==1),] ,pch=19,col="red")
abline(0,1)

plot(Y1J[which(conf[,1]==1),1],Y1J[which(conf[,1]==1),2] ,pch=19,col="orange", xlim=c(-0.06,0.06), ylim=c(-0.06,0.06) )
points(Y1J)					#, 
points(Y1J[which(IUTchi2[,1]==1),] ,pch=19,col="red")
abline(0,1)

plot(SJ[,1])
points(which(IUTchi2[,1]==1),SJ[which(IUTchi2[,1]==1),1],pch=19,col="green" )
points(which(conf[,1]==1),SJ[which(conf[,1]==1),1],pch=19,col="red" )

plot(SJ[,2])
points(which(IUTchi2[,2]==1),SJ[which(IUTchi2[,2]==1),2],pch=19,col="green" )
points(which(conf[,2]==1),SJ[which(conf[,2]==1),2],pch=19,col="red" )

### Save for combined figure chi2

Y2J_BEI_MO<-Y2J
Y2J_BEI_MO_BoJ<-Y2J[which(IUTBon[,2]==1),] 
Y2J_BEI_MO_conf<-Y2J[which(conf[,2]==1),]

Y2J_S_MO<-Y2J
Y2J_S_MO_BoJ<-Y2J[which(IUTBon[,2]==1),] 
Y2J_S_MO_conf<-Y2J[which(conf[,2]==1),]

Y2J_BEI_INF<-Y2J
Y2J_BEI_INF_BoJ<-Y2J[which(IUTBon[,2]==1),] 
Y2J_BEI_INF_conf<-Y2J[which(conf[,2]==1),]

Y2J_S_INF<-Y2J
Y2J_S_INF_BoJ<-Y2J[which(IUTBon[,2]==1),] 
Y2J_S_INF_conf<-Y2J[which(conf[,2]==1),]

Y2J_BEI_NFP<-Y2J
Y2J_BEI_NFP_BoJ<-Y2J[which(IUTBon[,2]==1),] 
Y2J_BEI_NFP_conf<-Y2J[which(conf[,2]==1),]

Y2J_S_NFP<-Y2J
Y2J_S_NFP_BoJ<-Y2J[which(IUTBon[,2]==1),] 
Y2J_S_NFP_conf<-Y2J[which(conf[,2]==1),]

z<-1.2
plot(Y2J_BEI_MO,pch=5, xlim=c(-0.1,0.1), ylim=c(-0.1,0.1), cex=z,cex.axis=1.2 )
abline(0,1)
points(Y2J_BEI_INF,pch=1, cex=z)
points(Y2J_BEI_NFP,pch=2, cex=z)
		
points(Y2J_BEI_MO_BoJ,pch=18, col="green", cex=z)
points(Y2J_BEI_INF_BoJ,pch=16, col="green", cex=z)
points(Y2J_BEI_NFP_BoJ,pch=17, col="green", cex=z)

points(Y2J_BEI_MO_conf,pch=18, col="red", cex=z)
points(Y2J_BEI_INF_conf,pch=16, col="red", cex=z)
points(Y2J_BEI_NFP_conf,pch=17, col="red", cex=z)
			
legend("topleft", legend=c("Monetary policy", "Inflation", "Employment", "Spread jump", "Incoherent"),
       col=c("black", "black","black", "green", "red"), lty=0, pch=c(5,1,2,15,15,15), lwd=1, cex=1.3,box.lty=1)

z<-1.2
plot(Y2J_S_MO,pch=5, xlim=c(-0.1,0.1), ylim=c(-0.1,0.1), cex=z,cex.axis=1.2 )
abline(0,1)
points(Y2J_S_INF,pch=1, cex=z)
points(Y2J_S_NFP,pch=2, cex=z)
		
points(Y2J_S_MO_BoJ,pch=18, col="green", cex=z)
points(Y2J_S_INF_BoJ,pch=16, col="green", cex=z)
points(Y2J_S_NFP_BoJ,pch=17, col="green", cex=z)

points(Y2J_S_MO_conf,pch=18, col="red", cex=z)
points(Y2J_S_INF_conf[1],Y2J_S_INF_conf[2],pch=16, col="red", cex=z)
points(Y2J_S_NFP_conf[1],Y2J_S_NFP_conf[2],pch=17, col="red", cex=z)
			
legend("topleft", legend=c("Monetary policy", "Inflation", "Employment", "Spread jump", "Incoherent"),
       col=c("black", "black","black", "green", "red"), lty=0, pch=c(5,1,2,15,15,15), lwd=1, cex=1.3,box.lty=1)