setwd() 

library(stringr)

#Jan2020--June2020: 124 days		#kick-out 4.Jun20 trading until 11am
# 1-3Y, 3-7Y, 7-10Y, 10-20Y, 20Y+

no<-125

P<-array(NA,dim=c(390*125,5))					#1 minute mid quotes of prices of 
files<-seq(1,125)

bonds<-list.files(path=getwd())

# extract message files of all stocks across days
for(i in 1:length(bonds)){
	
	csvfiles<-list.files(path=paste(getwd(),bonds[i],sep="/"))
	
	for (day in 1:length(csvfiles)){ 
		S<-read.csv(file=paste(getwd(),bonds[i],csvfiles[day],sep="/"),header=TRUE) #dat
		#P[((day-1)*390+1):(day*390) ,i]<-(S[1:390,2]-S[1:390,4])/2+S[1:390,4]/10000
		P[((day-1)*390+1):(day*390) ,i]<-S[1:390,2]/10000 # askprice only
	}
}	

###################################################################
# 13 minutes integrated covola
#27 weeks in 2020 Jan 2--June 30
#390 1 minute prices per day; 

weekdays<-c(2,5,5,4,5,5,5,4,5,5,5,5,5,5,4,5,5,5,5,5,5,4,4,5,5,5,2) 
block<-13*2 # smallest block length
bpd<-390/block #blocks per day with 13 min length

#Pick specific bonds
N<-P[,c(2,1,4)]	#3-7Y, 7-10Y,10-20Y
#N<-P[,c(3,2,1)]
#N<-P[,c(3,2,1,4,5)]		#1-3Y, 3-7Y, 7-10Y,10-20Y, 20+
dim<-dim(N)[2]		#number of bonds

#estimate vol of vol for rank test and estimate
#!!! Set week befor running the code below !!! gives dwN
#K2<-1950/26=75 	#50-1#floor((3600*13)^(1/3)*14.43)	#n^(1/3)*c1
opb2<-65	# 39 obs per block, 5 week days
bpd2<-390/opb2

#weekly integrated eigenvalues (13min) and eigenvalues of weekly integrated covola 
inteig<-array(NA,dim=c(dim,length(weekdays))); eigint<-array(NA,dim=c(dim,length(weekdays)))
intSig<-array(NA,dim=c(dim^2,length(weekdays)))
NQV_hat<-seq(1,length(weekdays));varNQV<-seq(1,length(weekdays)); kappa_gap<-array(NA,dim=c(2,length(weekdays)));

