The code are as below
#########################################################################################
# Step 1 - Loading related data
#########################################################################################
install.packages("lubridate")
install.packages("RPostgres")
install.packages("zoo")
install.packages("broom")
install.packages("readr")
install.packages("data.table")
uninstall.packages("dplyr")
install.packages("moments")
install.packages("tidyquant")
install.packages("timetk")
install.packages("rlang")
install.packages("tidyverse")
install.packages("dplyr")
install.packages("installr")
detach(package:tidyverse, unload = TRUE)
packageDescription("dplyr", fields = "Depends")
updateR()
######################
## start
######################
library(tidyverse)
library(rlang)
library(tidyquant) # To download the data
library(plotly) # To create interactive charts
library(timetk) # To manipulate the data series
library(tidyr)
library(dplyr)
# Related libraries
library(tidyverse); library(tidyr); library(zoo);
library(lubridate); library(readr); library(moments);
library(RPostgres); library(scales); library(broom);
require(data.table)
portfolio <- "FF"
# calculate BE
AUT.BE <- read.csv("AUT_yly.csv")
AUT.BE <- AUT.BE %>%
filter(!is.na(cequ))
AUT.BE[,5][is.na(AUT.BE[,5])] <- 0
AUT.BE <- AUT.BE %>%
mutate(BE = cequ + dtax)
# create portfolio
AUT.ret <- read.csv("AUT_mly.csv")
AUT.ret <- AUT.ret %>%
select(-c(ret,mv))
AUT.ret <- AUT.ret %>%
filter(ret_usd<= 890) %>%
filter(!is.na(ret_usd))
AUT.ret <- AUT.ret %>%
mutate(YEAR = year(Date))
AUT.ret <- AUT.ret %>%
left_join(AUT.BE, by=c("YEAR","Id"))
AUT.ret <- AUT.ret %>%
select(-?..country.y)
AUT.ret <- AUT.ret %>%
mutate(BEME = BE/mv_usd)
AUT.ret <- AUT.ret %>%
mutate(ym = as.yearmon(Date))
AUT_quantiles <- AUT.ret %>%
filter(month(ym) == 6) %>%
group_by(ym) %>%
summarize(BSmed = median(mv_usd, na.rm = TRUE),
HNLq30 = quantile(BEME, probs = 0.30, na.rm = TRUE),
HNLq70 = quantile(BEME, probs = 0.70, na.rm = TRUE)) %>%
mutate(ym = as.Date(ym) %m+% months(1)) %>%
mutate(ym = as.yearmon(ym))
AUT.ret.QUANTILE <- AUT.ret %>%
left_join(AUT_quantiles, by=c("ym")) %>%
arrange(ym) %>%
mutate_at(vars(BSmed,HNLq30,HNLq70), funs(na.locf(.,na.rm = FALSE)))
AUT.ret.QUANTILE <- AUT.ret.QUANTILE %>%
mutate(BS = case_when( mv_usd < BSmed ~ "S",
mv_usd >= BSmed ~ "B"),
HNL = case_when( BEME < HNLq30 ~ "L",
BEME >= HNLq30 & BEME < HNLq70 ~ "N",
BEME >= HNLq70 ~ "H"))
AUT.ret.QUANTILE <- AUT.ret.QUANTILE %>%
mutate(LABEL = paste0(BS,HNL))
AUT.ret.QUANTILE <- AUT.ret.QUANTILE %>%
select(-c(BS,HNL))
AUT.ret.QUANTILE <- AUT.ret.QUANTILE %>%
filter(!is.na(ret_usd))
setDT(AUT.ret.QUANTILE)
AUT_mly_FF3 <- AUT.ret.QUANTILE %>%
group_by(LABEL, ym) %>%
summarize(ret_VW = weighted.mean(ret_usd, mv_usd, na.rm = TRUE),
ret_EW = mean(ret_usd, na.rm = TRUE)) %>%
ungroup() %>% as.data.table() %>%
melt(id.vars = c("ym", "LABEL")) %>%
reshape(idvar= c("ym", "variable"),
timevar = "LABEL",
direction ="wide") %>%
mutate(SMB = ((value.SH + value.SN + value.SL)/3) - ((value.BH + value.BN + value.BL)/3),
HML = ((value.SH + value.BH)/2) - ((value.SL + value.BL)/2)) %>%
select(ym,variable,SMB,HML)
AUT.retANALYZE <- AUT.ret %>%
filter(!is.na(ret_usd)) %>%
left_join(AUT_mly_FF3, by = "ym") %>%
filter(variable == "ret_EW")
REG_MODEL <- lm("ret_usd ~ SMB + HML", data = AUT.retANALYZE)
############################ start again
## load factor
factors_ff = read.csv("Europe_5_Factors.csv")
factors_subset_month = factors_ff %>% filter(yearmonth >= 200001) %>% filter(yearmonth < 201800)
factors_subset_month
rf = (factors_subset_month %>% colMeans())['RF']/100
rf
#%>% mean()
#####################################
AUT.BE <- read.csv("AUT_yly.csv")
AUT.BE <- AUT.BE %>%
filter(!is.na(cequ))
AUT.BE[,5][is.na(AUT.BE[,5])] <- 0
AUT.BE <- AUT.BE %>%
mutate(BE = cequ + dtax) %>% mutate(OLD_YEAR = YEAR , YEAR = YEAR +1 )
AUT.retANALYZE2 <- AUT.retANALYZE %>% select(Id , Date , ret_usd , up , mv_usd, YEAR , ym , variable , SMB, HML)
Aut.joined <- AUT.retANALYZE2 %>% left_join(AUT.BE , by =c("YEAR", "Id") )
list_11 <- Aut.joined %>% filter(YEAR >= 2000& YEAR < 2018) %>% filter(is.na(BE) | is.na(ret_usd)) %>% select(Id) %>% unique()
df= Aut.joined %>% anti_join(list_11,by="Id")
df %>% filter(YEAR >= 2000 & YEAR < 2018) %>% count()
df3 <- df %>% filter(YEAR >= 2000& YEAR < 2018) %>% group_by(Id) %>% summarize(cnt = n())
view(df3)
n_cnt = ((2017-2000)+1) *12
n_cnt
# view(cnt)
list_stock2 = df3 %>% filter(cnt == n_cnt) %>% select(Id) %>% unique()
df_joined <- df %>% inner_join(list_stock2 , by =c("Id") ) %>% filter(YEAR >= 2000& YEAR < 2018)
## number = 216
df_joined %>% select(Id) %>% unique() %>% count()
log_ret_tidy <- df_joined %>% mutate(ret_usd = ret_usd/100) %>%
select(Id , ym,ret_usd)
log_ret_xts <- log_ret_tidy %>%
spread(Id, value = ret_usd) %>%
tk_xts()
log_ret_xts
mean_ret <- colMeans(log_ret_xts)
mean_ret
cov_mat <- cov(log_ret_xts) * 12
cov_mat
print(round(cov_mat,4))
wts = runif(n = length(mean_ret ))
sum(wts)
wts <- wts/sum(wts)
wts
sum(wts)
port_returns <- (sum(wts * mean_ret) + 1)^12 - 1
port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts))
print(port_returns)
print(port_risk)
# eq weight
weight = 1/length(mean_ret )
wts = rep(weight, length(mean_ret )) %>% as.matrix()
length(mean_ret )
mat = log_ret_xts %>% as.matrix()
retport_eqw = mat %*% wts
## MV weight
log_ret_tidy_mv <- df_joined %>% mutate(ret_usd = ret_usd/100) %>%
select(Id , ym, mv_usd)
log_ret_xts_mv <- log_ret_tidy_mv %>%
spread(Id, value = mv_usd) %>%
tk_xts()
dfsum <- log_ret_xts_mv %>% as.data.frame() %>% rowSums() %>% as.data.frame()
names(dfsum) <- "mv"
log_ret_xts_mv %>% as.data.frame()
inv_mv_df = dfsum %>% mutate(inv_mv = 1/mv) %>% select(inv_mv)
x = log_ret_xts_mv %>% as.matrix()
cbind(x[1], prop.table(as.matrix(x[-1])), margin = 1)
weighted_x <- x/rowSums(x)
names(weighted_x %>% as.data.frame())
names(mat %>% as.data.frame())
rowSums(weighted_x)
## st
dim(mat)
dim(weighted_x)
mv_mat <- rowSums(mat * weighted_x)
retport_MV <- mv_mat
retport_MV
### book value
names(df_joined)
log_ret_tidy_be <- df_joined %>% mutate(ret_usd = ret_usd/100) %>%
select(Id , ym, BE)
log_ret_xts_be <- log_ret_tidy_be %>%
spread(Id, value = BE) %>%
tk_xts()
dfsum <- log_ret_xts_be %>% as.data.frame() %>% rowSums() %>% as.data.frame()
names(dfsum) <- "be"
log_ret_xts_be
x = log_ret_xts_be %>% as.matrix()
cbind(x[1], prop.table(as.matrix(x[-1])), margin = 1)
weighted_x <- x/rowSums(x)
names(weighted_x %>% as.data.frame())
names(mat %>% as.data.frame())
rowSums(weighted_x)
## st
dim(mat)
dim(weighted_x)
be_mat <- rowSums(mat * weighted_x)
retport_BE <- be_mat
## mean
retport_eqw %>% mean()
retport_MV %>% mean()
retport_BE%>% mean()
## sd
retport_eqw %>% sd()
retport_MV %>% sd()
retport_BE%>% sd()
##########################################
#################### start replacement : run install package only 1st time
##########################################
#install.packages("timeSeries")
#install.packages("PortfolioAnalytics")
#install.packages("ROI")
#install.packages(c("ROI.plugin.glpk", "ROI.plugin.quadprog" , "ROI.plugin.symphony"))
#######################################
library(timeSeries)
library(PortfolioAnalytics)
library(PerformanceAnalytics)
library(ROI)
library(zoo)
library(ROI.plugin.glpk)
library(ROI.plugin.quadprog)
library(ROI.plugin.symphony)
# convert zoo into timeSeries
index(log_ret_xts) <- as.Date(index(log_ret_xts))
data_ret <- as.timeSeries(log_ret_xts)
data_ret
assets <- colnames(log_ret_xts)
portfolio.init <- portfolio.spec(assets)
# portfolio.init <- add.constraint(portfolio.init, type = "full_investment")
portfolio.init <- add.constraint(portfolio.init, type="long_only")
portfolio.init <- add.constraint(portfolio.init, type="box",min=0, max=0.1)
# calculate minimum std. dev. portfolio
portfolio.minSD <- add.objective(portfolio = portfolio.init, type="risk", name="StdDev")
portfolio.minSD.opt <- optimize.portfolio(data_ret, portfolio = portfolio.minSD, optimize_method = "ROI", trace = TRUE)
portfolio.minSD.weights <- portfolio.minSD.opt$weights
portfolio.minSD.weights
##excess ret
excess_ret_data = data_ret - rf
## max sharp
sharpe.portf <- add.objective(portfolio=portfolio.init, type="risk", name="StdDev")
sharpe.portf <- add.objective(portfolio=sharpe.portf, type="return", name="mean")
# Optimization to maximize Sharpe Ratio
max_sharpe_opt <- optimize.portfolio(excess_ret_data, portfolio=sharpe.portf, optimize_method="ROI", maxSR=TRUE)
max_sharpe_opt
max_sharpe_opt.weights <- max_sharpe_opt$weights
max_sharpe_opt.weights
##########################################
### end replacement
##########################################
#min_var_weight
#max_sharp_weight
# w_min = t(portfolio.minSD.weights %>% as.matrix())
retport_minvar = mat %*% portfolio.minSD.weights
retport_minvar
# w_max = t(max_sharp_weight %>% as.matrix())
retport_maxshp = mat %*% max_sharpe_opt.weights
retport_maxshp
retport_minvar %>% mean()
retport_minvar %>% sd()
retport_maxshp %>% mean()
retport_maxshp %>% sd()
### gather all return
retport_eqw %>% mean()
retport_MV %>% as.matrix() %>% mean()
retport_BE %>% as.matrix() %>% mean()
retport_maxshp %>% mean()
retport_minvar %>% mean()
retport_eqw %>% sd()
retport_MV %>% as.matrix() %>% sd()
retport_BE %>% as.matrix() %>% sd()
retport_maxshp %>% sd()
retport_minvar %>% sd()
retport_eqw %>% str()
retport_MV %>% str()
retport_BE%>% str()
retport_maxshp%>% str()
retport_minvar%>% str()
retport_MV = retport_MV %>% as.data.frame()
retport_BE = retport_BE %>% as.data.frame()
ret_combine = cbind(rownames(retport_eqw) , retport_eqw , retport_MV , retport_BE , retport_maxshp , retport_minvar)
names(ret_combine) <- c("ym","ret_equal_weight","ret_market_value" , "ret_book_value" , "ret_max_sharpe" ,"ret_min_var")
### all return and sd
ret_combine[c("ret_equal_weight","ret_market_value" , "ret_book_value" , "ret_max_sharpe" ,"ret_min_var")] %>% colSums()
apply(ret_combine[c("ret_equal_weight","ret_market_value" , "ret_book_value" , "ret_max_sharpe" ,"ret_min_var")], 2, sd)
df_factor_1 <- df_joined %>% select(ym,SMB ,HML)
df_factor_2 <- factors_subset_month %>% select(yearmonth , Mkt.RF , RMW , CMA, RF)
df_combined_premod = cbind(ret_combine , df_factor_1 , df_factor_2)
tail(df_combined_premod)
names(df_combined_premod)
###################################
## new apply
###############################
factor_mom <- read.csv("factor_mom.csv")
factor_mom <- factor_mom %>% mutate( ym = str_pad( ym , width =6 , side= c("left") , pad = "0" )) %>% select(ym, return)
names(factor_mom) <- c("yearmonth" , "WML")
names(factor_mom)
names(df_combined_premod) <- c("ym1", "ret_equal_weight" , "ret_market_value" , "ret_book_value" , "ret_max_sharpe","ret_min_var", "ym", "SMB", "HML", "yearmonth","Mkt.RF","RMW","CMA","RF")
factor_mom = transform(factor_mom, yearmonth = as.numeric(yearmonth))
df_combined_premod2 <- df_combined_premod %>%
left_join(factor_mom, by=c("yearmonth"))
ff4_equal_weight <- lm(ret_equal_weight ~ Mkt.RF+ SMB + HML + WML , data = df_combined_premod2)
ff4_market_value <- lm(ret_market_value ~ Mkt.RF + SMB + HML + WML , data = df_combined_premod2)
ff4_book_value <- lm(ret_book_value ~ Mkt.RF + SMB + HML + WML , data = df_combined_premod2)
ff4_max_sharpe <- lm(ret_max_sharpe ~ Mkt.RF + SMB + HML + WML , data = df_combined_premod2)
ff4_min_var <- lm(ret_min_var ~ Mkt.RF + SMB + HML + WML , data = df_combined_premod2)
### modeling
# CAPM with alpha
capm_equal_weight <- lm(ret_equal_weight ~ Mkt.RF , data = df_combined_premod)
capm_market_value <- lm(ret_market_value ~ Mkt.RF , data = df_combined_premod)
capm_book_value <- lm(ret_book_value ~ Mkt.RF , data = df_combined_premod)
capm_max_sharpe <- lm(ret_max_sharpe ~ Mkt.RF , data = df_combined_premod)
capm_min_var <- lm(ret_min_var ~ Mkt.RF , data = df_combined_premod)
# summary
summary(capm_equal_weight)
summary(capm_market_value)
summary(capm_book_value)
summary(capm_max_sharpe)
summary(capm_min_var)
## wrrite csv
write.csv(tidy(capm_equal_weight), "capm_equal_weight.csv")
write.csv(tidy(capm_market_value), "capm_market_value.csv")
write.csv(tidy(capm_book_value), "capm_book_value.csv")
write.csv(tidy(capm_max_sharpe), "capm_max_sharpe.csv")
write.csv(tidy(capm_min_var), "capm_min_var.csv")
## FF 4 factor model
ff4_equal_weight <- lm(ret_equal_weight ~ Mkt.RF+ SMB + HML + RMW , data = df_combined_premod)
ff4_market_value <- lm(ret_market_value ~ Mkt.RF + SMB + HML + RMW , data = df_combined_premod)
ff4_book_value <- lm(ret_book_value ~ Mkt.RF + SMB + HML + RMW , data = df_combined_premod)
ff4_max_sharpe <- lm(ret_max_sharpe ~ Mkt.RF + SMB + HML + RMW , data = df_combined_premod)
ff4_min_var <- lm(ret_min_var ~ Mkt.RF + SMB + HML + RMW , data = df_combined_premod)
summary(ff4_equal_weight)
summary(ff4_market_value)
summary(ff4_book_value)
summary(ff4_max_sharpe)
summary(ff4_min_var)
write.csv(tidy(ff4_equal_weight), "ff4_equal_weight.csv")
write.csv(tidy(ff4_market_value), "ff4_market_value.csv")
write.csv(tidy(ff4_book_value), "ff4_book_value.csv")
write.csv(tidy(ff4_max_sharpe), "ff4_max_sharpe.csv")
write.csv(tidy(ff4_min_var), "ff4_min_var.csv")
### FF 5 factor model
ff_equal_weight <- lm(ret_equal_weight ~ Mkt.RF+ SMB + HML + RMW + CMA , data = df_combined_premod)
ff_market_value <- lm(ret_market_value ~ Mkt.RF + SMB + HML + RMW + CMA, data = df_combined_premod)
ff_book_value <- lm(ret_book_value ~ Mkt.RF + SMB + HML + RMW + CMA, data = df_combined_premod)
ff_max_sharpe <- lm(ret_max_sharpe ~ Mkt.RF + SMB + HML + RMW + CMA, data = df_combined_premod)
ff_min_var <- lm(ret_min_var ~ Mkt.RF + SMB + HML + RMW + CMA, data = df_combined_premod)
## summary
summary(ff_equal_weight)
summary(ff_market_value)
summary(ff_book_value)
summary(ff_max_sharpe)
summary(ff_min_var)
write.csv(tidy(ff_equal_weight), "ff_equal_weight.csv")
write.csv(tidy(ff_market_value), "ff_market_value.csv")
write.csv(tidy(ff_book_value), "ff_book_value.csv")
write.csv(tidy(ff_max_sharpe), "ff_max_sharpe.csv")
write.csv(tidy(ff_min_var), "ff_min_var.csv")




Leave a comment