library(readxl)
library(quantreg)
library(DataCombine)
library(tidyr)
library(quantreg)
library(pastecs)

#1. Import data------------------------------------------------------------------
spot <- read_excel("data_paper2.xlsx", sheet = 1)
forecast3<-read_excel("data_paper2.xlsx",sheet=2)
forecast6<-read_excel("data_paper2.xlsx",sheet=3)

ir <- read_excel("data_paper2.xlsx", sheet =4)
EPU <- read_excel("data_paper2.xlsx", sheet =5)                 

ir[ir$ir_AUD3m=="NA","ir_AUD3m"]<-NA
ir$ir_AUD3m <-as.numeric(ir$ir_AUD3m)
ir[ir$ir_EUR3m=="NA","ir_EUR3m"]<-NA
ir$ir_EUR3m <-as.numeric(ir$ir_EUR3m)

ir[ir$ir_AUD6m=="NA","ir_AUD6m"]<-NA
ir$ir_AUD6m <-as.numeric(ir$ir_AUD6m)
ir[ir$ir_EUR6m=="NA","ir_EUR6m"]<-NA
ir$ir_EUR6m <-as.numeric(ir$ir_EUR6m)


ir$fw_GBPUSD3m <-(log(1+ir$ir_GBP3m/(100*365/90))-log(1+ir$ir_USD3m/(100*365/90)))*100
ir$fw_EURUSD3m <-(log(1+ir$ir_EUR3m/(100*365/90))-log(1+ir$ir_USD3m/(100*365/90)))*100
ir$fw_CHFUSD3m <-(log(1+ir$ir_CHF3m/(100*365/90))-log(1+ir$ir_USD3m/(100*365/90)))*100
ir$fw_AUDUSD3m <-(log(1+ir$ir_AUD3m/(100*365/90))-log(1+ir$ir_USD3m/(100*365/90)))*100
ir$fw_JPYUSD3m <-(log(1+ir$ir_JPY3m/(100*365/90))-log(1+ir$ir_USD3m/(100*365/90)))*100
ir$fw_CADUSD3m <-(log(1+ir$ir_CAD3m/(100*365/90))-log(1+ir$ir_USD3m/(100*365/90)))*100

ir$fw_GBPUSD6m <-(log(1+ir$ir_GBP6m/(100*365/180))-log(1+ir$ir_USD6m/(100*365/180)))*100
ir$fw_EURUSD6m <-(log(1+ir$ir_EUR6m/(100*365/180))-log(1+ir$ir_USD6m/(100*365/180)))*100
ir$fw_CHFUSD6m <-(log(1+ir$ir_CHF6m/(100*365/180))-log(1+ir$ir_USD6m/(100*365/180)))*100
ir$fw_AUDUSD6m <-(log(1+ir$ir_AUD6m/(100*365/180))-log(1+ir$ir_USD6m/(100*365/180)))*100
ir$fw_JPYUSD6m <-(log(1+ir$ir_JPY6m/(100*365/180))-log(1+ir$ir_USD6m/(100*365/180)))*100
ir$fw_CADUSD6m <-(log(1+ir$ir_CAD6m/(100*365/180))-log(1+ir$ir_USD6m/(100*365/180)))*100

df1<-data.frame(index=spot$date,spot_GBPUSD=log(spot$GBPUSD),forecast_GBPUSD3m=log(forecast3$GBPUSD3m),
               spot_EURUSD=log(spot$EURUSD),forecast_EURUSD3m=log(forecast3$EURUSD3m),
               spot_CHFUSD=log(spot$CHFUSD),forecast_CHFUSD3m=log(forecast3$CHFUSD3m),
               spot_AUDUSD=log(spot$AUDUSD),forecast_AUDUSD3m=log(forecast3$AUDUSD3m),
               spot_JPYUSD=log(spot$JPYUSD),forecast_JPYUSD3m=log(forecast3$JPYUSD3m),
               spot_CADUSD=log(spot$CADUSD),forecast_CADUSD3m=log(forecast3$CADUSD3m),
               fw_GBPUSD3m=ir$fw_GBPUSD3m,
               fw_EURUSD3m=ir$fw_EURUSD3m,
               fw_CHFUSD3m=ir$fw_CHFUSD3m,
               fw_AUDUSD3m=ir$fw_AUDUSD3m,
               fw_JPYUSD3m=ir$fw_JPYUSD3m,
               fw_CADUSD3m=ir$fw_CADUSD3m,
               US_EPU=EPU$US_monetary_uncertainty)

df2<-data.frame(index=spot$date,spot_GBPUSD=log(spot$GBPUSD),forecast_GBPUSD6m=log(forecast6$GBPUSD6m),
               spot_EURUSD=log(spot$EURUSD),forecast_EURUSD6m=log(forecast6$EURUSD6m),
               spot_CHFUSD=log(spot$CHFUSD),forecast_CHFUSD6m=log(forecast6$CHFUSD6m),
               spot_AUDUSD=log(spot$AUDUSD),forecast_AUDUSD6m=log(forecast6$AUDUSD6m),
               spot_JPYUSD=log(spot$JPYUSD),forecast_JPYUSD6m=log(forecast6$JPYUSD6m),
               spot_CADUSD=log(spot$CADUSD),forecast_CADUSD6m=log(forecast6$CADUSD6m),
               fw_GBPUSD6m=ir$fw_GBPUSD6m,
               fw_EURUSD6m=ir$fw_EURUSD6m,
               fw_CHFUSD6m=ir$fw_CHFUSD6m,
               fw_AUDUSD6m=ir$fw_AUDUSD6m,
               fw_JPYUSD6m=ir$fw_JPYUSD6m,
               fw_CADUSD6m=ir$fw_CADUSD6m,
               US_EPU=EPU$US_monetary_uncertainty)

#Create 3-month forecast return series
df1$GBPUSD_return<-(lead(df1$forecast_GBPUSD3m,3)-df1$spot_GBPUSD)*100
df1$EURUSD_return<-(lead(df1$forecast_EURUSD3m,3)-df1$spot_EURUSD)*100
df1$CHFUSD_return<-(lead(df1$forecast_CHFUSD3m,3)-df1$spot_CHFUSD)*100
df1$AUDUSD_return<-(lead(df1$forecast_AUDUSD3m,3)-df1$spot_AUDUSD)*100
df1$JPYUSD_return<-(lead(df1$forecast_JPYUSD3m,3)-df1$spot_JPYUSD)*100
df1$CADUSD_return<-(lead(df1$forecast_CADUSD3m,3)-df1$spot_CADUSD)*100
#Create 6-month forecast return series
df2$GBPUSD_return<-(lead(df2$forecast_GBPUSD6m,6)-df2$spot_GBPUSD)*100
df2$EURUSD_return<-(lead(df2$forecast_EURUSD6m,6)-df2$spot_EURUSD)*100
df2$CHFUSD_return<-(lead(df2$forecast_CHFUSD6m,6)-df2$spot_CHFUSD)*100
df2$AUDUSD_return<-(lead(df2$forecast_AUDUSD6m,6)-df2$spot_AUDUSD)*100
df2$CNYUSD_return<-(lead(df2$forecast_CNYUSD6m,6)-df2$spot_CNYUSD)*100
df2$JPYUSD_return<-(lead(df2$forecast_JPYUSD6m,6)-df2$spot_JPYUSD)*100
df2$CADUSD_return<-(lead(df2$forecast_CADUSD6m,6)-df2$spot_CADUSD)*100

#2. Group database for each pair of currencies----------------------------------------------------------
#3-month horizon
#GBPUSD
df1_GBP<-data.frame(index=df1$date,GBPUSD_return=df1$GBPUSD_return,fw_GBPUSD3m=df1$fw_GBPUSD3m,EPU=df1$US_EPU)
df1_GBP<-df1_GBP[complete.cases(df1_GBP),]
y <- as.matrix(df1_GBP$GBPUSD_return)
q <- as.matrix(df1_GBP$EPU)
x <- as.matrix(df1_GBP$fw_GBPUSD3m)
#EURUSD
df1_EUR<-data.frame(index=df1$date,EURUSD_return=df1$EURUSD_return,fw_EURUSD3m=df1$fw_EURUSD3m,EPU=df1$US_EPU)
df1_EUR<-df1_EUR[complete.cases(df1_EUR),]
y <- as.matrix(df1_EUR$EURUSD_return)
q <- as.matrix(df1_EUR$EPU)
x <- as.matrix(df1_EUR$fw_EURUSD3m)
#CHFUSD
df1_CHF<-data.frame(index=df1$date,CHFUSD_return=df1$CHFUSD_return,fw_CHFUSD3m=df1$fw_CHFUSD3m,EPU=df1$US_EPU)
df1_CHF<-df1_CHF[complete.cases(df1_CHF),]
y <- as.matrix(df1_CHF$CHFUSD_return)
q <- as.matrix(df1_CHF$EPU)
x <- as.matrix(df1_CHF$fw_CHFUSD3m)
#AUDUSD
df1_AUD<-data.frame(index=df1$date,AUDUSD_return=df1$AUDUSD_return,fw_AUDUSD3m=df1$fw_AUDUSD3m,EPU=df1$US_EPU)
df1_AUD<-df1_AUD[complete.cases(df1_AUD),]
y <- as.matrix(df1_AUD$AUDUSD_return)
q <- as.matrix(df1_AUD$EPU)
x <- as.matrix(df1_AUD$fw_AUDUSD3m)

