# Load libraries
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(tidyquant)
## Loading required package: lubridate
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
## Loading required package: PerformanceAnalytics
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
library(lubridate)
library(timetk)
library(purrr)
library(tibble)
library(PerformanceAnalytics)
library(magrittr)
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
##
## set_names
library(readr)
library(xts)
library(dplyr)
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
##
## first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ tidyr 1.3.1
## ✔ stringr 1.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks magrittr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks xts::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks xts::last()
## ✖ magrittr::set_names() masks purrr::set_names()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Define the ETF symbols
symbols <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
# Download the data from Yahoo Finance
prices <- getSymbols(symbols, src = 'yahoo', from = "2010-01-01", auto.assign = TRUE, warnings = FALSE) %>%
map(~ Ad(get(.))) %>%
reduce(merge) %>%
`colnames<-`(symbols)
# Convert xts to tibble
prices_tibble <- tk_tbl(prices, rename_index = "date")
# Ensure date column is in character format, then convert to Date type
prices_tibble <- prices_tibble %>%
mutate(date = as.character(date)) %>%
mutate(date = as.Date(date))
# Calculate weekly log returns
weekly_prices <- prices_tibble %>%
group_by(week = floor_date(date, "week")) %>%
summarize(across(all_of(symbols), last))
weekly_returns <- weekly_prices %>%
mutate(across(all_of(symbols), ~ log(. / lag(.)))) %>%
drop_na()
# Calculate monthly log returns
monthly_prices <- prices_tibble %>%
group_by(month = floor_date(date, "month")) %>%
summarize(across(all_of(symbols), last))
monthly_returns <- monthly_prices %>%
mutate(across(all_of(symbols), ~ log(. / lag(.)))) %>%
drop_na()
# Convert to tibble format
monthly_returns_tibble <- as_tibble(monthly_returns)
# Rename 'month' column to 'date'
monthly_returns_tibble <- monthly_returns_tibble %>%
rename(date = month)
# Load the Fama/French data from URL
FFdata <- read.csv("https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors.CSV", skip = 3)
# Clean and format the data, handling parsing errors
FFdata <- FFdata %>%
mutate(date = parse_date_time(X, orders = "Y%m", quiet = TRUE)) %>%
filter(!is.na(date)) %>%
select(date, Mkt.RF, SMB, HML, RF)
# Convert character columns to numeric
FFdata <- FFdata %>%
mutate(across(-date, ~ as.numeric(gsub(" ", "", .))))
# Ensure date columns are present and in Date format
FFdata <- FFdata %>%
mutate(date = as.Date(date))
# Merge the ETF returns with the FF factors data
merged_data <- left_join(monthly_returns_tibble, FFdata, by = 'date')
# Check the first few rows of the merged data
head(merged_data)
## # A tibble: 6 × 13
## date SPY QQQ EEM IWM EFA TLT IYR GLD
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010-02-01 0.0307 0.0450 0.0176 0.0438 0.00266 -0.00343 0.0531 0.0322
## 2 2010-03-01 0.0591 0.0743 0.0780 0.0791 0.0619 -0.0208 0.0930 -0.00440
## 3 2010-04-01 0.0154 0.0222 -0.00166 0.0552 -0.0284 0.0327 0.0619 0.0572
## 4 2010-05-01 -0.0828 -0.0768 -0.0986 -0.0784 -0.119 0.0498 -0.0585 0.0301
## 5 2010-06-01 -0.0531 -0.0616 -0.0141 -0.0806 -0.0208 0.0564 -0.0478 0.0233
## 6 2010-07-01 0.0661 0.0701 0.104 0.0651 0.110 -0.00951 0.0899 -0.0522
## # ℹ 4 more variables: Mkt.RF <dbl>, SMB <dbl>, HML <dbl>, RF <dbl>
# Extract the relevant returns for the CAPM model
capm_returns <- merged_data %>%
filter(date >= "2010-02-01" & date <= "2015-01-31") %>%
select(-date, -Mkt.RF, -SMB, -HML, -RF)
# Compute the covariance matrix
capm_cov_matrix <- cov(capm_returns)
print(capm_cov_matrix)
## SPY QQQ EEM IWM EFA
## SPY 0.0013789585 0.0014370261 0.0017213335 0.0017634421 0.0015961626
## QQQ 0.0014370261 0.0017632689 0.0017694609 0.0018116951 0.0016582437
## EEM 0.0017213335 0.0017694609 0.0033901742 0.0023127699 0.0024955577
## IWM 0.0017634421 0.0018116951 0.0023127699 0.0026629398 0.0019648695
## EFA 0.0015961626 0.0016582437 0.0024955577 0.0019648695 0.0024253990
## TLT -0.0009663049 -0.0009584967 -0.0011821325 -0.0013080775 -0.0011039475
## IYR 0.0011816473 0.0012041857 0.0017966318 0.0015876886 0.0015471687
## GLD 0.0002043555 0.0003904126 0.0009310469 0.0005394849 0.0004417113
## TLT IYR GLD
## SPY -0.0009663049 0.0011816473 0.0002043555
## QQQ -0.0009584967 0.0012041857 0.0003904126
## EEM -0.0011821325 0.0017966318 0.0009310469
## IWM -0.0013080775 0.0015876886 0.0005394849
## EFA -0.0011039475 0.0015471687 0.0004417113
## TLT 0.0015464559 -0.0003632221 0.0001797673
## IYR -0.0003632221 0.0019945691 0.0005353329
## GLD 0.0001797673 0.0005353329 0.0029025882
# Extract the relevant returns and FF factors for the FF model
ff_returns <- merged_data %>%
filter(date >= "2010-02-01" & date <= "2015-01-31")
# Compute the residuals for each ETF by regressing on the FF factors
residuals <- ff_returns %>%
select(all_of(symbols)) %>%
map_dfc(~ residuals(lm(. ~ Mkt.RF + SMB + HML, data = ff_returns)))
# Compute the covariance matrix of the residuals
ff_cov_matrix <- cov(residuals)
print(ff_cov_matrix)
## SPY QQQ EEM IWM EFA
## SPY 3.259897e-06 3.622633e-06 1.785161e-06 8.352802e-07 1.903229e-06
## QQQ 3.622633e-06 2.067136e-04 -5.274408e-05 -1.248065e-05 -3.610005e-05
## EEM 1.785161e-06 -5.274408e-05 1.213627e-03 4.639690e-05 5.052318e-04
## IWM 8.352802e-07 -1.248065e-05 4.639690e-05 1.656415e-05 3.183300e-05
## EFA 1.903229e-06 -3.610005e-05 5.052318e-04 3.183300e-05 5.327978e-04
## TLT 8.487542e-06 -4.927408e-06 1.964234e-05 2.785540e-05 -2.701238e-05
## IYR 2.136786e-05 -1.638223e-05 3.328146e-04 5.932737e-05 2.082022e-04
## GLD -8.117752e-06 5.140160e-05 5.668786e-04 5.843982e-05 1.954893e-04
## TLT IYR GLD
## SPY 8.487542e-06 2.136786e-05 -8.117752e-06
## QQQ -4.927408e-06 -1.638223e-05 5.140160e-05
## EEM 1.964234e-05 3.328146e-04 5.668786e-04
## IWM 2.785540e-05 5.932737e-05 5.843982e-05
## EFA -2.701238e-05 2.082022e-04 1.954893e-04
## TLT 7.800182e-04 4.560201e-04 2.593838e-04
## IYR 4.560201e-04 1.008805e-03 3.079627e-04
## GLD 2.593838e-04 3.079627e-04 2.510801e-03
# Function to compute GMV portfolio weights
compute_gmv_weights <- function(cov_matrix) {
inv_cov <- solve(cov_matrix)
ones <- rep(1, nrow(inv_cov))
weights <- inv_cov %*% ones / sum(inv_cov %*% ones)
return(weights)
}
# Compute the weights
gmv_weights_capm <- compute_gmv_weights(capm_cov_matrix)
gmv_weights_ff <- compute_gmv_weights(ff_cov_matrix)
print(gmv_weights_capm)
## [,1]
## SPY 0.94117398
## QQQ -0.18088557
## EEM -0.01437269
## IWM -0.05253263
## EFA -0.01342355
## TLT 0.48106080
## IYR -0.21174024
## GLD 0.05071992
print(gmv_weights_ff)
## [,1]
## SPY 0.8299917745
## QQQ 0.0052263248
## EEM 0.0001148441
## IWM 0.1868608954
## EFA 0.0009323804
## TLT 0.0032730421
## IYR -0.0286072359
## GLD 0.0022079745