library(quantreg)
library(dplyr)
library(readxl)
library(pastecs)


df<-read_excel("data_paper3.xlsx", sheet ="data")
names(df)[names(df) == 'BTC'] <- "Bitcoin_price"
df[df == ""] = NA
cols.num <- c(4:12)
df[cols.num] <- sapply(df[cols.num],as.numeric)
df$return <- log(df$Bitcoin_price) - lag(log(df$Bitcoin_price),1)
df <- df[complete.cases(df),]

#Descriptive statistics--------------------------------------------------------------
stat.desc(df[-1])

#Graphical description----------------------------------------------------------------
library('ggplot2')
plot_return <- ggplot(df, aes(x = Time)) +
  geom_line(aes(y = return), colour = "blue", linetype = 1, size = 0.5) 
  +ylab("") + xlab("") + ggtitle("Bitcoin's returns")
  +theme_bw()+theme_linedraw()
  +theme(plot.title=element_text(size=12,hjust=0.5,face="bold"))

visual1 <- gather(df,variable,EPU,df.Canada:df.Singapore)

ggplot(visual1,aes(x=df.Dateformat,y=EPU,colour=variable,group=variable))+geom_line()

#Estimation#######################################################################

#Quantile on quantile Function--------------------------------------------------------------------
install.packages("quantreg")
library(quantreg)
taus<-seq(0.1,0.9,0.1)

lprq <-  function(x, y, h, m=9 , tau=0.5)
{
  xx <- seq(min(x),max(x),length=m)
  fv <- xx
  dv <- xx
  for(i in 1:length(xx)) {
    z <- x - xx[i]
    wx <- dnorm(z/h)
    r <- rq(y~z, weights=wx, tau=tau, ci=FALSE)
    fv[i] <- r$coef[1.]
    dv[i] <- r$coef[2.]
  }
  data.frame(dv = dv)
}

#The impact of EPU on Bitcoin returns-----------------------------------------------------------
#Create a matrix to save the QQR estimates
QQR.coef_Canada <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_France <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_Germany <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_Japan <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_UK <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_US <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_China <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_Russia <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_Netherlands <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.coef_Singapore <- as.data.frame(matrix(0, ncol = 9, nrow = 9))

# Run QQR in a loop and save estimates in matrix "QQR.coef"
#Note: 0.05 in below loop is the bandwidth that can be adjusted
#Canada
library(tidyr)
library(plot3D)

for (i in 1:9){
  a <- lprq(log(df$Canada),df$return,0.06,tau=taus[i])
  QQR.coef_Canada[,i]<-a
}

#C2: persp3D()
library(plot3D)
library(rgl)
library(RColorBrewer)
require(RColorBrewer)

x= seq(0.1,0.9,0.1)
y=seq(0.1,0.9,0.1)
#Canada
z <- as.matrix(QQR.coef_Canada)
nbcol = 100
color = rev(rainbow(nbcol, start = 0/6, end = 4/6))
zcol  = cut(z, nbcol)
persp3D(x,y,z, theta=310, phi=20, shade=0,d =2, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="Canada")
#US
for (i in 1:9){
  a<-lprq(log(df$US),df$return,0.1,tau=taus[i])
  QQR.coef_US[,i]<-a
}
z <- as.matrix(QQR.coef_US)
persp3D(x,y,z, theta=325, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="US")
#UK
for (i in 1:9){
  a<-lprq(log(df$UK),df$return,0.1,tau=taus[i])
  QQR.coef_UK[,i]<-a
}
z <- as.matrix(QQR.coef_UK)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="UK")
#France
for (i in 1:9){
  a<-lprq(log(df$France),df$return,0.05,tau=taus[i])
  QQR.coef_France[,i]<-a
}
z <- as.matrix(QQR.coef_France)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="France")
#Germany
for (i in 1:9){
  a<-lprq(log(df$Germany),df$return,0.05,tau=taus[i])
  QQR.coef_Germany[,i]<-a
}
z <- as.matrix(QQR.coef_Germany)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="Germany")
#Russia
for (i in 1:9){
  a<-lprq(log(df$Russia),df$return,0.09,tau=taus[i])
  QQR.coef_Russia[,i]<-a
}
z <- as.matrix(QQR.coef_Russia)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="Russia")

