library(readxl)
library(vtable)
library(tidyr)
library(ggplot2)
library(neverhpfilter)
library(xts)
library(corx)
#The Impact of EPU on US bond markets------------------------------------------
#Import data-------------------------------------------
data <- read_excel("USA_finaldata.xlsx",sheet="Data")
data <- data[complete.cases(data),]

# data <- read_excel("UK_finaldata.xlsx",sheet="Data")
# data <- read_excel("Japan_finaldata.xlsx",sheet="Data")
#Graph of short and long term yields-------------------------------------------
yield1 <- data.frame(data$date,data$yield_1y,data$yield_2y,data$yield_3y,
                      data$yield_5y,data$yield_7y,data$yield_10y,
                     data$yield_20y,data$yield_30y)

visual1 <- gather(yield1,variable,yield,data.yield_1y:data.yield_30y)

ggplot(visual1,aes(x=data.date,y=yield,colour=variable,group=variable))+geom_line()

#Graph of short and long term bond volatilities--------------------------------
vol1 <- data.frame(data$date,data$vol_1y,data$vol_2y,data$vol_3y,data$vol_5y,
                     data$vol_7y,data$vol_10y,data$vol_20y,data$vol_30y,data$EPU)
visual2 <- gather(vol1,variable,yield,data.vol_1y:data.EPU)
ggplot(visual2,aes(x=data.date,y=yield,colour=variable,group=variable))+geom_line()
#model estimation
#1.Unit root test--------------------------------------------------------------
install.packages("remotes")
remotes::install_github("chaoguo14/crowbar")
library('remotes')
adf_test_summary <- function(df, p_select = "BIC") {
  #' Run ADF test for unit roots on all columns of a data frame.
  #'
  #' @param p_select Method for selecting number of lags to include in the test regression. Possible values: "AIC", "BIC"
  #'
  #' @return A data frame with 3 rows and ncol(df) columns.
  #'
  #' @examples
  #' # Example 1: We test it on iid white noises and a random walk.
  #' # As expected, white_noise's p-value is smaller than 0.01, so it is stationary. But random walk is not.
  #' set.seed(10)
  #' df <- data.frame(white_noise = rnorm(100), random_walk = cumsum(rnorm(100)))
  #' adf_test_summary(df)
  #'
  #' # Example 2: We test it on a "trend stationary" time series.
  #' set.seed(10)
  #' df <- data.frame(trend_s = 0.8 * seq(1:20) + rnorm(20))
  #' adf_test_summary(df)
  type_explanation <- c("no intercept, no time trend",
                        "with intercept, no time trend",
                        "with intercept, with linear time trend")
  names(type_explanation) <- c("none", "drift", "trend")
  
  urca_to_adftest <- c("nc", "c", "ct")
  names(urca_to_adftest) <- c("none","drift","trend")
  
  cat(paste("Augmented Unit-Root Test.\n"))
  # cat(paste("Type: ", type_explanation[type], "\n"))
  cat("  [1]. p-value < 0.05 implies no unit root exists (i.e. stationary).\n")
  cat("  [2]. Use domain knowledge to determine if including intercept/trend is appropriate.\n\n")
  
  data.frame(
    lapply(df, FUN = function(x) c("None" = fUnitRoots::adfTest(x,
                                                                lags = urca::ur.df(x, type = "none", selectlags = p_select)@lags,
                                                                type = urca_to_adftest["none"])@test$p.value,
                                   "Intercept" = fUnitRoots::adfTest(x,
                                                                     lags = urca::ur.df(x, type = "drift", selectlags = p_select)@lags,
                                                                     type = urca_to_adftest["drift"])@test$p.value,
                                   "Trend" = fUnitRoots::adfTest(x,
                                                                 lags = urca::ur.df(x, type = "trend", selectlags = p_select)@lags,
                                                                 type = urca_to_adftest["trend"])@test$p.value)
    )
  )
}
adf_test_summary(data, p_select = "BIC")
#2. Hamilton filter----------------------------------------------------------
data$lnEPU <- ln(data$EPU)
data$EPU <- xts(x=data$EPU, order.by=data$date)
EPU_cyc <- yth_filter(data$EPU,h=24,p=12,output="cycle",family = gaussian)

