library(tidyquant) library(timetk) library(lubridate) library(dplyr) library(tidyr) library(purrr)
tickers <- c(“SPY”, “QQQ”, “EEM”, “IWM”, “EFA”, “TLT”, “IYR”, “GLD”)
etf_data <- tq_get(tickers, from = “2010-01-01”, to = Sys.Date(), get = “stock.prices”)
etf_prices <- etf_data %>% select(date, symbol, adjusted) %>% pivot_wider(names_from = symbol, values_from = adjusted)
monthly_returns <- etf_prices %>% tk_xts(date_var = date) %>% Return.calculate(method = “simple”) %>% to.monthly(indexAt = “lastof”, OHLC = FALSE) %>% tk_tbl(preserve_index = TRUE, rename_index = “date”) %>% na.omit()
weekly_returns <- etf_prices %>% tk_xts(date_var = date) %>% to.weekly(indexAt = “lastof”, OHLC = FALSE) %>% Return.calculate(method = “simple”) %>% tk_tbl(preserve_index = TRUE, rename_index = “date”) %>% na.omit()
tk_tbl()monthly_returns_tbl <- monthly_returns
ff_factors <- read.csv(“F-F_Research_Data_Factors.csv”, skip = 3) %>% rename(date = X) %>% mutate(date = ymd(paste0(date, “01”))) %>% filter(!is.na(Mkt.RF)) %>% mutate(across(Mkt.RF:HML, ~ . / 100)) # Convert from percentage
merged_data <- left_join(monthly_returns_tbl, ff_factors, by = “date”)
returns_window <- merged_data %>% filter(date >= ymd(“2010-02-01”) & date <= ymd(“2015-01-01”))
capm_data <- returns_window %>% mutate(across(SPY:GLD, ~ . - RF))
capm_cov_matrix <- capm_data %>% select(SPY:GLD) %>% cov()
ff_cov_matrix <- returns_window %>% mutate(across(SPY:GLD, ~ . - RF)) %>% select(SPY:GLD) %>% cov()
library(quadprog)
cov_mat <- capm_cov_matrix n_assets <- ncol(cov_mat)
Dmat <- cov_mat dvec <- rep(0, n_assets) Amat <- cbind(rep(1, n_assets)) # weights sum to 1 bvec <- 1 meq <- 1
gmv_result <- solve.QP(Dmat, dvec, Amat, bvec, meq) gmv_weights <- gmv_result$solution names(gmv_weights) <- colnames(cov_mat)
library(PerformanceAnalytics)
backtest_gmv <- function(returns, window = 60) { n <- nrow(returns) weights_list <- list() gmv_returns <- numeric(n - window)
for (i in 1:(n - window)) { ret_window <- returns[i:(i + window - 1), ] cov_mat <- cov(ret_window) sol <- solve.QP(Dmat = cov_mat, dvec = rep(0, ncol(cov_mat)), Amat = cbind(rep(1, ncol(cov_mat))), bvec = 1, meq = 1) weights <- sol$solution next_ret <- as.numeric(returns[i + window, ]) %*% weights gmv_returns[i] <- next_ret }
return(xts(gmv_returns, order.by = index(returns)[(window + 1):n])) }
monthly_xts <- monthly_returns_tbl %>% select(date, SPY:GLD) %>% tk_xts(date_var = date)
gmv_capm_returns <- backtest_gmv(monthly_xts)
chart.CumReturns(gmv_capm_returns, main = “GMV Portfolio (CAPM Model)”, wealth.index = TRUE)