#Japan
for (i in 1:9){
  a<-lprq(log(df$Japan),df$return,0.05,tau=taus[i])
  QQR.coef_Japan[,i]<-a
}
z <- as.matrix(QQR.coef_Japan)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="Japan")
#China
for (i in 1:9){
  a<-lprq(log(df$`Mainland China`),df$return,0.05,tau=taus[i])
  QQR.coef_China[,i]<-a
}
z <- as.matrix(QQR.coef_China)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="China")
# Other countries (Singapore and Netherlands)
df2<-read_excel("data_paper3.xlsx", sheet ="data2")
names(df2)[names(df2) == 'BTC'] <- "Bitcoin_price"
df2$return<-log(df2$Bitcoin_price)-log(lag(df2$Bitcoin_price,1))
df2 <- df2[complete.cases(df2),]
#Singapore
for (i in 1:9){
  a<-lprq(log(df2$Singapore),df2$return,0.06,tau=taus[i])
  QQR.coef_Singapore[,i]<-a
}
x = seq(0.1,0.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.coef_Singapore)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="Singapore")
#Netherlands
for (i in 1:9){
  a<-lprq(log(df2$Netherlands),df2$return,0.05,tau=taus[i])
  QQR.coef_Netherlands[,i]<-a
}
z <- as.matrix(QQR.coef_Netherlands)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of returns", zlab= "beta",main="Netherlands")