#convert xts to dataframe
xts2data.frame = function(x) {
  df = data.frame(x=index(x), coredata(x))
  colnames(df) = c("date", colnames(x))
  return(df)
}
EPU_cyc <- xts2data.frame(EPU_cyc)

#merge by date
data$my <-format(as.Date(data$date),"%Y-%m")
EPU_cyc$my <- format(as.Date(EPU_cyc$date),"%Y-%m")
data_regression <- merge(data,EPU_cyc, by ='my')
#convert into time series
data_regression <- ts(data_regression)
#3. Uncertainty and the yield curve----------------------------------------------
result <- rbind(data_regression %>% lm(data=.,yield_1y ~ .lnEPU) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_2y ~ .lnEPU ) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_3y ~ .lnEPU) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_5y ~ .lnEPU) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_7y ~ .lnEPU) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_10y ~ .lnEPU) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_20y ~ .lnEPU) %>% coef() %>% .[2],
                data_regression %>% lm(data=.,yield_30y ~ .lnEPU) %>% coef() %>% .[2])
maturity <- c(1,2,3,5,7,10,20,30)
install.packages('ggplot2')
library(ggplot2)
ggplot(data.frame(maturity,result),aes(x=maturity,y=result))+geom_line()+geom_point()

#control for macro and financial variable
result2 <- rbind(data_regression %>% lm(data=.,yield_1y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_2y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_3y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_5y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_7y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_10y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_20y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,yield_30y ~ .lnEPU + Log_DY+ Log_IP + Log_CPI ) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result2),aes(x=maturity,y=result2))+geom_line()+geom_point()
#data %>% lm(data=.,yield_2y ~ EPU + Log_DY+ Log_IP + Log_CPI ) %>% summary()
#data %>% lm(data=.,yield_30y ~ EPU  ) %>% summary()
#NeweyWest
fit <- data %>% lm(data=.,yield_1y ~ EPU  )
NeweyWest(fit, lag = 5, order.by = NULL, prewhite = TRUE, adjust = FALSE, 
          diagnostics = FALSE, sandwich = TRUE, ar.method = "ols", data = list(),
          verbose = FALSE)
plot(density(data$EPU))
install.packages('sandwich')
library(sandwich)

#4. Uncertaity and Bond volatility-----------------------------------------------
result3 <- rbind(data_regression %>% lm(data=.,vol_1y ~ .lnEPU) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_2y ~ .lnEPU ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_3y ~ .lnEPU) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_5y ~ .lnEPU) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_7y ~ .lnEPU) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_10y ~ .lnEPU) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_20y ~ .lnEPU) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_30y ~ .lnEPU) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result3),aes(x=maturity,y=result3))+geom_line()+geom_point()

#with control
result4 <- rbind(data_regression %>% lm(data=.,vol_1y ~ .lnEPU +Log_DY+Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_2y ~ .lnEPU +Log_DY+Log_IP + Log_CPI ) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_3y ~ .lnEPU +Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_5y ~ .lnEPU +Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_7y ~ .lnEPU +Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_10y ~ .lnEPU +Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_20y ~ .lnEPU +Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 data_regression %>% lm(data=.,vol_30y ~ .lnEPU +Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result4),aes(x=maturity,y=result4))+geom_line()+geom_point()


#5. Uncertainty and Bond Excess return-----------------------------------------------
#5.1. 1-month excess return-----------------------------------------------
#graph
data1m <- read_excel("USA_finaldata.xlsx",sheet="1-month-data")
basicStats(data1m[,c(2,3,4,5,6,7,8,9)])

excess1 <- data.frame(data1m$date,data1m$year_1,data1m$year_2,data1m$year_3,
                      data1m$year_5,data1m$year_7,data1m$year_10,data1m$year_20,
                      data1m$year_30)
visual1 <- gather(excess1,variable,yield,data1m.year_1:data1m.year_30)
ggplot(visual1,aes(x=data1m.date,y=yield,colour=variable,group=variable))+geom_line()+facet_wrap(.~variable)+theme_grey()
#Correlation matrix
corx(data=data1m[,c(-1)])