for(week in 1:length(weekdays)){								
		wN<-N[(((sum(weekdays[1:week])-weekdays[week])*390)+1):(sum(weekdays[1:week])*390),]
		for(i in 1:dim(wN)[1]){for(j in 1:dim(wN)[2]){if(wN[i,j]<0){wN[i,j]<-wN[i-1,j]}}}#there are negative prices
		#dwN<-diff(log(wN))*100			#return
		dwN<-t(t(diff(wN))*c(1/6,1/7,1/20))		#c(1/3,1/3,1/7,1/20,1/25))				#yield chage within day in percent 3Bonds: c(1/3,1/7,1/20))

			#trunkate jumps with weekly vol
			for(i in 1:dim(dwN)[2]){					#kick out jumps
				#dwN[which(abs(dwN[,i])>1.25),i]<-0
				dwN[which(abs(dwN[,i])>mad(dwN[which(abs(dwN[,i])>0),i])*4),i]<-0 #replace >4xmeadian absolute deviation
				#dwN[,i]<-dwN[,i]/sqrt(sum(dwN[,i]^2))	#standardized returns
			}

		#########################

		K<-weekdays[week]*390/block-1 # by taking first difference wee lose one obs; we throw away the last block on Fridy
		
		Sig<-array(NA,dim=c(dim^2,K))
		beig<-array(NA,dim=c(dim,K))
		bvec1<-array(NA,dim=c(dim,K)); bvec2<-array(NA,dim=c(dim,K)); bvec3<-array(NA,dim=c(dim,K)); 
		for(k in 1:K){	# covola on 13 min block  
			dNb<-dwN[((k-1)*block+1):(k*block), ];
			Sigb<-(5*bpd)*t(dNb)%*%dNb; 			#scaled to 1 week
			Sig[,k]<-c(Sigb);  
			beig[,k]<-eigen(Sigb)$values	#eigenvalues on smallest (13min) blocks
			bvec1[,k]<-(eigen(Sigb)$vectors)[,1]	
			bvec2[,k]<-(eigen(Sigb)$vectors)[,2]	
			bvec3[,k]<-(eigen(Sigb)$vectors)[,3]	
		}
		
		inteig[,week]<-rowMeans(beig)					#integrated eigenvalues
		eigint[,week]<-eigen(matrix(rowMeans(Sig),nrow=dim))$values #eigenvalues of integrated covol
		intSig[,week]<-rowMeans(Sig)

#estimate vol of vol for rank test and estimate
K2<-weekdays[week]*390/(opb2)-1
#Realized covola on large n^1/3 blocks
	Sig2<-array(NA,dim=c(dim^2,K2))	
	for(k in 1:K2){	
			dNb<-dwN[((k-1)*opb2+1):(k*opb2), ];
			Sigb<-(5*bpd2)*t(dNb)%*%dNb; 
			Sig2[,k]<-c(Sigb);  
	}
	
	#Sig2[,1]<-rep(0,dim)

	# volofvol NQV (quadratic variation), non overlapping blocks
	NQVk<-seq(1,K2/3)
	for(i in 1:(K2/3)){
			NQVk[i]<-max(abs(eigen(matrix((Sig2[,(i-1)*3+3]-2*Sig2[,(i-1)*3+2]+Sig2[,(i-1)*3+1]),nrow=dim))$values))^2  # Alternative zu sum d\sig^2	
	}
	NQV_hat[week]<-3/2*sum(NQVk)		#/55

	# covariance of the volatility of volatility volofvol^2
		rho_4<- K2*3/4*sum(NQVk^2)
		rho_22<-seq(1,K2/6)
		for(i in 1:(K2/6)){
			rho_22[i]<-max(abs(eigen(matrix(Sig2[,(i-1)*6+3]-2*Sig2[,(i-1)*6+2]+Sig2[,(i-1)*6+1],nrow=dim))$values))^2*
					max(abs(eigen(matrix(Sig2[,(i-1)*6+6]-2*Sig2[,(i-1)*6+5]+Sig2[,(i-1)*6+4],nrow=dim))$values))^2		}
		rho_22<-6*K2/4*sum(rho_22)

		# asyptotic variance of quadratic covariation of covolatility
		varNQV[week]<-3/K2*(rho_4-rho_22)
	# asymptotic critical value
	kappa_gap[1,week]<-1/(12*K*max(median(beig[2,],na.rm=TRUE),0.0001))*(NQV_hat[week]+varNQV[week]^(1/2)*qnorm(0.95,mean=0, sd = 1) )
	kappa_gap[2,week]<-1/(12*K*max(median(beig[3,],na.rm=TRUE),0.0001))*(NQV_hat[week]+varNQV[week]^(1/2)*qnorm(0.95,mean=0, sd = 1) )
	#kappa<-sqrt(NQV_hat)*K^(-1/2)*(1+(1/K2*varNQV)^(1/2)*qnorm(0.95,mean=0, sd = 1))

}

#######################################################################

#normalize eigenvector (sign of first entry is >0)
for(i in 1:dim(bvec1)[2]){if(bvec1[1,i]<0){bvec1[,i]<-bvec1[,i]*(-1)}}
for(i in 1:dim(bvec2)[2]){if(bvec1[1,i]<0){bvec2[,i]<-bvec2[,i]*(-1)}}
for(i in 1:dim(bvec3)[2]){if(bvec3[1,i]<0){bvec3[,i]<-bvec3[,i]*(-1)}}

# Eigenvector plots CYCLE
theta=runif(1000,min=0, max=2*pi)
x=cos(theta)
y=sin(theta)
plot(x,y, type="l",col="gray")
 xspline(bvec1[2,],bvec1[3,],shape=-1,border="blue")