#JPYUSD
df1_JPY<-data.frame(index=df1$date,JPYUSD_return=df1$JPYUSD_return,fw_JPYUSD3m=df1$fw_JPYUSD3m,EPU=df1$US_EPU)
df1_JPY<-df1_JPY[complete.cases(df1_JPY),]
y <- as.matrix(df1_JPY$JPYUSD_return)
q <- as.matrix(df1_JPY$EPU)
x <- as.matrix(df1_JPY$fw_JPYUSD3m)

#CADUSD
df1_CAD<-data.frame(index=df1$index,CADUSD_return=df1$CADUSD_return,fw_CADUSD3m=df1$fw_CADUSD3m,EPU=df1$US_EPU)
df1_CAD<-df1_CAD[complete.cases(df1_CAD),]
y <- as.matrix(df1_CAD$CADUSD_return)
q <- as.matrix(df1_CAD$EPU)
x <- as.matrix(df1_CAD$fw_CADUSD3m)

#6-month horizon
#GBPUSD
df2_GBP<-data.frame(index=df2$date,GBPUSD_return=df2$GBPUSD_return,fw_GBPUSD6m=df2$fw_GBPUSD6m,EPU=df2$US_EPU)
df2_GBP<-df2_GBP[complete.cases(df2_GBP),]
y <- as.matrix(df2_GBP$GBPUSD_return)
q <- as.matrix(df2_GBP$EPU)
x <- as.matrix(df2_GBP$fw_GBPUSD6m)

#EURUSD
df2_EUR<-data.frame(index=df2$date,EURUSD_return=df2$EURUSD_return,fw_EURUSD6m=df2$fw_EURUSD6m,EPU=df2$US_EPU)
df2_EUR<-df2_EUR[complete.cases(df2_EUR),]
y <- as.matrix(df2_EUR$EURUSD_return)
q <- as.matrix(df2_EUR$EPU)
x <- as.matrix(df2_EUR$fw_EURUSD6m)

#CHFUSD
df2_CHF<-data.frame(index=df2$date,CHFUSD_return=df2$CHFUSD_return,fw_CHFUSD6m=df2$fw_CHFUSD6m,EPU=df2$US_EPU)
df2_CHF<-df2_CHF[complete.cases(df2_CHF),]
y <- as.matrix(df2_CHF$CHFUSD_return)
q <- as.matrix(df2_CHF$EPU)
x <- as.matrix(df2_CHF$fw_CHFUSD6m)

#AUDUSD
df2_AUD<-data.frame(index=df2$date,AUDUSD_return=df2$AUDUSD_return,fw_AUDUSD6m=df2$fw_AUDUSD6m,EPU=df2$US_EPU)
y <- as.matrix(df2_AUD$AUDUSD_return)
q <- as.matrix(df2_AUD$EPU)
x <- as.matrix(df2_AUD$fw_AUDUSD6m)

#JPYUSD
df2_JPY<-data.frame(index=df2$date,JPYUSD_return=df2$JPYUSD_return,fw_JPYUSD6m=df2$fw_JPYUSD6m,EPU=df2$US_EPU)
df2_JPY<-df2_JPY[complete.cases(df2_JPY),]
y <- as.matrix(df2_JPY$JPYUSD_return)
q <- as.matrix(df2_JPY$EPU)
x <- as.matrix(df2_JPY$fw_JPYUSD6m)

#CADUSD
df2_CAD<-data.frame(index=df2$date,CADUSD_return=df2$CADUSD_return,fw_CADUSD6m=df2$fw_CADUSD6m,EPU=df2$US_EPU)
df2_CAD<-df2_CAD[complete.cases(df2_CAD),]
y <- as.matrix(df2_CAD$CADUSD_return)
q <- as.matrix(df2_CAD$US_EPU2)
x <- as.matrix(df2_CAD$fw_CADUSD6m)