#Hamilton filter
data1m$lnEPU <- ln(data1$EPU)
data1m$EPU <- xts(x=data1m$EPU, order.by=data1m$date)
EPU_cyc <- yth_filter(data1m$EPU,h=24,p=12,output="cycle",family = gaussian)

#convert xts to dataframe
EPU_cyc <- xts2data.frame(EPU_cyc)
data1m$EPU <- xts2data.frame(data1m$EPU)
data1m$EPU <- data1m$EPU[,2]
#Merging by date
data1m$my <-format(as.Date(data1m$date),"%Y-%m")
EPU_cyc$my <- format(as.Date(EPU_cyc$date),"%Y-%m")
excess1m_regression <- merge(data1m,EPU_cyc, by ='my')
#Convert into time series
excess1m_regression<-ts(excess1m_regression)
#result
result1m <- rbind(excess1m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess1m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result1m),aes(x=maturity,y=result1m))+geom_line()+geom_point()
excess1m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature) %>% summary()
#Control for macro and financial variables
result1m <- rbind(excess1m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                 excess1m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result1m),aes(x=maturity,y=result1m))+geom_line()+geom_point()
excess1m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% summary()
#5.2. 3-month excess return--------------------------------------------------
#Graph
data3m <- read_excel("USA_finaldata.xlsx",sheet="3-month-data")
basicStats(data3m[,c(2,3,4,5,6,7,8,9)])
excess3 <- data.frame(data3m$date,data3m$year_1,data3m$year_2,data3m$year_3,
                      data3m$year_5,data3m$year_7,data3m$year_10,data3m$year_20,
                      data3m$year_30)
visual3 <- gather(excess3,variable,yield,data3m.year_1:data3m.year_30)
ggplot(visual3,aes(x=data3m.date,y=yield,colour=variable,group=variable))+geom_line()+facet_wrap(.~variable)
#Hamilton Filter
data3m$lnEPU <- ln(data3m$EPU)
data3m$EPU <- xts(x=data3m$EPU, order.by=data3m$date)
EPU_cyc <- yth_filter(data3m$EPU,h=24,p=12,output="cycle",family = gaussian)

#convert xts to dataframe
EPU_cyc <- xts2data.frame(EPU_cyc)

#Merging by date
data3m$my <-format(as.Date(data3m$date),"%Y-%m")
EPU_cyc$my <- format(as.Date(EPU_cyc$date),"%Y-%m")
excess3m_regression <- merge(data3m,EPU_cyc, by ='my')
#Convert into time series
excess3m_regression<-ts(excess3m_regression)
#result
result3m <- rbind(excess3m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result3m),aes(x=maturity,y=result3m))+geom_line()+geom_point()
#Control for macro and financial variables
result3m <- rbind(excess3m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess3m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2])
excess3m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% summary()
#5.3. 6-month excess return--------------------------------------------------------
#Graph
data6m <- read_excel("USA_finaldata.xlsx",sheet="6-month-data")
basicStats(data6m[,c(2,3,4,5,6,7,8,9)])
excess6 <- data.frame(data6m$date,data6m$year_1,data6m$year_2,data6m$year_3,
                      data6m$year_5,data6m$year_7,data6m$year_10,data6m$year_20,
                      data6m$year_30)
visual6 <- gather(excess6,variable,yield,data6m.year_1:data6m.year_30)
ggplot(visual6,aes(x=data6m.date,y=yield,colour=variable,group=variable))+geom_line()+facet_wrap(.~variable)
#Hamilton Filter
data6m$lnEPU <- ln(data6m$EPU)
data6m$EPU <- xts(x=data6m$EPU, order.by=data6m$date)
EPU_cyc <- yth_filter(data6m$EPU,h=24,p=12,output="cycle",family = gaussian)