lines(bvec1[2,],bvec1[3,],type="p", col="blue")

 xspline(bvec2[2,],bvec2[3,],shape=-1,border="red")
lines(bvec2[2,],bvec2[3,],type="p", col="red")

 xspline(bvec1[2,],bvec1[3,],shape=-1,border="blue")

#old 
#plot(inteig[3,]/colSums(inteig)*100,ylim=c(2,30),xlim=c(1,27),pch=24, lwd=2, bg="red",col="black",cex=1.5)
#par(new=TRUE)
#plot(c(1,3,4,7,9,11:25,27),inteig[2,c(1,3,4,7,9,11:25,27)]/colSums(inteig[,c(1,3,4,7,9,11:25,27)])*100, ylim=c(2,30),xlim=c(1,27), pch=21, lwd=2, bg="green",col="black",cex=1.5 , ann=FALSE,xaxt='n'
,yaxt='n', bty="n")
#par(new=TRUE)
#plot(c(2,5,6,8,10,26),inteig[2,c(2,5,6,8,10,26)]/colSums(inteig[,c(2,5,6,8,10,26)])*100, ylim=c(2,30),xlim=c(1,27), pch=21, lwd=2, bg="red",col="black",cex=1.5 , ann=FALSE,xaxt='n'
,yaxt='n', bty="n")

plot(inteig[3,]/colSums(inteig)*100,ylim=c(2,30),xlim=c(1,27),pch=24, lwd=2, bg="red",col="black",cex=2)
par(new=TRUE)
plot(c(3,4,7,12,12:23,27),inteig[2,c(3,4,7,12,12:23,27)]/colSums(inteig[,c(3,4,7,12,12:23,27)])*100, ylim=c(2,30),xlim=c(1,27), pch=21, lwd=2, bg="green",col="black",cex=2 , ann=FALSE,xaxt='n'
,yaxt='n', bty="n")
par(new=TRUE)
plot(c(1,2,5,6,8:11,24:26),inteig[2,c(1,2,5,6,8:11,24:26)]/colSums(inteig[,c(1,2,5,6,8:11,24:26)])*100, ylim=c(2,30),xlim=c(1,27), pch=21, lwd=2, bg="red",col="black",cex=2 , ann=FALSE,xaxt='n'
,yaxt='n', bty="n")

lines(eigint[3,]/colSums(eigint)*100,type="l", lty=3,lwd=2,col="gray")
lines(eigint[2,]/colSums(eigint)*100,type="l", lty=3,col="gray",lwd=2)

plot(beig[2,]/colSums(beig)*100,type="o")

plot((eigint[3,]-inteig[3,])/inteig[3,]*100,type='l',ylim=c(-10,110),col='red')
lines((eigint[2,]-inteig[2,])/inteig[2,]*100,col="blue")
legend("topleft", legend=c("2nd eigenvalue","3rd eigenvalue"),
       col=c("blue", "red"), lty=c(1,1), cex=1)

plot(inteig[1,],log='y',type='l',ylim=c(0.0004,0.55))
lines(eigint[1,],col='darkgray')
lines(inteig[2,],col='blue')
lines(eigint[2,],col='darkgray')
lines(inteig[3,],col='red')
lines(eigint[3,],col='darkgray')

eigint[3,]/inteig[3,]	# take week 10 (1st week of March) faktor 1.24; week 12 (3rd week of March) factor 3.01

plot(inteig[,12])
lines(eigint[,12])
inteig[,12]/sum(inteig[,12])
eigint[,12]/sum(eigint[,12])

# run above with week<-23

plot(rowMeans(beig))
lines(eigen(matrix(rowMeans(Sig),nrow=dim))$values)

plot(inteig[,12]/sum(inteig[,12]))
lines(eigint[,12]/sum(eigint[,12]))

plot(t(bvec1[2:3,]),type='o')
plot(t(bvec2[2:3,]))

plot(inteig[,1],type="o", ylim=c(0,0.1))
for(i in 2:27){lines(inteig[,i],type="o",col=i)}