#Quantile regression estimation-------------------------------------------------------------
fit_us <-rq(df$return~log(df$US),tau=taus)
plot(summary(fit_us), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "US", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_uk <- rq(df$return~log(df$UK),tau=taus)
plot(summary(fit_uk), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "UK", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_cad <- rq(df$return~log(df$Canada),tau=taus)
plot(summary(fit_cad), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Canada", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_france <- rq(df$return~log(df$France),tau=taus)
plot(summary(fit_france), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "France", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_germany <- rq(df$return~log(df$Germany),tau=taus)
plot(summary(fit_germany), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Germany", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_russia <- rq(df$return~log(df$Russia),tau=taus)
plot(summary(fit_russia), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Russia", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_japan <- rq(df$return~log(df$Japan),tau=taus)
plot(summary(fit_japan), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Japan", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_china <- rq(df$return~log(df$`Mainland China`),tau=taus)
plot(summary(fit_china), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "China", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_sing <- rq(df2$return~log(df2$Singapore),tau=taus)
plot(summary(fit_sing), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Singapore", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
fit_netherlands <- rq(df2$return~log(df2$Netherlands),tau=taus)
plot(summary(fit_netherlands), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Netherlands", xlab = "quantile of Bitcoin return", 
     ylab = "coefficient", cex = 1, pch = 19)
# Impact of EPU on Bitcoin's volatility-------------------------------------------------------------------
df_vol<-read_excel("data_paper3.xlsx", sheet ="volatility")
df_vol1 <- df_vol[,c(-11,-12)]
df_vol2 <- df_vol[,c(1,2,11,12)]
df_vol2 <- df_vol2[complete.cases(df_vol2),]
#volaitlity graph
plot_volatility <- ggplot(df_vol1, aes(x = Time)) +
  geom_line(aes(y = volatility), colour = "blue", linetype = 1, size = 0.5) +
  ylab("") + xlab("") +ggtitle("Bitcoin's volatility")+ theme_bw()+theme_linedraw()+theme(plot.title=element_text(size=12,hjust=0.6,face="bold"))
#Quantile estimation
vol_us <- rq(df_vol1$volatility~log(df_vol1$US),tau=taus)
plot(summary(vol_us), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "US", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_uk <- rq(df_vol1$volatility~log(df_vol1$UK),tau=taus)
plot(summary(vol_uk), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "UK", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_canada <- rq(df_vol1$volatility~log(df_vol1$Canada),tau=taus)
plot(summary(vol_canada), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Canada", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_france <- rq(df_vol1$volatility~log(df_vol1$France),tau=taus)
plot(summary(vol_france), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "France", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_germany <- rq(df_vol1$volatility~log(df_vol1$Germany),tau=taus)
plot(summary(vol_germany), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Germany", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_russia <- rq(df_vol1$volatility~log(df_vol1$Russia),tau=taus)
plot(summary(vol_canada), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Russia", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_japan <- rq(df_vol1$volatility~log(df_vol1$Japan),tau=taus)
plot(summary(vol_japan), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Japan", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_china <- rq(df_vol1$volatility~log(df_vol1$`Mainland China`),tau=taus)
plot(summary(vol_china), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "China", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_singapore <- rq(df_vol2$volatility~log(df_vol2$Singapore),tau=taus)
plot(summary(vol_singapore), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Singapore", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
vol_netherlands <- rq(df_vol2$volatility~log(df_vol2$Netherlands),tau=taus)
plot(summary(vol_netherlands), parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "Netherlands", xlab = "quantile of Bitcoin volatility", 
     ylab = "coefficients", cex = 1, pch = 19)
#Quantile on quantile estimation

QQR.vol_US <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_UK <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_Canada <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_France <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_Germany <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_Russia <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_Japan <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_China <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_Singapore <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
QQR.vol_Netherlands <- as.data.frame(matrix(0, ncol = 9, nrow = 9))
#US
for (i in 1:9){
  a<-lprq(log(df_vol1$US),df_vol1$volatility,0.05,tau=taus[i])
  QQR.vol_US[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_US)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="USA")
#UK
for (i in 1:9){
  a<-lprq(log(df_vol1$UK),df_vol1$volatility,0.1,tau=taus[i])
  QQR.vol_UK[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_UK)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="UK")
#Canada
for (i in 1:9){
  a<-lprq(log(df_vol1$Canada),df_vol1$volatility,0.1,tau=taus[i])
  QQR.vol_Canada[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_Canada)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="Canada")
#France
for (i in 1:9){
  a<-lprq(log(df_vol1$France),df_vol1$volatility,0.09,tau=taus[i])
  QQR.vol_France[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_France)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="France")
#Germany
for (i in 1:9){
  a<-lprq(log(df_vol1$Germany),df_vol1$volatility,0.07,tau=taus[i])
  QQR.vol_Germany[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_Germany)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="Germany")
#Russia
for (i in 1:9){
  a<-lprq(log(df_vol1$Russia),df_vol1$volatility,0.1,tau=taus[i])
  QQR.vol_Russia[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_Russia)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="Russia")
#Japan
for (i in 1:9){
  a<-lprq(log(df_vol1$Japan),df_vol1$volatility,0.1,tau=taus[i])
  QQR.vol_Japan[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_Japan)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="Japan")
#China
for (i in 1:9){
  a<-lprq(log(df_vol1$`Mainland China`),df_vol1$volatility,0.09,tau=taus[i])
  QQR.vol_China[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_China)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="China")
#Singapore
for (i in 1:9){
  a<-lprq(log(df_vol2$Singapore),df_vol2$volatility,0.1,tau=taus[i])
  QQR.vol_Singapore[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_Singapore)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="Singapore")
#Netherlands
for (i in 1:9){
  a<-lprq(log(df_vol2$Netherlands),df_vol2$volatility,0.05,tau=taus[i])
  QQR.vol_Netherlands[,i]<-a
}
x = seq(0.1,.9,0.1)
y = seq(0.1,0.9,0.1)
z <- as.matrix(QQR.vol_Netherlands)
persp3D(x,y,z, theta=320, phi=20, shade=0,d =5, r=9, bty="b2", border ="black", lwd=0.3, lty=1,axes = TRUE, facets = TRUE,
        colkey=list(side = 4, plot = TRUE, length = 1, width = 0.6, dist = 0.03, shift = 0, addlines = FALSE, col.clab = NULL, cex.clab = par("cex.lab"), side.clab = NULL, line.clab = NULL, adj.clab = NULL, font.clab = NULL),
        ticktype= "detailed", xlab = "quantiles of EPU", ylab = "quantiles of volatility", zlab= "beta",main="Netherlands")


#Causality in quantiles (Example here is GERMANY, apply the same for other countries)------------------------------------
data<-read_excel("data_paper3.xlsx", sheet ="data")
data2 <-read_excel("data_paper3.xlsx", sheet ="data2")
data3<-read_excel("data_paper3.xlsx", sheet ="volatility")
data4 <-read_excel("data_paper3.xlsx", sheet ="volatility2")
library(quantreg) # For quantile regression
library(qpcR) # for the Hannan Quinn Criterion
library(sandwich) #for HAC Estimator
# Set up the orders of lagged dependent variables and lagged causal
#####################################################################
px = 0 #order of the lagged independenr variable has to be specified
py = 1 #order of the lagged dependent variable has to be specified

build_arp = function(x,maxorder){ # maxorder is the maximum lag number
  ordmax = maxorder  
  lagmat = NULL
  l1 = length(x)
  for (i in 1:ordmax){
    lagmat = cbind(lagmat[-1,],x[(1):(l1-i)]) # lagged matrix
  }
  return(lagmat)
}

J <- px ## Number of linear hypothesis in the null under testing
R1 <- matrix(0:0, J, py+1) ## py+1 is the number of explanatory variables,
R2 <- diag(px)
R <- cbind(R1, R2)
r <- matrix(0:0,px,1)

#Causality in quantile - (Bitcoin returns ~ Germnay EPU)---------------------------------------------------
close <- data$BTC
n <- length(close)
ret <- 100*(log(close[2:n]) - log(close[1:n-1]))
nret <- length(ret)
EPU <- log(data$Germany[2:n])
t <- seq(1, nret, 1)
t2 <- t^2
fit <- lm(EPU ~ t + t2)
dEPU <- fit$residuals

# Quantile Regression Non-Granger Causlaity Test at all taus
tau = 5:95/100
ntaus <- length(tau)
Wtau1<-rep(0,ntaus)
for (i in 1:ntaus){
  fit <- rq(y~ lagy + lagx,tau[i])
  getcov<-summary.rq(fit, se="nid", covariance=TRUE, hs = FALSE)
  covm<-(getcov$cov)
  betahat<-coef(fit)
  Rbetahat <- R %*% betahat
  Wtau1[i] <- t(Rbetahat - r) %*% solve(R %*% covm %*% t(R)) %*% (Rbetahat)
}
print (max(Wtau1))
(length(Wtau1))
Wtau1
print ("---------Output from the Max Wald test:---------")
print ("the quantile of interest is:")
print ("the maximum value of the statistic is:")
print (max(Wtau1))

cau_Germany <- Wtau1
#Graph
plot.new()
rect(par("usr")[1], par("usr")[3],
     par("usr")[2], par("usr")[4],
     col = "#f7f7f7") # Color
par(new = TRUE)
plot(tau,cau_Germany, type = "l",lty=1, lwd=2, col="blue",xlab = "", ylab = "",
     mar = c(5.1, 4.1, 2.1, 2.1))
title("Germany", line = -1.2)
abline(h=7.19,lwd=2, lty=2, col="red")


#Bitcoin volatility ~ Germany EPU------------------------------------------------
n <- length(data3$volatility)
vola <- data3$volatility[2:n]
nret <- length(vola)
EPU <- log(data3$Germany[2:n])
t <- seq(1, nret, 1)
t2 <- t^2
fit <- lm(EPU ~ t + t2)
dEPU <- fit$residuals

#Quantile Regression Non-Granger Causlaity Test at all taus
tau = 5:95/100
ntaus <- length(tau)
Wtau1<-rep(0,ntaus)
for (i in 1:ntaus){
  fit <- rq(y~ lagy + lagx,tau[i])
  getcov<-summary.rq(fit, se="nid", covariance=TRUE, hs = FALSE)
  covm<-(getcov$cov)
  betahat<-coef(fit)
  Rbetahat <- R %*% betahat
  Wtau1[i] <- t(Rbetahat - r) %*% solve(R %*% covm %*% t(R)) %*% (Rbetahat)
}
print (max(Wtau1))
(length(Wtau1))
Wtau1
print ("---------Output from the Max Wald test:---------")
print ("the quantile of interest is:")
print ("the maximum value of the statistic is:")
print (max(Wtau1))
#print("--------------------------------------------")
cau_vol_Germany <-Wtau1

#Graph
plot.new()
rect(par("usr")[1], par("usr")[3],
     par("usr")[2], par("usr")[4],
     col = "#f7f7f7") # Color
par(new = TRUE)
plot(tau,cau_vol_Germany, type = "l",lty=1, lwd=2, col="blue",xlab = "", ylab = "",
     mar = c(5.1, 4.1, 2.1, 2.1))
title("Germany", line = -1.2)
abline(h=7.19,lwd=2, lty=2, col="red")

#DCC-GARCH estimation--------------------------------
library(tseries)
library(rugarch)
library(rmgarch)
library(xts)

dcc<-read_excel("data_paper3.xlsx", sheet ="robustness")
NCSKEW <- xts(dcc$NCSKEW,dcc$Date)
Canada <- xts(dcc$Canada,dcc$Date)
France <- xts(dcc$France,dcc$Date)
Germany <- xts(dcc$Germany,dcc$Date)
Japan <- xts(dcc$Japan,dcc$Date)
Russia <- xts(dcc$Russia,dcc$Date)
US <- xts(dcc$US,dcc$Date)
UK <- xts(dcc$UK,dcc$Date)
China <- xts(dcc$China,dcc$Date)
Netherland <- xts(dcc$Netherlands,dcc$Date)
Singapore <- xts(dcc$Singapore,dcc$Date)

model1=ugarchspec(mean.model=list(armaOrder=c(0,0)),variance.model=list(garchOrder=c(1,1),model="sGARCH"),distribution.model="norm")
modelspec=dccspec(uspec = multispec(replicate(2,model1)),dccOrder=c(1,1),distribution = "mvnorm")                  
modelfit=dccfit(modelspec,data=data.frame(NCSKEW,US))
plot(modelfit)
modelfit1=dccfit(modelspec,data=data.frame(NCSKEW,UK))
plot(modelfit1)
modelfit2=dccfit(modelspec,data=data.frame(NCSKEW,Canada))
plot(modelfit2)
modelfit3=dccfit(modelspec,data=data.frame(NCSKEW,France))
plot(modelfit3)
modelfit4=dccfit(modelspec,data=data.frame(NCSKEW,Germany))
plot(modelfit4)
modelfit5=dccfit(modelspec,data=data.frame(NCSKEW,Japan))
plot(modelfit5)
modelfit6=dccfit(modelspec,data=data.frame(NCSKEW,Russia))
plot(modelfit6)
modelfit7=dccfit(modelspec,data=data.frame(NCSKEW,China))
plot(modelfit7)
modelfit8=dccfit(modelspec,data=data.frame(NCSKEW,Netherland))
plot(modelfit8)
modelfit9=dccfit(modelspec,data=data.frame(NCSKEW,Singapore))
plot(modelfit9)