RiskPremium

Measuring the market risk premium
Log | Files | Refs

rp_measure.R (4287B)

      1 #!/usr/bin/env Rscript  
      2 #
      3 # rp_measure.R
      4 #
      5 # This code runs the predictive regression
      6 # We save in output the expected excess return estimate
      7 # 
      8 #
      9 # (c) Valentin Haddad, Erik Loualiche & Matthew Plosser
     10 #
     11 # Last updated on June 4th 2019
     12 # 
     13 ##################################################################################
     14 
     15 ##################################################################################
     16 message("Log file for code executed at\n")
     17 message(format(Sys.time(), "%a %b %d %X %Y"))
     18 ##################################################################################
     19 
     20 
     21 ##################################################################################
     22 # APPEND REQUIRED PACKAGES
     23 
     24 # See this https://stackoverflow.com/questions/4090169/elegant-way-to-check-for-missing-packages-and-install-them
     25 using<-function(...) {
     26     libs<-unlist(list(...))
     27     req<-unlist(lapply(libs,require,character.only=TRUE))
     28     need<-libs[req==FALSE]
     29     if(length(need)>0){ 
     30         install.packages(need)
     31         lapply(need,require,character.only=TRUE)
     32     }
     33 }
     34 
     35 package_to_load <- c("crayon", "devtools", "wesanderson", "ggplot2", "statar", 
     36 	"stringr", "lubridate", "lmtest", "sandwich", "stargazer", "data.table")
     37 
     38 using(package_to_load)
     39 
     40 
     41 check_file = file.exists("log/R-session-info.log.R")
     42 sink("log/R-session-info.log.R", append=check_file)
     43 cat(bold("\n\n# -----\n# Session info for rp_measure.csv\n\n")) 
     44 session_info()
     45 sink()
     46 ##################################################################################
     47 
     48 
     49 ##################################################################################
     50 dt_predict <- fread("./tmp/predict.csv")
     51 dt_predict[, datem := as.monthly(ISOdate(str_sub(dateym, 1, 4), str_sub(dateym, 5, 6), 1)) ]
     52 ##################################################################################
     53 
     54 
     55 ##################################################################################
     56 r1_3 = lm(rmrf_y3 ~ dp + cay + rf, data = dt_predict[ year(datem) < 2011])
     57 nw_r1_3 = coeftest(r1_3, df = Inf, vcov = NeweyWest(r1_3, lag = 12, prewhite = FALSE) )
     58 r2_3 = lm(rmrf_y3 ~ dp + cay + rf, data = dt_predict)
     59 nw_r2_3 = coeftest(r2_3, df = Inf, vcov = NeweyWest(r2_3, lag = 12, prewhite = FALSE) )
     60 
     61 stargazer(r1_3, nw_r1_3, r2_3, nw_r2_3, type="text")
     62 
     63 # --- output for the readme
     64 star = stargazer(r2_3, type="text",  style = "aer",
     65 	covariate.labels = c("D/P ratio", "cay", "T-bill (three-month)"),
     66 	dep.var.labels   = "Future Excess Returns",
     67 	omit.stat = c("ser", "adj.rsq") )
     68 
     69 star[1] = "~~~R"
     70 star[length(star)+1] = "~~~"
     71 cat(star, sep = '\n', file = './tmp/reg_update.txt')
     72 ##################################################################################
     73 
     74 
     75 ##################################################################################
     76 # OUTPUT PREDICTED VALUE
     77 dt_exp_rmrf <- cbind(dt_predict[!is.na(rmrf_y3), -c("datem")], exp_rmrf = predict(r2_3))
     78 
     79 fwrite(dt_exp_rmrf, "./output/predict.csv")
     80 ##################################################################################
     81 
     82 
     83 ##################################################################################
     84 # PLOT
     85 dt_plot <- dt_exp_rmrf[, .(
     86 	date=as.Date(ISOdate(str_sub(dateym,1, 4), as.integer(str_sub(dateym, 5, 6)), 1)), 
     87 	dp, cay, rf, rmrf_y3, exp_rmrf)]
     88 dt_plot[]
     89 
     90 
     91 p0 <- dt_plot[, .(date, dp, cay, rf, rmrf_y3) ] %>% 
     92     melt(id.vars="date") %>%
     93 	ggplot(aes(date, value, colour = variable)) + 
     94 	geom_line(alpha=0.75, size=0.25) + geom_point(shape=1, size = 1, alpha=0.5) + 
     95 	theme_bw()
     96 # p0
     97 
     98 p1 <- dt_plot[, .(date, exp_rmrf, rmrf_y3) ] %>% 
     99  	melt(id.vars="date") %>%
    100 	ggplot(aes(date, 100*value, colour = variable)) + 
    101 	geom_line(alpha=0.75, size=0.25) + geom_point(shape=1, size = 1, alpha=0.5) + 
    102 	xlab("") + ylab("Returns (percent)") + 
    103 	theme_bw() +
    104 	theme(legend.position = c(0.3, 0.9)) + 
    105 	scale_colour_manual(name  = "",
    106                         breaks = c("exp_rmrf", "rmrf_y3"),
    107                         values = c(wes_palette("Zissou1")[1], wes_palette("Zissou1")[5]),
    108                         labels=c("Expected", "Realized")) + 
    109 	guides(colour = guide_legend(nrow = 1))
    110 ggsave("./output/predict.png", p1, width = 8, height=6)
    111 ##################################################################################
    112 
    113 
    114 ##################################################################################