#3. Graphical description of 3-month expected exchange rate changes and 3-month interest rate differential--------------------
plot1_UK <- ggplot(df1_GBP, aes(x = index)) +
  geom_line(aes(y = GBPUSD_return), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("GBP/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot2_UK <- ggplot(df1_GBP, aes(x = index)) +
  geom_line(aes(y = fw_GBPUSD3m), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("GBP/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot1_EUR <- ggplot(df1_EUR, aes(x = index)) +
  geom_line(aes(y = EURUSD_return), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("EUR/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot2_EUR <- ggplot(df1_EUR, aes(x = index)) +
  geom_line(aes(y = fw_EURUSD3m), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("EUR/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot1_CHF <- ggplot(df1_CHF, aes(x = index)) +
  geom_line(aes(y = CHFUSD_return), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("CHF/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot2_CHF <- ggplot(df1_CHF, aes(x = index)) +
  geom_line(aes(y = fw_CHFUSD3m), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("CHF/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot1_AUD <- ggplot(df1_AUD, aes(x = index)) +
  geom_line(aes(y = AUDUSD_return), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("AUD/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot2_AUD <- ggplot(df1_AUD, aes(x = index)) +
  geom_line(aes(y = fw_AUDUSD3m), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("AUD/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot1_JPY <- ggplot(df1_JPY, aes(x = index)) +
  geom_line(aes(y = JPYUSD_return), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("JPY/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot2_JPY <- ggplot(df1_JPY, aes(x = index)) +
  geom_line(aes(y = fw_JPYUSD3m), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("JPY/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot1_CAD <- ggplot(df1_CAD, aes(x = index)) +
  geom_line(aes(y = CADUSD_return), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("CAD/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))
plot2_CAD <- ggplot(df1_CAD, aes(x = index)) +
  geom_line(aes(y = fw_CADUSD3m), colour = "blue", linetype = 1, size = 0.2) +
  ylab("") + xlab("time") + ggtitle("CAD/USD") + theme_bw()+theme(plot.title=element_text(size=10,hjust=0.5))

#4. Summary statistics--------------------------------------------------------------------------------------
#3-month exchange rate changes and forward premia
stat.desc(df1_GBP$GBPUSD_return,df1_EUR$EURUSD_return,df1_CHF$CHFUSD_return,df1_AUD$AUDUSD_return,df1_JPY$JPYUSD_return,df1_CAD$CADUSD_return)
mean_sta<-sapply(df1[,c("fw_GBPUSD3m","fw_EURUSD3m","fw_CHFUSD3m","fw_AUDUSD3m","fw_JPYUSD3m","fw_CADUSD3m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")], mean, na.rm=TRUE)*4
std_sta<-sapply(df1[,c("fw_GBPUSD3m","fw_EURUSD3m","fw_CHFUSD3m","fw_AUDUSD3m","fw_JPYUSD3m","fw_CADUSD3m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")],sd,na.rm=TRUE)*sqrt(4)
max_sta<-sapply(df1[,c("fw_GBPUSD3m","fw_EURUSD3m","fw_CHFUSD3m","fw_AUDUSD3m","fw_JPYUSD3m","fw_CADUSD3m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")],max,na.rm=TRUE)*2
min_sta<-sapply(df1[,c("fw_GBPUSD3m","fw_EURUSD3m","fw_CHFUSD3m","fw_AUDUSD3m","fw_JPYUSD3m","fw_CADUSD3m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")],min,na.rm=TRUE)*2
sta3m<-data.frame(mean_sta,std_sta,max_sta,min_sta)
View(sta3m)
#6-month exchange rate changes and forward premia
stat.desc(df2_GBP$GBPUSD_return,df2_EUR$EURUSD_return,df2_CHF$CHFUSD_return,df2_AUD$AUDUSD_return,df2_JPY$JPYUSD_return,df2_CAD$CADUSD_return)
mean_sta<-sapply(df2[,c("fw_GBPUSD6m","fw_EURUSD6m","fw_CHFUSD6m","fw_AUDUSD6m","fw_JPYUSD6m","fw_CADUSD6m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")], mean, na.rm=TRUE)*2
std_sta<-sapply(df2[,c("fw_GBPUSD6m","fw_EURUSD6m","fw_CHFUSD6m","fw_AUDUSD6m","fw_JPYUSD6m","fw_CADUSD6m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")],sd,na.rm=TRUE)*sqrt(2)
max_sta<-sapply(df2[,c("fw_GBPUSD6m","fw_EURUSD6m","fw_CHFUSD6m","fw_AUDUSD6m","fw_JPYUSD6m","fw_CADUSD6m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")],max,na.rm=TRUE)*2
min_sta<-sapply(df2[,c("fw_GBPUSD6m","fw_EURUSD6m","fw_CHFUSD6m","fw_AUDUSD6m","fw_JPYUSD6m","fw_CADUSD6m","GBPUSD_return","EURUSD_return","CHFUSD_return","AUDUSD_return","JPYUSD_return","CADUSD_return")],min,na.rm=TRUE)*2
sta6m<-data.frame(mean_sta,std_sta,max_sta,min_sta)
View(sta6m)

#5. Quantile Threshold Estimation ----------------------------------------------------
joint_thresh <- function(y,x,q,graph){
  n=nrow(y)
  k=ncol(x)
  e=y-x%*%rq(y~x,tau)$coefficients[2]-rq(y~x,tau)$coefficients[1]
  s0 <- det(t(e)%*%e)    
  n1 <- round(.05*n)+k
  n2 <- round(.95*n)-k
  qs <- sort(q)
  qs <- qs[n1:n2]
  qs <- as.matrix(unique(qs))
  qn <- nrow(qs)
  sn <- matrix(0,qn,1)
  for (r in 1:qn){
    d <- (q<=qs[r])
    xx <- (x)*(d%*%matrix(1,1,k))
    xx <- xx-x%*%rq(xx~x,tau)$coefficients[2]-rq(xx~x,tau)$coefficients[1]
    ex <- e-xx%*%rq(e~xx,tau)$coefficients[2]-rq(e~xx,tau)$coefficients[1]
    exw <- ex*(tau-(ex<0))
    sn[r] <- sum(exw)   
  }
  r <- which.min(sn)
  smin <- sn[r] 
  qhat <- qs[r]
  d <- (q<=qhat)
  x1 <- x*(d%*%matrix(1,1,k))
  x2 <- x*((1-d)%*%matrix(1,1,k))
  beta1 <- rq(y~x1,tau)$coefficients[2]
  beta2 <- rq(y~x2,tau)$coefficients[2]
  yhat <- x1%*%beta1+x2%*%beta2
  result1 <- summary(rq(y~x1,tau))$coefficients
  result2 <- summary(rq(y~x2,tau))$coefficients
  result1_full<-summary(rq(y~x1,tau),se="iid")$coefficients
  result2_full<-summary(rq(y~x2,tau),se="iid")$coefficients
  lr<-(sn-smin)/(tau*(1-tau))
  if(tau==0.1&tau==0.9){c0<-12.254
  } else if (tau==0.2&tau==0.8){
    c0 <-9.19
  }else if (tau==0.3&tau==0.7){
    c0 <- 8.022
  }else if (tau==0.4&tau==0.6){
    c0<-7.504
  }else{
    c0<-7.352
  }
  clr <- matrix(1,qn,1)*c0
  if (graph==1){
  xxlim <- range(qs)          
  yylim <- range(rbind(lr,c0))
  mtit <- paste(c("Confidence Interval Construction \nfor Threshold at tau="),tau,sep="")
  ytit <- "Likelihood Ratio Sequence in gamma" 
  xtit <- paste(c("Threshold Variable"),c(" "),c("(tau="),tau,c(")"),sep="")
  plot(qs,lr,lty=1,col=1,xlim=xxlim,ylim=yylim,type="l",ann=0)
  lines(qs,clr,lty=2,col=2)
  title(main=mtit,ylab=ytit,xlab=xtit)
  tit1 <- "LRn(gamma)"
  tit2 <- "95% Critical"
  legend("topright",c(tit1,tit2),lty=c(1,2),col=c(1,2))
  }
  list(yhat=yhat,qhat=qhat,regime1=result1,regime2=result2,result1_full=result1_full,result2_full=result2_full)
}

#Results
#tau=0.1
tau=0.1
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.1<-out1$regime1
regime2_UK_0.1<-out1$regime2
regime1_UK_full_0.1<-out1$result1_full
regime2_UK_full_0.1<-out1$result2_full
qhat_UK_0.1<-out1$qhat
#tau=0.2
tau=0.2
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.2<-out1$regime1
regime2_UK_0.2<-out1$regime2
regime1_UK_full_0.2<-out1$result1_full
regime2_UK_full_0.2<-out1$result2_full
qhat_UK_0.2<-out1$qhat
#tau=0.3
tau=0.3
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.3<-out1$regime1
regime2_UK_0.3<-out1$regime2
regime1_UK_full_0.3<-out1$result1_full
regime2_UK_full_0.3<-out1$result2_full
qhat_UK_0.3<-out1$qhat
#tau0.4
tau=0.4
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.4<-out1$regime1
regime2_UK_0.4<-out1$regime2
regime1_UK_full_0.4<-out1$result1_full
regime2_UK_full_0.4<-out1$result2_full
qhat_UK_0.4<-out1$qhat
#tau=0.5
tau=0.5
out1 <- joint_thresh(y=y,x=x,q=q,1)
qhat_UK_0.5<-out1$qhat
regime1_UK_0.5<-out1$regime1
regime2_UK_0.5<-out1$regime2
regime1_UK_full_0.5<-out1$result1_full
regime2_UK_full_0.5<-out1$result2_full
#tau=0.6
tau=0.6
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.6<-out1$regime1
regime2_UK_0.6<-out1$regime2
regime1_UK_full_0.6<-out1$result1_full
regime2_UK_full_0.6<-out1$result2_full
qhat_UK_0.6<-out1$qhat
#tau=0.7
tau=0.7
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.7<-out1$regime1
regime2_UK_0.7<-out1$regime2
regime1_UK_full_0.7<-out1$result1_full
regime2_UK_full_0.7<-out1$result2_full
qhat_UK_0.7<-out1$qhat
#tau=0.8
tau=0.8
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.8<-out1$regime1
regime2_UK_0.8<-out1$regime2
regime1_UK_full_0.8<-out1$result1_full
regime2_UK_full_0.8<-out1$result2_full
qhat_UK_0.8<-out1$qhat
#tau=0.9
tau=0.9
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_UK_0.9<-out1$regime1
regime2_UK_0.9<-out1$regime2
regime1_UK_full_0.9<-out1$result1_full
regime2_UK_full_0.9<-out1$result2_full
qhat_UK_0.9<-out1$qhat
#Estimation results
#GBPUSD
regime1_UK_sta<-rbind(t(regime1_UK_full_0.1),t(regime1_UK_full_0.2),t(regime1_UK_full_0.3),t(regime1_UK_full_0.4),
                      t(regime1_UK_full_0.5),t(regime1_UK_full_0.6),t(regime1_UK_full_0.7),t(regime1_UK_full_0.8),
                      t(regime1_UK_full_0.9))
regime2_UK_sta<-rbind(t(regime2_UK_full_0.1),t(regime2_UK_full_0.2),t(regime2_UK_full_0.3),t(regime2_UK_full_0.4),
                      t(regime2_UK_full_0.5),t(regime2_UK_full_0.6),t(regime2_UK_full_0.7),t(regime2_UK_full_0.8),
                      t(regime2_UK_full_0.9))
UK_sta<-cbind(regime1_UK_sta,regime2_UK_sta)
#Graph of results
TAU<-seq(0.1,0.9,0.1)
qhat_UK<-c(qhat_UK_0.1,qhat_UK_0.2,qhat_UK_0.3,qhat_UK_0.4,qhat_UK_0.5,qhat_UK_0.6,qhat_UK_0.7,qhat_UK_0.8,qhat_UK_0.9)
bangthamchieu_UK<-cbind(tau=TAU,qhat=qhat_UK)
heso_regime1_UK<-rbind(regime1_UK_0.1[2,],regime1_UK_0.2[2,],regime1_UK_0.3[2,],regime1_UK_0.4[2,],
                       regime1_UK_0.5[2,],regime1_UK_0.6[2,],regime1_UK_0.7[2,],regime1_UK_0.8[2,],regime1_UK_0.9[2,])
heso_regime1_UK<-cbind(TAU,heso_regime1_UK)
heso_regime2_UK<-rbind(regime2_UK_0.1[2,],regime2_UK_0.2[2,],regime2_UK_0.3[2,],regime2_UK_0.4[2,],
                       regime2_UK_0.5[2,],regime2_UK_0.6[2,],regime2_UK_0.7[2,],regime2_UK_0.8[2,],regime2_UK_0.9[2,])
heso_regime2_UK<-cbind(TAU,heso_regime2_UK)
#1st regime
plot(rep(TAU,2),c(heso_regime1_UK[,3],heso_regime1_UK[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime1_UK[,3], rev(heso_regime1_UK[,4])), col = "LightSkyBlue")
points(TAU, heso_regime1_UK[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime1_UK[,2], col = "blue")
lr_UK<-summary(lm(y~x))$coefficients[2,1]
lr_set_UK<-summary(lm(y~x))$coefficients[2,2]
abline(h=lr_UK,col="red")
abline(h=lr_UK-1.96*lr_set_UK,col="red",lty="dotted")
abline(h=lr_UK+1.96*lr_set_UK,col="red",lty="dotted")
title(bquote("UK_"~alpha*2))
#2nd regime
plot(rep(TAU,2),c(heso_regime2_UK[,3],heso_regime2_UK[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime2_UK[,3], rev(heso_regime2_UK[,4])), col = "LightSkyBlue")
points(TAU, heso_regime2_UK[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime2_UK[,2], col = "blue")
abline(h=lr_UK,col="red")
abline(h=lr_UK-1.96*lr_set_UK,col="red",lty="dotted")
abline(h=lr_UK+1.96*lr_set_UK,col="red",lty="dotted")
title(bquote("UK_"~beta*2))

#JPYUSD
#tau=0.1
tau=0.1
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.1<-out1$regime1
regime2_JPY_0.1<-out1$regime2
regime1_JPY_full_0.1<-out1$result1_full
regime2_JPY_full_0.1<-out1$result2_full
qhat_JPY_0.1<-out1$qhat
#tau=0.2
tau=0.2
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.2<-out1$regime1
regime2_JPY_0.2<-out1$regime2
regime1_JPY_full_0.2<-out1$result1_full
regime2_JPY_full_0.2<-out1$result2_full
qhat_JPY_0.2<-out1$qhat
#tau=0.3
tau=0.3
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.3<-out1$regime1
regime2_JPY_0.3<-out1$regime2
regime1_JPY_full_0.3<-out1$result1_full
regime2_JPY_full_0.3<-out1$result2_full
qhat_JPY_0.3<-out1$qhat
#tau0.4
tau=0.4
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.4<-out1$regime1
regime2_JPY_0.4<-out1$regime2
regime1_JPY_full_0.4<-out1$result1_full
regime2_JPY_full_0.4<-out1$result2_full
qhat_JPY_0.4<-out1$qhat
#tau=0.5
tau=0.5
out1 <- joint_thresh(y=y,x=x,q=q,1)
qhat_JPY_0.5<-out1$qhat
regime1_JPY_0.5<-out1$regime1
regime2_JPY_0.5<-out1$regime2
regime1_JPY_full_0.5<-out1$result1_full
regime2_JPY_full_0.5<-out1$result2_full
#tau=0.6
tau=0.6
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.6<-out1$regime1
regime2_JPY_0.6<-out1$regime2
regime1_JPY_full_0.6<-out1$result1_full
regime2_JPY_full_0.6<-out1$result2_full
qhat_JPY_0.6<-out1$qhat
#tau=0.7
tau=0.7
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.7<-out1$regime1
regime2_JPY_0.7<-out1$regime2
regime1_JPY_full_0.7<-out1$result1_full
regime2_JPY_full_0.7<-out1$result2_full
qhat_JPY_0.7<-out1$qhat
#tau=0.8
tau=0.8
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.8<-out1$regime1
regime2_JPY_0.8<-out1$regime2
regime1_JPY_full_0.8<-out1$result1_full
regime2_JPY_full_0.8<-out1$result2_full
qhat_JPY_0.8<-out1$qhat
#tau=0.9
tau=0.9
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_JPY_0.9<-out1$regime1
regime2_JPY_0.9<-out1$regime2
regime1_JPY_full_0.9<-out1$result1_full
regime2_JPY_full_0.9<-out1$result2_full
qhat_JPY_0.9<-out1$qhat
#Estimation results
regime1_JPY_sta<-rbind(t(regime1_JPY_full_0.1),t(regime1_JPY_full_0.2),t(regime1_JPY_full_0.3),t(regime1_JPY_full_0.4),
                       t(regime1_JPY_full_0.5),t(regime1_JPY_full_0.6),t(regime1_JPY_full_0.7),t(regime1_JPY_full_0.8),
                       t(regime1_JPY_full_0.9))
regime2_JPY_sta<-rbind(t(regime2_JPY_full_0.1),t(regime2_JPY_full_0.2),t(regime2_JPY_full_0.3),t(regime2_JPY_full_0.4),
                       t(regime2_JPY_full_0.5),t(regime2_JPY_full_0.6),t(regime2_JPY_full_0.7),t(regime2_JPY_full_0.8),
                       t(regime2_JPY_full_0.9))
JPY_sta<-cbind(regime1_JPY_sta,regime2_JPY_sta)
#Graph of results
TAU<-seq(0.1,0.9,0.1)
qhat_JPY<-c(qhat_JPY_0.1,qhat_JPY_0.2,qhat_JPY_0.3,qhat_JPY_0.4,qhat_JPY_0.5,qhat_JPY_0.6,qhat_JPY_0.7,qhat_JPY_0.8,qhat_JPY_0.9)
bangthamchieu_JPY<-cbind(tau=TAU,qhat=qhat_JPY)
heso_regime1_JPY<-rbind(regime1_JPY_0.1[2,],regime1_JPY_0.2[2,],regime1_JPY_0.3[2,],regime1_JPY_0.4[2,],
                        regime1_JPY_0.5[2,],regime1_JPY_0.6[2,],regime1_JPY_0.7[2,],regime1_JPY_0.8[2,],regime1_JPY_0.9[2,])
heso_regime1_JPY<-cbind(TAU,heso_regime1_JPY)
heso_regime2_JPY<-rbind(regime2_JPY_0.1[2,],regime2_JPY_0.2[2,],regime2_JPY_0.3[2,],regime2_JPY_0.4[2,],
                        regime2_JPY_0.5[2,],regime2_JPY_0.6[2,],regime2_JPY_0.7[2,],regime2_JPY_0.8[2,],regime2_JPY_0.9[2,])
heso_regime2_JPY<-cbind(TAU,heso_regime2_JPY)
#1st regime
plot(rep(TAU,2),c(heso_regime1_JPY[,3],heso_regime1_JPY[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime1_JPY[,3], rev(heso_regime1_JPY[,4])), col = "LightSkyBlue")
points(TAU, heso_regime1_JPY[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime1_JPY[,2], col = "blue")
lr_JPY<-summary(lm(y~x))$coefficients[2,1]
lr_set_JPY<-summary(lm(y~x))$coefficients[2,2]
abline(h=lr_JPY,col="red")
abline(h=lr_JPY-1.96*lr_set_JPY,col="red",lty="dotted")
abline(h=lr_JPY+1.96*lr_set_JPY,col="red",lty="dotted")
title(bquote("JPY_"~alpha*2))
#2nd regime
plot(rep(TAU[-9],2),c(heso_regime2_JPY[-9,][,3],heso_regime2_JPY[-9,][,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU[-9], rev(TAU[-9])), c(heso_regime2_JPY[-9,][,3], rev(heso_regime2_JPY[-9,][,4])), col = "LightSkyBlue")
points(TAU[-9], heso_regime2_JPY[-9,][,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU[-9], heso_regime2_JPY[-9,][,2], col = "blue")
abline(h=lr_JPY,col="red")
abline(h=lr_JPY-1.96*lr_set_JPY,col="red",lty="dotted")
abline(h=lr_JPY+1.96*lr_set_JPY,col="red",lty="dotted")
title(bquote("JPY_"~beta*2))

#CADUSD
#tau=0.1
tau=0.1
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.1<-out1$regime1
regime2_CAD_0.1<-out1$regime2
regime1_CAD_full_0.1<-out1$result1_full
regime2_CAD_full_0.1<-out1$result2_full
qhat_CAD_0.1<-out1$qhat
#tau=0.2
tau=0.2
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.2<-out1$regime1
regime2_CAD_0.2<-out1$regime2
regime1_CAD_full_0.2<-out1$result1_full
regime2_CAD_full_0.2<-out1$result2_full
qhat_CAD_0.2<-out1$qhat
#tau=0.3
tau=0.3
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.3<-out1$regime1
regime2_CAD_0.3<-out1$regime2
regime1_CAD_full_0.3<-out1$result1_full
regime2_CAD_full_0.3<-out1$result2_full
qhat_CAD_0.3<-out1$qhat
#tau0.4
tau=0.4
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.4<-out1$regime1
regime2_CAD_0.4<-out1$regime2
regime1_CAD_full_0.4<-out1$result1_full
regime2_CAD_full_0.4<-out1$result2_full
qhat_CAD_0.4<-out1$qhat
#tau=0.5
tau=0.5
out1 <- joint_thresh(y=y,x=x,q=q,1)
qhat_CAD_0.5<-out1$qhat
regime1_CAD_0.5<-out1$regime1
regime2_CAD_0.5<-out1$regime2
regime1_CAD_full_0.5<-out1$result1_full
regime2_CAD_full_0.5<-out1$result2_full
#tau=0.6
tau=0.6
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.6<-out1$regime1
regime2_CAD_0.6<-out1$regime2
regime1_CAD_full_0.6<-out1$result1_full
regime2_CAD_full_0.6<-out1$result2_full
qhat_CAD_0.6<-out1$qhat
#tau=0.7
tau=0.7
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.7<-out1$regime1
regime2_CAD_0.7<-out1$regime2
regime1_CAD_full_0.7<-out1$result1_full
regime2_CAD_full_0.7<-out1$result2_full
qhat_CAD_0.7<-out1$qhat
#tau=0.8
tau=0.8
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.8<-out1$regime1
regime2_CAD_0.8<-out1$regime2
regime1_CAD_full_0.8<-out1$result1_full
regime2_CAD_full_0.8<-out1$result2_full
qhat_CAD_0.8<-out1$qhat
#tau=0.9
tau=0.9
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CAD_0.9<-out1$regime1
regime2_CAD_0.9<-out1$regime2
regime1_CAD_full_0.9<-out1$result1_full
regime2_CAD_full_0.9<-out1$result2_full
qhat_CAD_0.9<-out1$qhat
#Estimation results
regime1_CAD_sta<-rbind(t(regime1_CAD_full_0.1),t(regime1_CAD_full_0.2),t(regime1_CAD_full_0.3),t(regime1_CAD_full_0.4),
                       t(regime1_CAD_full_0.5),t(regime1_CAD_full_0.6),t(regime1_CAD_full_0.7),t(regime1_CAD_full_0.8),
                       t(regime1_CAD_full_0.9))
regime2_CAD_sta<-rbind(t(regime2_CAD_full_0.1),t(regime2_CAD_full_0.2),t(regime2_CAD_full_0.3),t(regime2_CAD_full_0.4),
                       t(regime2_CAD_full_0.5),t(regime2_CAD_full_0.6),t(regime2_CAD_full_0.7),t(regime2_CAD_full_0.8),
                       t(regime2_CAD_full_0.9))
CAD_sta<-cbind(regime1_CAD_sta,regime2_CAD_sta)
#Graph of results
TAU<-seq(0.1,0.9,0.1)
qhat_CAD<-c(qhat_CAD_0.1,qhat_CAD_0.2,qhat_CAD_0.3,qhat_CAD_0.4,qhat_CAD_0.5,qhat_CAD_0.6,qhat_CAD_0.7,qhat_CAD_0.8,qhat_CAD_0.9)
bangthamchieu_CAD<-cbind(tau=TAU,qhat=qhat_CAD)
heso_regime1_CAD<-rbind(regime1_CAD_0.1[2,],regime1_CAD_0.2[2,],regime1_CAD_0.3[2,],regime1_CAD_0.4[2,],
                        regime1_CAD_0.5[2,],regime1_CAD_0.6[2,],regime1_CAD_0.7[2,],regime1_CAD_0.8[2,],regime1_CAD_0.9[2,])
heso_regime1_CAD<-cbind(TAU,heso_regime1_CAD)
heso_regime2_CAD<-rbind(regime2_CAD_0.1[2,],regime2_CAD_0.2[2,],regime2_CAD_0.3[2,],regime2_CAD_0.4[2,],
                        regime2_CAD_0.5[2,],regime2_CAD_0.6[2,],regime2_CAD_0.7[2,],regime2_CAD_0.8[2,],regime2_CAD_0.9[2,])
heso_regime2_CAD<-cbind(TAU,heso_regime2_CAD)
#1st regime
plot(rep(TAU,2),c(heso_regime1_CAD[,3],heso_regime1_CAD[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime1_CAD[,3], rev(heso_regime1_CAD[,4])), col = "LightSkyBlue")
points(TAU, heso_regime1_CAD[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime1_CAD[,2], col = "blue")
lr_CAD<-summary(lm(y~x))$coefficients[2,1]
lr_set_CAD<-summary(lm(y~x))$coefficients[2,2]
abline(h=lr_CAD,col="red")
abline(h=lr_CAD-1.96*lr_set_CAD,col="red",lty="dotted")
abline(h=lr_CAD+1.96*lr_set_CAD,col="red",lty="dotted")
title("CAD_beta_01")
#2nd regime
plot(rep(TAU,2),c(heso_regime2_CAD[,3],heso_regime2_CAD[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime2_CAD[,3], rev(heso_regime2_CAD[,4])), col = "LightSkyBlue")
points(TAU, heso_regime2_CAD[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime2_CAD[,2], col = "blue")
abline(h=lr_CAD,col="red")
abline(h=lr_CAD-1.96*lr_set_CAD,col="red",lty="dotted")
abline(h=lr_CAD+1.96*lr_set_CAD,col="red",lty="dotted")
title("CAD_beta_02")

#AUDUSD
#tau=0.1
tau=0.1
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.1<-out1$regime1
regime2_AUD_0.1<-out1$regime2
regime1_AUD_full_0.1<-out1$result1_full
regime2_AUD_full_0.1<-out1$result2_full
qhat_AUD_0.1<-out1$qhat
#tau=0.2
tau=0.2
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.2<-out1$regime1
regime2_AUD_0.2<-out1$regime2
regime1_AUD_full_0.2<-out1$result1_full
regime2_AUD_full_0.2<-out1$result2_full
qhat_AUD_0.2<-out1$qhat
#tau=0.3
tau=0.3
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.3<-out1$regime1
regime2_AUD_0.3<-out1$regime2
regime1_AUD_full_0.3<-out1$result1_full
regime2_AUD_full_0.3<-out1$result2_full
qhat_AUD_0.3<-out1$qhat
#tau0.4
tau=0.4
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.4<-out1$regime1
regime2_AUD_0.4<-out1$regime2
regime1_AUD_full_0.4<-out1$result1_full
regime2_AUD_full_0.4<-out1$result2_full
qhat_AUD_0.4<-out1$qhat
#tau=0.5
tau=0.5
out1 <- joint_thresh(y=y,x=x,q=q,1)
qhat_AUD_0.5<-out1$qhat
regime1_AUD_0.5<-out1$regime1
regime2_AUD_0.5<-out1$regime2
regime1_AUD_full_0.5<-out1$result1_full
regime2_AUD_full_0.5<-out1$result2_full
#tau=0.6
tau=0.6
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.6<-out1$regime1
regime2_AUD_0.6<-out1$regime2
regime1_AUD_full_0.6<-out1$result1_full
regime2_AUD_full_0.6<-out1$result2_full
qhat_AUD_0.6<-out1$qhat
#tau=0.7
tau=0.7
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.7<-out1$regime1
regime2_AUD_0.7<-out1$regime2
regime1_AUD_full_0.7<-out1$result1_full
regime2_AUD_full_0.7<-out1$result2_full
qhat_AUD_0.7<-out1$qhat
#tau=0.8
tau=0.8
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.8<-out1$regime1
regime2_AUD_0.8<-out1$regime2
regime1_AUD_full_0.8<-out1$result1_full
regime2_AUD_full_0.8<-out1$result2_full
qhat_AUD_0.8<-out1$qhat
#tau=0.9
tau=0.9
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_AUD_0.9<-out1$regime1
regime2_AUD_0.9<-out1$regime2
regime1_AUD_full_0.9<-out1$result1_full
regime2_AUD_full_0.9<-out1$result2_full
qhat_AUD_0.9<-out1$qhat
#Estimation results
regime1_AUD_sta<-rbind(t(regime1_AUD_full_0.1),t(regime1_AUD_full_0.2),t(regime1_AUD_full_0.3),t(regime1_AUD_full_0.4),
                       t(regime1_AUD_full_0.5),t(regime1_AUD_full_0.6),t(regime1_AUD_full_0.7),t(regime1_AUD_full_0.8),
                       t(regime1_AUD_full_0.9))
regime2_AUD_sta<-rbind(t(regime2_AUD_full_0.1),t(regime2_AUD_full_0.2),t(regime2_AUD_full_0.3),t(regime2_AUD_full_0.4),
                       t(regime2_AUD_full_0.5),t(regime2_AUD_full_0.6),t(regime2_AUD_full_0.7),t(regime2_AUD_full_0.8),
                       t(regime2_AUD_full_0.9))
AUD_sta<-cbind(regime1_AUD_sta,regime2_AUD_sta)
#Graph of results
TAU<-seq(0.1,0.9,0.1)
qhat_AUD<-c(qhat_AUD_0.1,qhat_AUD_0.2,qhat_AUD_0.3,qhat_AUD_0.4,qhat_AUD_0.5,qhat_AUD_0.6,qhat_AUD_0.7,qhat_AUD_0.8,qhat_AUD_0.9)
bangthamchieu_AUD<-cbind(tau=TAU,qhat=qhat_AUD)
heso_regime1_AUD<-rbind(regime1_AUD_0.1[2,],regime1_AUD_0.2[2,],regime1_AUD_0.3[2,],regime1_AUD_0.4[2,],
                        regime1_AUD_0.5[2,],regime1_AUD_0.6[2,],regime1_AUD_0.7[2,],regime1_AUD_0.8[2,],regime1_AUD_0.9[2,])
heso_regime1_AUD<-cbind(TAU,heso_regime1_AUD)
heso_regime2_AUD<-rbind(regime2_AUD_0.1[2,],regime2_AUD_0.2[2,],regime2_AUD_0.3[2,],regime2_AUD_0.4[2,],
                        regime2_AUD_0.5[2,],regime2_AUD_0.6[2,],regime2_AUD_0.7[2,],regime2_AUD_0.8[2,],regime2_AUD_0.9[2,])
heso_regime2_AUD<-cbind(TAU,heso_regime2_AUD)
#1st regime
plot(rep(TAU,2),c(heso_regime1_AUD[,3],heso_regime1_AUD[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n",ylim=c(-5,max(heso_regime1_AUD[,4])))
polygon(c(TAU, rev(TAU)), c(heso_regime1_AUD[,3], rev(heso_regime1_AUD[,4])), col = "LightSkyBlue")
points(TAU, heso_regime1_AUD[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime1_AUD[,2], col = "blue")
lr_AUD<-summary(lm(y~x))$coefficients[2,1]
lr_set_AUD<-summary(lm(y~x))$coefficients[2,2]
abline(h=lr_AUD,col="red")
abline(h=lr_AUD-1.96*lr_set_AUD,col="red",lty="dotted")
abline(h=lr_AUD+1.96*lr_set_AUD,col="red",lty="dotted")
title(bquote("AUD_"~alpha*2))
#2nd regime
plot(rep(TAU,2),c(heso_regime2_AUD[,3],heso_regime2_AUD[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n",ylim=c(-10,lr_AUD+1.96*lr_set_AUD))
polygon(c(TAU, rev(TAU)), c(heso_regime2_AUD[,3], rev(heso_regime2_AUD[,4])), col = "LightSkyBlue")
points(TAU, heso_regime2_AUD[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime2_AUD[,2], col = "blue")
abline(h=lr_AUD,col="red")
abline(h=lr_AUD-1.96*lr_set_AUD,col="red",lty="dotted")
abline(h=lr_AUD+1.96*lr_set_AUD,col="red",lty="dotted")
title(bquote("AUD_"~beta*2))

#EURUSD
#tau=0.1
tau=0.1
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.1<-out1$regime1
regime2_EUR_0.1<-out1$regime2
regime1_EUR_full_0.1<-out1$result1_full
regime2_EUR_full_0.1<-out1$result2_full
qhat_EUR_0.1<-out1$qhat
#tau=0.2
tau=0.2
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.2<-out1$regime1
regime2_EUR_0.2<-out1$regime2
regime1_EUR_full_0.2<-out1$result1_full
regime2_EUR_full_0.2<-out1$result2_full
qhat_EUR_0.2<-out1$qhat
#tau=0.3
tau=0.3
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.3<-out1$regime1
regime2_EUR_0.3<-out1$regime2
regime1_EUR_full_0.3<-out1$result1_full
regime2_EUR_full_0.3<-out1$result2_full
qhat_EUR_0.3<-out1$qhat
#tau0.4
tau=0.4
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.4<-out1$regime1
regime2_EUR_0.4<-out1$regime2
regime1_EUR_full_0.4<-out1$result1_full
regime2_EUR_full_0.4<-out1$result2_full
qhat_EUR_0.4<-out1$qhat
#tau=0.5
tau=0.5
out1 <- joint_thresh(y=y,x=x,q=q,1)
qhat_EUR_0.5<-out1$qhat
regime1_EUR_0.5<-out1$regime1
regime2_EUR_0.5<-out1$regime2
regime1_EUR_full_0.5<-out1$result1_full
regime2_EUR_full_0.5<-out1$result2_full
#tau=0.6
tau=0.6
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.6<-out1$regime1
regime2_EUR_0.6<-out1$regime2
regime1_EUR_full_0.6<-out1$result1_full
regime2_EUR_full_0.6<-out1$result2_full
qhat_EUR_0.6<-out1$qhat
#tau=0.7
tau=0.7
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.7<-out1$regime1
regime2_EUR_0.7<-out1$regime2
regime1_EUR_full_0.7<-out1$result1_full
regime2_EUR_full_0.7<-out1$result2_full
qhat_EUR_0.7<-out1$qhat
#tau=0.8
tau=0.8
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.8<-out1$regime1
regime2_EUR_0.8<-out1$regime2
regime1_EUR_full_0.8<-out1$result1_full
regime2_EUR_full_0.8<-out1$result2_full
qhat_EUR_0.8<-out1$qhat
#tau=0.9
tau=0.9
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_EUR_0.9<-out1$regime1
regime2_EUR_0.9<-out1$regime2
regime1_EUR_full_0.9<-out1$result1_full
regime2_EUR_full_0.9<-out1$result2_full
qhat_EUR_0.9<-out1$qhat
#Estimation results
regime1_EUR_sta<-rbind(t(regime1_EUR_full_0.1),t(regime1_EUR_full_0.2),t(regime1_EUR_full_0.3),t(regime1_EUR_full_0.4),
                       t(regime1_EUR_full_0.5),t(regime1_EUR_full_0.6),t(regime1_EUR_full_0.7),t(regime1_EUR_full_0.8),
                       t(regime1_EUR_full_0.9))
regime2_EUR_sta<-rbind(t(regime2_EUR_full_0.1),t(regime2_EUR_full_0.2),t(regime2_EUR_full_0.3),t(regime2_EUR_full_0.4),
                       t(regime2_EUR_full_0.5),t(regime2_EUR_full_0.6),t(regime2_EUR_full_0.7),t(regime2_EUR_full_0.8),
                       t(regime2_EUR_full_0.9))
EUR_sta<-cbind(regime1_EUR_sta,regime2_EUR_sta)
#Graph of results
TAU<-seq(0.1,0.9,0.1)
qhat_EUR<-c(qhat_EUR_0.1,qhat_EUR_0.2,qhat_EUR_0.3,qhat_EUR_0.4,qhat_EUR_0.5,qhat_EUR_0.6,qhat_EUR_0.7,qhat_EUR_0.8,qhat_EUR_0.9)
bangthamchieu_EUR<-cbind(tau=TAU,qhat=qhat_EUR)
heso_regime1_EUR<-rbind(regime1_EUR_0.1[2,],regime1_EUR_0.2[2,],regime1_EUR_0.3[2,],regime1_EUR_0.4[2,],
                        regime1_EUR_0.5[2,],regime1_EUR_0.6[2,],regime1_EUR_0.7[2,],regime1_EUR_0.8[2,],regime1_EUR_0.9[2,])
heso_regime1_EUR<-cbind(TAU,heso_regime1_EUR)
heso_regime2_EUR<-rbind(regime2_EUR_0.1[2,],regime2_EUR_0.2[2,],regime2_EUR_0.3[2,],regime2_EUR_0.4[2,],
                        regime2_EUR_0.5[2,],regime2_EUR_0.6[2,],regime2_EUR_0.7[2,],regime2_EUR_0.8[2,],regime2_EUR_0.9[2,])
heso_regime2_EUR<-cbind(TAU,heso_regime2_EUR)
#1st regime
plot(rep(TAU,2),c(heso_regime1_EUR[,3],heso_regime1_EUR[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime1_EUR[,3], rev(heso_regime1_EUR[,4])), col = "LightSkyBlue")
points(TAU, heso_regime1_EUR[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime1_EUR[,2], col = "blue")
lr_EUR<-summary(lm(y~x))$coefficients[2,1]
lr_set_EUR<-summary(lm(y~x))$coefficients[2,2]
abline(h=lr_EUR,col="red")
abline(h=lr_EUR-1.96*lr_set_EUR,col="red",lty="dotted")
abline(h=lr_EUR+1.96*lr_set_EUR,col="red",lty="dotted")
title("EUR_beta_01")
#2nd regime
plot(rep(TAU,2),c(heso_regime2_EUR[,3],heso_regime2_EUR[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime2_EUR[,3], rev(heso_regime2_EUR[,4])), col = "LightSkyBlue")
points(TAU, heso_regime2_EUR[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime2_EUR[,2], col = "blue")
abline(h=lr_EUR,col="red")
abline(h=lr_EUR-1.96*lr_set_EUR,col="red",lty="dotted")
abline(h=lr_EUR+1.96*lr_set_EUR,col="red",lty="dotted")
title("EUR_beta_02")

#CHFUSD
#tau=0.1
tau=0.1
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.1<-out1$regime1
regime2_CHF_0.1<-out1$regime2
regime1_CHF_full_0.1<-out1$result1_full
regime2_CHF_full_0.1<-out1$result2_full
qhat_CHF_0.1<-out1$qhat
#tau=0.2
tau=0.2
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.2<-out1$regime1
regime2_CHF_0.2<-out1$regime2
regime1_CHF_full_0.2<-out1$result1_full
regime2_CHF_full_0.2<-out1$result2_full
qhat_CHF_0.2<-out1$qhat
#tau=0.3
tau=0.3
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.3<-out1$regime1
regime2_CHF_0.3<-out1$regime2
regime1_CHF_full_0.3<-out1$result1_full
regime2_CHF_full_0.3<-out1$result2_full
qhat_CHF_0.3<-out1$qhat
#tau0.4
tau=0.4
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.4<-out1$regime1
regime2_CHF_0.4<-out1$regime2
regime1_CHF_full_0.4<-out1$result1_full
regime2_CHF_full_0.4<-out1$result2_full
qhat_CHF_0.4<-out1$qhat
#tau=0.5
tau=0.5
out1 <- joint_thresh(y=y,x=x,q=q,1)
qhat_CHF_0.5<-out1$qhat
regime1_CHF_0.5<-out1$regime1
regime2_CHF_0.5<-out1$regime2
regime1_CHF_full_0.5<-out1$result1_full
regime2_CHF_full_0.5<-out1$result2_full
#tau=0.6
tau=0.6
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.6<-out1$regime1
regime2_CHF_0.6<-out1$regime2
regime1_CHF_full_0.6<-out1$result1_full
regime2_CHF_full_0.6<-out1$result2_full
qhat_CHF_0.6<-out1$qhat
#tau=0.7
tau=0.7
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.7<-out1$regime1
regime2_CHF_0.7<-out1$regime2
regime1_CHF_full_0.7<-out1$result1_full
regime2_CHF_full_0.7<-out1$result2_full
qhat_CHF_0.7<-out1$qhat
#tau=0.8
tau=0.8
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.8<-out1$regime1
regime2_CHF_0.8<-out1$regime2
regime1_CHF_full_0.8<-out1$result1_full
regime2_CHF_full_0.8<-out1$result2_full
qhat_CHF_0.8<-out1$qhat
#tau=0.9
tau=0.9
out1 <- joint_thresh(y=y,x=x,q=q,1)
regime1_CHF_0.9<-out1$regime1
regime2_CHF_0.9<-out1$regime2
regime1_CHF_full_0.9<-out1$result1_full
regime2_CHF_full_0.9<-out1$result2_full
qhat_CHF_0.9<-out1$qhat
#Estimation results
regime1_CHF_sta<-rbind(t(regime1_CHF_full_0.1),t(regime1_CHF_full_0.2),t(regime1_CHF_full_0.3),t(regime1_CHF_full_0.4),
                       t(regime1_CHF_full_0.5),t(regime1_CHF_full_0.6),t(regime1_CHF_full_0.7),t(regime1_CHF_full_0.8),
                       t(regime1_CHF_full_0.9))
regime2_CHF_sta<-rbind(t(regime2_CHF_full_0.1),t(regime2_CHF_full_0.2),t(regime2_CHF_full_0.3),t(regime2_CHF_full_0.4),
                       t(regime2_CHF_full_0.5),t(regime2_CHF_full_0.6),t(regime2_CHF_full_0.7),t(regime2_CHF_full_0.8),
                       t(regime2_CHF_full_0.9))
CHF_sta<-cbind(regime1_CHF_sta,regime2_CHF_sta)
#Graph of results
TAU<-seq(0.1,0.9,0.1)
qhat_CHF<-c(qhat_CHF_0.1,qhat_CHF_0.2,qhat_CHF_0.3,qhat_CHF_0.4,qhat_CHF_0.5,qhat_CHF_0.6,qhat_CHF_0.7,qhat_CHF_0.8,qhat_CHF_0.9)
bangthamchieu_CHF<-cbind(tau=TAU,qhat=qhat_CHF)
heso_regime1_CHF<-rbind(regime1_CHF_0.1[2,],regime1_CHF_0.2[2,],regime1_CHF_0.3[2,],regime1_CHF_0.4[2,],
                        regime1_CHF_0.5[2,],regime1_CHF_0.6[2,],regime1_CHF_0.7[2,],regime1_CHF_0.8[2,],regime1_CHF_0.9[2,])
heso_regime1_CHF<-cbind(TAU,heso_regime1_CHF)
heso_regime2_CHF<-rbind(regime2_CHF_0.1[2,],regime2_CHF_0.2[2,],regime2_CHF_0.3[2,],regime2_CHF_0.4[2,],
                        regime2_CHF_0.5[2,],regime2_CHF_0.6[2,],regime2_CHF_0.7[2,],regime2_CHF_0.8[2,],regime2_CHF_0.9[2,])
heso_regime2_CHF<-cbind(TAU,heso_regime2_CHF)
#1st regime
plot(rep(TAU,2),c(heso_regime1_CHF[,3],heso_regime1_CHF[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(heso_regime1_CHF[,3], rev(heso_regime1_CHF[,4])), col = "LightSkyBlue")
points(TAU, heso_regime1_CHF[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime1_CHF[,2], col = "blue")
lr_CHF<-summary(lm(y~x))$coefficients[2,1]
lr_set_CHF<-summary(lm(y~x))$coefficients[2,2]
abline(h=lr_CHF,col="red")
abline(h=lr_CHF-1.96*lr_set_CHF,col="red",lty="dotted")
abline(h=lr_CHF+1.96*lr_set_CHF,col="red",lty="dotted")
title(bquote("CHF_"~alpha*2))
#2nd regime
plot(rep(TAU,2),c(heso_regime2_CHF[,3],heso_regime2_CHF[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n",ylim = c(lr_CHF-1.96*lr_set_CHF,max(heso_regime2_CHF[,4])))
polygon(c(TAU, rev(TAU)), c(heso_regime2_CHF[,3], rev(heso_regime2_CHF[,4])), col = "LightSkyBlue")
points(TAU, heso_regime2_CHF[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, heso_regime2_CHF[,2], col = "blue")
abline(h=lr_CHF,col="red")
abline(h=lr_CHF-1.96*lr_set_CHF,col="red",lty="dotted")
abline(h=lr_CHF+1.96*lr_set_CHF,col="red",lty="dotted")
title(bquote("CHF_"~beta*2))

#6. Quantile regression results--------------------------------------------------------
coef_0.1<-summary(rq(y~x,tau=0.1))$coefficients
coef_0.2<-summary(rq(y~x,tau=0.2))$coefficients
coef_0.3<-summary(rq(y~x,tau=0.3))$coefficients
coef_0.4<-summary(rq(y~x,tau=0.4))$coefficients
coef_0.5<-summary(rq(y~x,tau=0.5))$coefficients
coef_0.6<-summary(rq(y~x,tau=0.6))$coefficients
coef_0.7<-summary(rq(y~x,tau=0.7))$coefficients
coef_0.8<-summary(rq(y~x,tau=0.8))$coefficients
coef_0.9<-summary(rq(y~x,tau=0.9))$coefficients
TAU<-seq(0.1,0.9,0.1)
coef_all<-rbind(coef_0.1[3,],coef_0.2[3,],coef_0.3[3,],coef_0.4[3,],coef_0.5[3,],
                coef_0.6[3,],coef_0.7[3,],coef_0.8[3,],coef_0.9[3,])
coef_all<-cbind(TAU,coef_all)

plot(rep(TAU,2),c(coef_all[,3],coef_all[,4]),xlab = "quantiles", ylab = "Coefficient", type = "n")
polygon(c(TAU, rev(TAU)), c(coef_all[,3], rev(coef_all[,4])), col = "LightSkyBlue")
points(TAU, coef_all[,2], cex = 0.5, pch = "o", col = "blue")
lines(TAU, coef_all[,2], col = "blue")
lr<-summary(lm(y~x))$coefficients[3,1]
lr_set<-summary(lm(y~x))$coefficients[3,2]
abline(h=lr,col="red")
abline(h=lr-1.96*lr_set,col="red",lty="dotted")
abline(h=lr+1.96*lr_set,col="red",lty="dotted")

#7. Conditional GDP-Skewed T Student: Distribution Fitting and Risk Measure of Exchange rate movements-----------------------------------
library(quantreg)
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(pracma)
library(readr)
library(sn)
#data input (run a group of codes after another to get the results for each pair of currency)
#CADUSD
author_1_q5 <- rq(CADUSD_return~fw_CADUSD3m+EPU, tau = .05, data = df1_CAD)
author_1_q25 <- rq(CADUSD_return~fw_CADUSD3m+EPU, tau = .25, data = df1_CAD)
author_1_q75 <- rq(CADUSD_return~fw_CADUSD3m+EPU, tau = .75, data = df1_CAD)
author_1_q95 <- rq(CADUSD_return~fw_CADUSD3m+EPU, tau = .95, data = df1_CAD)
# #CHFUSD
# author_1_q5 <- rq(CHFUSD_return~fw_CHFUSD3m+EPU, tau = .05, data = df1_CHF)
# author_1_q25 <- rq(CHFUSD_return~fw_CHFUSD3m+EPU, tau = .25, data = df1_CHF)
# author_1_q75 <- rq(CHFUSD_return~fw_CHFUSD3m+EPU, tau = .75, data = df1_CHF)
# author_1_q95 <- rq(CHFUSD_return~fw_CHFUSD3m+EPU, tau = .95, data = df1_CHF)
# #GBPUSD
# author_1_q5 <- rq(GBPUSD_return~fw_GBPUSD3m+EPU, tau = .05, data = df1_GBP)
# author_1_q25 <- rq(GBPUSD_return~fw_GBPUSD3m+EPU, tau = .25, data = df1_GBP)
# author_1_q75 <- rq(GBPUSD_return~fw_GBPUSD3m+EPU, tau = .75, data = df1_GBP)
# author_1_q95 <- rq(GBPUSD_return~fw_GBPUSD3m+EPU, tau = .95, data = df1_GBP)
# #AUDUSD
# author_1_q5 <- rq(AUDUSD_return~fw_AUDUSD3m+EPU, tau = .05, data = df1_AUD)
# author_1_q25 <- rq(AUDUSD_return~fw_AUDUSD3m+EPU, tau = .25, data = df1_AUD)
# author_1_q75 <- rq(AUDUSD_return~fw_AUDUSD3m+EPU, tau = .75, data = df1_AUD)
# author_1_q95 <- rq(AUDUSD_return~fw_AUDUSD3m+EPU, tau = .95, data = df1_AUD)
# #JPYUSD
# author_1_q5 <- rq(JPYUSD_return~fw_JPYUSD3m+EPU, tau = .05, data = df1_JPY)
# author_1_q25 <- rq(JPYUSD_return~fw_JPYUSD3m+EPU, tau = .25, data = df1_JPY)
# author_1_q75 <- rq(JPYUSD_return~fw_JPYUSD3m+EPU, tau = .75, data = df1_JPY)
# author_1_q95 <- rq(JPYUSD_return~fw_JPYUSD3m+EPU, tau = .95, data = df1_JPY)

#Fitted values
fitted_valuesP <- data.frame("Q5" = author_1_q5$fitted.values,
                             "Q25" = author_1_q25$fitted.values,
                             "Q75" = author_1_q75$fitted.values,
                             "Q95" = author_1_q95$fitted.values)
#parameter storage
parameters <- data.frame("location" = rep(NA, nrow(fitted_valuesP)),
                         "scale" = rep(NA, nrow(fitted_valuesP)),
                         "shape" = rep(NA, nrow(fitted_valuesP)),
                         "freedom" = rep(NA, nrow(fitted_valuesP)),
                         "DATE" = df2_CAD$index)
#Loss function 5, 25, 75 and 95 percent quantiles
Loss <- function(x, q, freedom){
  
  distance <- q-sn::qst(c(0.05,0.25,0.75,0.95),
                        xi = x[1], 
                        omega = x[2],
                        alpha = x[3],
                        nu = freedom)
  
  return(sum(distance^2))
}

# One Quarter ahead Parameter Estimation
print("Estimating One Quarter Ahead Distribution")
for(i in 1:nrow(fitted_valuesP)){
  
  print(paste("Estimating case ", i, " out of ", nrow(fitted_valuesP)))
  
  e <- c(fitted_valuesP$Q5[i], fitted_valuesP$Q25[i], 
         fitted_valuesP$Q75[i], fitted_valuesP$Q95[i])
  
  sol <- 1e10
  
  for (f in 1:30){
    
    opt <- optim(c(0,1,0),
                 Loss, 
                 method = "L-BFGS-B",
                 lower = c(-20, 0, -30),
                 upper = c(20, 50, 30),
                 q=e,
                 freedom = f)
    
    if (as.numeric(opt$value) < sol){
      
      sol <- opt$value
      
      parameters$location[i] <- opt$par[1]
      parameters$scale[i] <- opt$par[2]
      parameters$shape[i] <- opt$par[3]
      parameters$ freedom[i] <- f
      
    }
  }
}

#storage the results of time-varying skew-t distributions' parameters for each pair of currency
parameters_CAD3month <- parameters
# parameters_UK3month <- parameters
# parameters_CHF3month <- parameters
# parameters_AUD3month <- parameters
# parameters_EUR3month <- parameters
# parameters_JPY3month <- parameters

#Calculate Expected Shortfall for GBPUSD
Esf <- c()
Elr <- c()
alpha <- 0.05
delta1 <- 0.01

for (i in 1:nrow(parameters_UK3month)){
  
  Esf[i] <- 1/alpha * sum( qst(c(0.01, 0.02, 0.03, 0.04, 0.05), 
                               xi = parameters_UK3month$location[i], 
                               omega = parameters_UK3month$scale[i], 
                               alpha = parameters_UK3month$shape[i], 
                               nu = parameters_UK3month$freedom[i]) 
                           * delta1)
  Elr[i] <- 1/alpha * sum( qst(c(0.95, 0.96, 0.97, 0.98, 0.99), 
                               xi = parameters_UK3month$location[i], 
                               omega = parameters_UK3month$scale[i], 
                               alpha = parameters_UK3month$shape[i], 
                               nu = parameters_UK3month$freedom[i]) 
                           * delta1)
}
risk_measures <- data.frame("Esf" = Esf, "Elr"= Elr, DATE = parameters_UK3month$DATE)

plot1 <- ggplot(risk_measures, aes(x = DATE)) +
  geom_line(aes(y = Esf, colour = "Expected shortfall"), linetype = 1, size = 0.2) +
  geom_line(aes(y = Elr, colour = "Expected longrise"), linetype = 1, size=0.2)+
  ylab("") + xlab("") + ggtitle("US dollar against UK pound sterling") +
  scale_colour_manual("", 
                      breaks = c("Expected shortfall","Expected longrise"),
                      values = c("blue","red"))+
  theme_bw()+ theme(legend.position="bottom",plot.title=element_text(size=10,hjust=0.5))

#Calculate Expected Shortfall CHF/USD
Esf <- c()
Elr <- c()
alpha <- 0.05
delta1 <- 0.01

for (i in 1:nrow(parameters_CHF3month)){
  
  Esf[i] <- 1/alpha * sum( qst(c(0.01, 0.02, 0.03, 0.04, 0.05), 
                               xi = parameters_CHF3month$location[i], 
                               omega = parameters_CHF3month$scale[i], 
                               alpha = parameters_CHF3month$shape[i], 
                               nu = parameters_CHF3month$freedom[i]) 
                           * delta1)
  Elr[i] <- 1/alpha * sum( qst(c(0.95, 0.96, 0.97, 0.98, 0.99), 
                               xi = parameters_CHF3month$location[i], 
                               omega = parameters_CHF3month$scale[i], 
                               alpha = parameters_CHF3month$shape[i], 
                               nu = parameters_CHF3month$freedom[i]) 
                           * delta1)
}
risk_measures <- data.frame("Esf" = Esf, "Elr"= Elr, DATE = parameters_CHF3month$DATE)

plot2 <- ggplot(risk_measures, aes(x = DATE)) +
  geom_line(aes(y = Esf, colour = "Expected shortfall"), linetype = 1, size = 0.2) +
  geom_line(aes(y = Elr, colour = "Expected longrise"), linetype = 1, size=0.2)+
  ylab("") + xlab("") + ggtitle("US dollar against Swiss franc") +
  scale_colour_manual("", 
                      breaks = c("Expected shortfall","Expected longrise"),
                      values = c("blue","red"))+
  theme_bw()+ theme(legend.position="bottom",plot.title=element_text(size=10,hjust=0.5))


#Calculate Expected Shortfall AUD/USD
Esf <- c()
Elr <- c()
alpha <- 0.05
delta1 <- 0.01

for (i in 1:nrow(parameters_AUD3month)){
  
  Esf[i] <- 1/alpha * sum( qst(c(0.01, 0.02, 0.03, 0.04, 0.05), 
                               xi = parameters_AUD3month$location[i], 
                               omega = parameters_AUD3month$scale[i], 
                               alpha = parameters_AUD3month$shape[i], 
                               nu = parameters_AUD3month$freedom[i]) 
                           * delta1)
  Elr[i] <- 1/alpha * sum( qst(c(0.95, 0.96, 0.97, 0.98, 0.99), 
                               xi = parameters_AUD3month$location[i], 
                               omega = parameters_AUD3month$scale[i], 
                               alpha = parameters_AUD3month$shape[i], 
                               nu = parameters_AUD3month$freedom[i]) 
                           * delta1)
}
risk_measures <- data.frame("Esf" = Esf, "Elr"= Elr, DATE = parameters_AUD3month$DATE)

plot4 <- ggplot(risk_measures, aes(x = DATE)) +
  geom_line(aes(y = Esf, colour = "Expected shortfall"), linetype = 1, size = 0.2) +
  geom_line(aes(y = Elr, colour = "Expected longrise"), linetype = 1, size=0.2)+
  ylab("") + xlab("") + ggtitle("US dollar against Australian dollar") +
  scale_colour_manual("", 
                      breaks = c("Expected shortfall","Expected longrise"),
                      values = c("blue","red"))+
  theme_bw()+ theme(legend.position="bottom",plot.title=element_text(size=10,hjust=0.5))

#Calculate Expected Shortfall JPY/USD
Esf <- c()
Elr <- c()
alpha <- 0.05
delta1 <- 0.01

for (i in 1:nrow(parameters_JPY3month)){
  
  Esf[i] <- 1/alpha * sum( qst(c(0.01, 0.02, 0.03, 0.04, 0.05), 
                               xi = parameters_JPY3month$location[i], 
                               omega = parameters_JPY3month$scale[i], 
                               alpha = parameters_JPY3month$shape[i], 
                               nu = parameters_JPY3month$freedom[i]) 
                           * delta1)
  Elr[i] <- 1/alpha * sum( qst(c(0.95, 0.96, 0.97, 0.98, 0.99), 
                               xi = parameters_JPY3month$location[i], 
                               omega = parameters_JPY3month$scale[i], 
                               alpha = parameters_JPY3month$shape[i], 
                               nu = parameters_JPY3month$freedom[i]) 
                           * delta1)
}
risk_measures <- data.frame("Esf" = Esf, "Elr"= Elr, DATE = parameters_JPY3month$DATE)

plot4 <- ggplot(risk_measures, aes(x = DATE)) +
  geom_line(aes(y = Esf, colour = "Expected shortfall"), linetype = 1, size = 0.2) +
  geom_line(aes(y = Elr, colour = "Expected longrise"), linetype = 1, size=0.2)+
  ylab("") + xlab("") + ggtitle("US dollar against Japanese yen") +
  scale_colour_manual("", 
                      breaks = c("Expected shortfall","Expected longrise"),
                      values = c("blue","red"))+
  theme_bw()+ theme(legend.position="bottom",plot.title=element_text(size=10,hjust=0.5))

#Calculate Expected Shortfall CAD/USD
Esf <- c()
Elr <- c()
alpha <- 0.05
delta1 <- 0.01

for (i in 1:nrow(parameters_CAD3month)){
  
  Esf[i] <- 1/alpha * sum( qst(c(0.01, 0.02, 0.03, 0.04, 0.05), 
                               xi = parameters_CAD3month$location[i], 
                               omega = parameters_CAD3month$scale[i], 
                               alpha = parameters_CAD3month$shape[i], 
                               nu = parameters_CAD3month$freedom[i]) 
                           * delta1)
  Elr[i] <- 1/alpha * sum( qst(c(0.95, 0.96, 0.97, 0.98, 0.99), 
                               xi = parameters_CAD3month$location[i], 
                               omega = parameters_CAD3month$scale[i], 
                               alpha = parameters_CAD3month$shape[i], 
                               nu = parameters_CAD3month$freedom[i]) 
                           * delta1)
}
risk_measures <- data.frame("Esf" = Esf, "Elr"= Elr, DATE = parameters_CAD3month$DATE)

plot5 <- ggplot(risk_measures, aes(x = DATE)) +
  geom_line(aes(y = Esf, colour = "Expected shortfall"), linetype = 1, size = 0.2) +
  geom_line(aes(y = Elr, colour = "Expected longrise"), linetype = 1, size=0.2)+
  ylab("") + xlab("") + ggtitle("US dollar against Canadian dollar") +
  scale_colour_manual("", 
                      breaks = c("Expected shortfall","Expected longrise"),
                      values = c("blue","red"))+
  theme_bw()+ theme(legend.position="bottom",plot.title=element_text(size=10,hjust=0.5))