#######################################################################
#PLOT week mit varying blocklength:
week<-7
opb<-13
		
		wN<-N[(((sum(weekdays[1:week])-weekdays[week])*390)+1):(sum(weekdays[1:week])*390),]
		for(i in 1:dim(wN)[1]){for(j in 1:dim(wN)[2]){if(wN[i,j]<0){wN[i,j]<-wN[i-1,j]}}}#there are negative prices
		#dwN<-diff(log(wN))*100			#return
		dwN<-t(t(diff(wN))*c(1/6,1/7,1/20))		#c(1/3,1/3,1/7,1/20,1/25))				#yield chage within day in percent 3Bonds: c(1/3,1/7,1/20))

			#trunkate jumps with weekly vol
			for(i in 1:dim(dwN)[2]){					#kick out jumps
				#dwN[which(abs(dwN[,i])>1.25),i]<-0
				dwN[which(abs(dwN[,i])>mad(dwN[which(abs(dwN[,i])>0),i])*4),i]<-0 #replace >4xmeadian absolute deviation
				#dwN[,i]<-dwN[,i]/sqrt(sum(dwN[,i]^2))	#standardized returns
			}

		dwN<-rbind(dwN,dwN[1949,])
		#########################

		K<-weekdays[week]*390/block # by taking first difference wee lose one obs; we throw away the last block on Fridy
		
		Sig<-array(NA,dim=c(dim^2,K))
		beig<-array(NA,dim=c(dim,K))
		for(k in 1:K){	# covola on 13 min block  
			dNb<-dwN[((k-1)*block+1):(k*block), ];
			Sigb<-(5*bpd)*t(dNb)%*%dNb; 			#scaled to 1 week
			Sig[,k]<-c(Sigb);  
			beig[,k]<-eigen(Sigb)$values	#eigenvalues on smallest (13min) blocks
		}

opb<-K/seq(1,K); opb<-opb[opb%%1==0];opb<-rev(opb[1:(length(opb)-1)]);
inteig<-array(NA,dim=c(dim,length(opb)))

for(i in 1:length(opb)){
	keig<-array(NA,dim=c(dim,K/opb[i]))
	for(j in 1:(K/opb[i])){
		Sigb<-matrix(rowMeans(Sig[,((j-1)*opb[i]+1):(j*opb[i])]),nrow=dim)
		keig[,j]<-eigen(Sigb)$values
	}
	
	inteig[,i]<-rowMeans(keig)
}
inteig<-cbind(rowMeans(beig),inteig)

#PERCENTAGE PCA
plot(c(1,opb[1:8])/150,inteig[1,1:9]/colSums(inteig[,1:9])*100,ylim=c(4,100),xlim=c(0,0.2),type="o",log="y",col="blue",lwd=2)
lines(c(1,opb[1:8])/150,inteig[2,1:9]/colSums(inteig[,1:9])*100,type="o",col="red",lwd=2)
lines(c(1,opb[1:8])/150,inteig[3,1:9]/colSums(inteig[,1:9])*100,type="o",col="green",lwd=2)

lines(c(1,opb[1:8])/150,inteig[1,1:9]/colSums(inteig[,1:9])*100,type="o",col="blue" ,lwd=2 )
lines(c(1,opb[1:8])/150,inteig[2,1:9]/colSums(inteig[,1:9])*100,type="o",col="red",lwd=2)
lines(c(1,opb[1:8])/150,inteig[3,1:9]/colSums(inteig[,1:9])*100,type="o",col="green",lwd=2)

#Absolute Eigval
plot(c(1,opb[1:8])/150,inteig[1,1:9],ylim=c(0.00025,0.02),xlim=c(0,0.2),type="l",log="y",col="blue")
lines(c(1,opb[1:8])/150,inteig[2,1:9],type="l",col="red")
lines(c(1,opb[1:8])/150,inteig[3,1:9],type="l",col="green")

lines(c(1,opb[1:8])/150,inteig[1,1:9],type="l",col="blue"  )
lines(c(1,opb[1:8])/150,inteig[2,1:9],type="l",col="red")
lines(c(1,opb[1:8])/150,inteig[3,1:9],type="l",col="green")