library(tidyverse)
library(lubridate)
library(gridExtra)
# library(vars)
load(file = "data/TPU.rdata")
# for loop
commodity.names = c("豆二")
# generate commodity.names[i]_RV.csv data file in the data folder
for (i in 1:length(commodity.names)){
read_csv(paste0("data/",commodity.names[i],".csv")) %>%
as_tibble() %>%
rename(date = '交易日期',prc = '收盘价') %>%
mutate(date = as.Date(date)) %>%
arrange(date) %>%
mutate(prc = as.numeric(prc)) %>%
mutate(ret_ = prc/lag(prc)-1) %>%
dplyr::select(date,prc,ret_) %>%
mutate(year = year(date),quarter = quarter(date)) %>%
filter(!is.na(prc)) %>%
group_by(year,quarter) %>%
summarise("RV" = 252*var(ret_)) %>%
write_csv(file = paste0("data\\",commodity.names[i],"_RV.csv"))
read_csv(paste0("data/",commodity.names[i],".csv")) %>%
as_tibble() %>%
rename(date = '交易日期',prc = '收盘价') %>%
mutate(date = as.Date(date)) %>%
arrange(date) %>%
mutate(prc = as.numeric(prc)) %>%
mutate(year = year(date),quarter = quarter(date)) %>%
filter(!is.na(prc)) %>%
group_by(year,quarter) %>%
filter(date == last(date)) %>%
ungroup() %>%
mutate(return = prc/lag(prc)-1) %>%
mutate(LR = log(1+return)) %>%
dplyr::select(year,quarter,LR) %>%
write_csv(file = paste0("data\\",commodity.names[i],"_LR.csv"))
}
data_gpr_export <- read_csv("data/data_gpr_export.csv") # 1461 rows
OECD_plus6_industrial_production <- read_csv("data/OECD_plus6_industrial_production.csv")
TPU_china <- read_csv("data/TPU_china.csv")
TPU_US <- read_csv("data/TPU_US.csv")
data_gpr_export %>%
dplyr::select(month,GPR) %>%
dplyr::filter(!is.na(GPR)) %>%
mutate(date = as.Date(month)) %>%
mutate(year = year(date),quarter = quarter(date)) %>%
dplyr::select(date,year,quarter,GPR) %>%
group_by(year,quarter) %>%
summarise(quarter_GPR = mean(GPR)) %>%
rename(GPR = quarter_GPR) -> GPR
TPU_china %>%
dplyr::select(year,month,TPU) %>%
rename(TPU_china = TPU) %>%
mutate(date = as.Date(paste0(year,"-",month,"-1"))) %>%
mutate(quarter = quarter(date)) %>%
group_by(year,quarter) %>%
mutate(mean_TPU_china = mean(TPU_china)) %>%
ungroup() %>%
dplyr::select(year,quarter,mean_TPU_china) %>%
rename(TPU_china = mean_TPU_china) %>%
distinct() -> TPUchina
TPU_US %>%
dplyr::select(Date,`9. Trade policy`) %>%
rename(TPU_US = `9. Trade policy`) %>%
mutate(Date = as.Date(Date)) %>%
mutate(year = year(Date),quarter = quarter(Date)) %>%
group_by(year,quarter) %>%
mutate(mean_TPU_US = mean(TPU_US)) %>%
dplyr::select(year,quarter,mean_TPU_US) %>%
rename(TPU_US = mean_TPU_US) %>%
distinct() -> TPUUS
lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C")
## [1] "C"
zoo::as.zoo(tpu_cimpr) %>%
zoo::index() %>%
as.character() %>%
cbind(as.matrix(zoo::as.zoo(tpu_cimpr))) %>%
as.data.frame() %>%
rename(date = ".",tpu_cimpr = "zoo::as.zoo(tpu_cimpr)") %>%
as_tibble() %>%
mutate(date = paste0(date," 1")) %>%
mutate(date = as.Date(date,tryFormats = "%b %Y %d")) %>%
mutate(year = year(date),quarter = quarter(date)) %>%
mutate(tpu_cimpr = as.numeric(tpu_cimpr)) %>%
group_by(year,quarter) %>%
mutate(mean_tpu_cimpr = mean(tpu_cimpr)) %>%
dplyr::select(year,quarter,mean_tpu_cimpr) %>%
rename(tpu_cimpr = mean_tpu_cimpr) %>%
distinct() %>%
dplyr::select(year,quarter,tpu_cimpr) -> tpu_cimpr
# read_csv("data/gdp.csv")
read_csv("data/gdp.csv") %>%
mutate(date = as.Date(date)) %>%
mutate(year = year(date),quarter = quarter(date)) %>%
rename(log_gdp_index = 'log gdp index') %>%
rename(gdp_qtr_to_qtr = '国内生产总值环比增长速度(%)') %>%
mutate(gdp_qtr_to_qtr = as.numeric(gdp_qtr_to_qtr)) %>%
dplyr::select(year,quarter,log_gdp_index,gdp_qtr_to_qtr) -> GDP
## Rows: 44 Columns: 4
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (2): date, 国内生产总值环比增长速度(%)
## dbl (2): gdp index, log gdp index
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Warning in mask$eval_all_mutate(quo): 强制改变过程中产生了NA
GDP
## # A tibble: 44 x 4
## year quarter log_gdp_index gdp_qtr_to_qtr
## <dbl> <int> <dbl> <dbl>
## 1 2010 4 0 NA
## 2 2011 1 0.0107 2.5
## 3 2011 2 0.0202 2.2
## 4 2011 3 0.0283 1.9
## 5 2011 4 0.0348 1.5
## 6 2012 1 0.0421 1.7
## 7 2012 2 0.0520 2.3
## 8 2012 3 0.0602 1.9
## 9 2012 4 0.0684 1.9
## 10 2013 1 0.0765 1.9
## # ... with 34 more rows
# GDP
# merge data
c("豆二")
## [1] "豆二"
read_csv(paste0("data\\","豆二","_LR.csv")) %>%
rename(LR_豆二 = LR) %>%
full_join(read_csv(paste0("data\\","豆二","_RV.csv"))) %>%
rename(RV_豆二 = RV) %>%
full_join(GPR) %>%
full_join(GDP) %>%
full_join(TPUchina) %>%
full_join(TPUUS) %>%
full_join(tpu_cimpr) %>%
mutate(LGPR = log(GPR)) %>%
dplyr::filter(year>=2011) %>%
mutate(Log_TPU_china = log(1+TPU_china),Log_TPU_US = log(1+TPU_US),Log_TPU_cimpr = log(1+tpu_cimpr)) %>%
arrange(year,quarter) -> data
ggplot_IRF <- function(data = data,commodity_type = commodity_type,ortho = ...,gdp = ...){
# test if the model is stationary
if(gdp == "log_gdp_index"){
data[c("year","quarter","LGPR","log_gdp_index","TPU_china","TPU_US","tpu_cimpr",paste0("LR_",commodity_type),paste0("RV_",commodity_type))] %>%
rename(gdp = log_gdp_index) ->data
}
else if(gdp == "gdp_qtr_to_qtr"){
data[c("year","quarter","LGPR","gdp_qtr_to_qtr","TPU_china","TPU_US","tpu_cimpr",paste0("LR_",commodity_type),paste0("RV_",commodity_type))] %>%
rename(gdp = gdp_qtr_to_qtr) ->data
}
data %>%
filter(!is.na(LGPR)) %>%
filter(!is.na(gdp)) %>%
filter(!is.na(TPU_china)) %>%
filter(!is.na(TPU_US)) %>%
filter(!is.na(tpu_cimpr)) %>%
filter(year<=2019) -> data
obj1 <- vars::irf(vars::VAR(data[c("gdp",paste0("LR_",commodity_type),paste0("RV_",commodity_type),"LGPR","TPU_china")],p = 1),impulse = "TPU_china", response = paste0("RV_",commodity_type),n.ahead = 8, ortho = ortho,boot = TRUE,runs = 1000)
obj2 <- vars::irf(vars::VAR(data[c("gdp",paste0("LR_",commodity_type),paste0("RV_",commodity_type),"LGPR","TPU_US")],p = 1),impulse = "TPU_US", response = paste0("RV_",commodity_type),n.ahead = 8, ortho = ortho,boot = TRUE,runs = 1000)
obj3 <- vars::irf(vars::VAR(data[c("gdp",paste0("LR_",commodity_type),paste0("RV_",commodity_type),"LGPR","tpu_cimpr")],p = 1),impulse = "tpu_cimpr", response = paste0("RV_",commodity_type),n.ahead = 8, ortho = ortho,boot = TRUE,runs = 1000)
obj1$irf %>%
as.data.frame() %>%
cbind(Lower = as.data.frame(obj1$Lower)) %>%
cbind(Upper = as.data.frame(obj1$Upper)) %>%
cbind(period = 0:8) -> df1
colnames(df1) <- c("estimate","lower","upper","period")
df1 %>%
ggplot()+
geom_line(aes(x = period,y = estimate),colour = "blue")+
geom_line(aes(x = period,y = lower),linetype = 2,colour = "red")+
geom_line(aes(x = period,y = upper),linetype = 2,colour = "red")+
ggtitle(paste0("Impulse:",obj1$impulse))+
theme(plot.title = element_text(hjust = 0.5))+
theme(legend.position = "bottom")+
ylab(paste0("RV_",commodity_type))+
theme_classic()+
geom_hline(yintercept = 0) -> p1
obj2$irf %>%
as.data.frame() %>%
cbind(Lower = as.data.frame(obj2$Lower)) %>%
cbind(Upper = as.data.frame(obj2$Upper)) %>%
cbind(period = 0:8) -> df2
colnames(df2) <- c("estimate","lower","upper","period")
df2 %>%
ggplot()+
geom_line(aes(x = period,y = estimate),colour = "blue")+
geom_line(aes(x = period,y = lower),linetype = 2,colour = "red")+
geom_line(aes(x = period,y = upper),linetype = 2,colour = "red")+
ggtitle(paste0("Impulse:",obj2$impulse),subtitle = paste0("ortho=",ortho))+
theme(plot.title = element_text(hjust = 0.5))+
# theme(legend.position = "bottom")+
ylab(paste0("RV_",commodity_type))+
theme_classic()+
geom_hline(yintercept = 0) -> p2
obj3$irf %>%
as.data.frame() %>%
cbind(Lower = as.data.frame(obj3$Lower)) %>%
cbind(Upper = as.data.frame(obj3$Upper)) %>%
cbind(period = 0:8) -> df3
colnames(df3) <- c("estimate","lower","upper","period")
df3 %>%
ggplot()+
geom_line(aes(x = period,y = estimate),colour = "blue")+
geom_line(aes(x = period,y = lower),linetype = 2,colour = "red")+
geom_line(aes(x = period,y = upper),linetype = 2,colour = "red")+
ggtitle(paste0("Impulse:",obj3$impulse))+
theme(plot.title = element_text(hjust = 0.5))+
theme(legend.position = "bottom")+
ylab(paste0("RV_",commodity_type))+
theme_classic()+
geom_hline(yintercept = 0) -> p3
return(gridExtra::grid.arrange(p1,p2,p3, ncol = 3))
}
# baseline
# ortho
ggplot_IRF(data = data,commodity_type = "豆二",ortho = F,gdp = "gdp_qtr_to_qtr")
ggplot_IRF(data = data,commodity_type = "豆二",ortho = T,gdp = "gdp_qtr_to_qtr")
ggplot_IRF(data = data,commodity_type = "豆二",ortho = F,gdp = "log_gdp_index")
ggplot_IRF(data = data,commodity_type = "豆二",ortho = T,gdp = "log_gdp_index")
data %>%
filter(year<=2017) -> data
ggplot_IRF(data = data,commodity_type = "豆二",ortho = F,gdp = "gdp_qtr_to_qtr")
ggplot_IRF(data = data,commodity_type = "豆二",ortho = T,gdp = "gdp_qtr_to_qtr")
ggplot_IRF(data = data,commodity_type = "豆二",ortho = F,gdp = "log_gdp_index")
ggplot_IRF(data = data,commodity_type = "豆二",ortho = T,gdp = "log_gdp_index")