#convert xts to dataframe
str(data6m)
EPU_cyc <- xts2data.frame(EPU_cyc)
#Merging by date
data6m$my <-format(as.Date(data6m$date),"%Y-%m")
EPU_cyc$my <- format(as.Date(EPU_cyc$date),"%Y-%m")
excess6m_regression <- merge(data6m,EPU_cyc, by ='my')
#Convert into time series
excess6m_regression<-ts(excess6m_regression)
#result
result6m <- rbind(excess6m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result6m),aes(x=maturity,y=result6m))+geom_line()+geom_point()
#Control for macro and financial variables
result6m <- rbind(excess6m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess6m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2])
excess6m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% summary()
#5.4. 9-month excess return---------------------------------------------------------
#Graph
data9m <- read_excel("USA_finaldata.xlsx",sheet="9-month-data")
basicStats(data9m[,c(2,3,4,5,6,7,8,9)])
excess9 <- data.frame(data9m$date,data9m$year_1,data9m$year_2,data9m$year_3,
                      data9m$year_5,data9m$year_7,data9m$year_10,data9m$year_20,
                      data9m$year_30)
visual9 <- gather(excess9,variable,yield,data9m.year_1:data9m.year_30)
ggplot(visual9,aes(x=data9m.date,y=yield,colour=variable,group=variable))+geom_line()+facet_wrap(.~variable)
#Hamilton Filter
data9m$lnEPU <- ln(data9m$EPU)
data9m$EPU <- xts(x=data9m$EPU, order.by=data9m$date)
EPU_cyc <- yth_filter(data9m$EPU,h=24,p=12,output="cycle",family = gaussian)
#convert xts to dataframe
EPU_cyc <- xts2data.frame(EPU_cyc)
#Merging by date
data9m$my <-format(as.Date(data9m$date),"%Y-%m")
EPU_cyc$my <- format(as.Date(EPU_cyc$date),"%Y-%m")
excess9m_regression <- merge(data9m,EPU_cyc, by ='my')
#Convert into time series
excess9m_regression<-ts(excess9m_regression)
#result
result9m <- rbind(excess9m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result9m),aes(x=maturity,y=result9m))+geom_line()+geom_point()
#Control for macro and financial variables
result9m <- rbind(excess9m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2])
excess9m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% summary()
#5.5. 12-month excess return----------------------------------------------------------
#Graph
data12m <- read_excel("USA_finaldata.xlsx",sheet="12-month-data")
basicStats(data12m[,c(2,3,4,5,6,7,8,9)])
excess12 <- data.frame(data12m$date,data12m$year_1,data12m$year_2,data12m$year_3,
                       data12m$year_5,data12m$year_7,data12m$year_10,data12m$year_20,
                       data12m$year_30)
visual12 <- gather(excess12,variable,yield,data12m.year_1:data12m.year_30)
ggplot(visual12,aes(x=data12m.date,y=yield,colour=variable,group=variable))+geom_line()+facet_wrap(.~variable)
#Hamilton Filter
data12m$lnEPU <- ln(data12m$EPU)
data12m$EPU <- xts(x=data12m$EPU, order.by=data12m$date)
EPU_cyc <- yth_filter(data12m$EPU,h=24,p=12,output="cycle",family = gaussian)
#convert xts to dataframe
EPU_cyc <- xts2data.frame(EPU_cyc)

#Merging by date
data12m$my <-format(as.Date(data12m$date),"%Y-%m")
EPU_cyc$my <- format(as.Date(EPU_cyc$date),"%Y-%m")
excess12m_regression <- merge(data12m,EPU_cyc, by ='my')
#Convert into time series
excess12m_regression<-ts(excess12m_regression)
#result
result12m <- rbind(excess12m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2],
                  excess12m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature) %>% coef() %>% .[2])
ggplot(data.frame(maturity,result12m),aes(x=maturity,y=result12m))+geom_line()+geom_point()
#Control for macro and financial variables
result12m <- rbind(excess9m_regression %>% lm(data=.,year_1 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_2 ~ .lnEPU + Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_3 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_5 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_7 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_10 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2],
                  excess9m_regression %>% lm(data=.,year_30 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% coef() %>% .[2])
excess12m_regression %>% lm(data=.,year_20 ~ .lnEPU +Level +Slope+ Curvature+Log_DY+Log_IP + Log_CPI) %>% summary()