library(tidyquant)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## ── Attaching core tidyquant packages ──────────────────────── tidyquant 1.0.9 ──
## ✔ PerformanceAnalytics 2.0.8 ✔ TTR 0.24.4
## ✔ quantmod 0.4.26 ✔ xts 0.14.1
## ── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
## ✖ zoo::as.Date() masks base::as.Date()
## ✖ zoo::as.Date.numeric() masks base::as.Date.numeric()
## ✖ PerformanceAnalytics::legend() masks graphics::legend()
## ✖ quantmod::summary() masks base::summary()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(quantmod)
library(lubridate)
##
## Attaching package: 'lubridate'
##
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(timetk)
##
## Attaching package: 'timetk'
##
## The following object is masked from 'package:tidyquant':
##
## FANG
library(purrr)
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(tibble)
library(xts)
library(PerformanceAnalytics)
library(tidyr)
# ====== Q1. Import Data ======
library(tidyquant)
library(quantmod)
library(lubridate)
library(timetk)
library(purrr)
library(dplyr)
library(tibble)
library(xts)
# ETF tickers
tickers <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
# Download adjusted daily prices from Yahoo
etf_prices <- tq_get(tickers,
get = "stock.prices",
from = "2010-01-01",
to = Sys.Date())
# Show first few rows of raw prices
cat(" Downloaded ETF Prices:\n")
## Downloaded ETF Prices:
# Wide format by symbol
etf_prices_wide <- etf_prices %>%
select(symbol, date, adjusted) %>%
pivot_wider(names_from = symbol, values_from = adjusted) %>%
arrange(date)
# Show wide format
cat("\n Wide Format ETF Prices:\n")
##
## Wide Format ETF Prices:
print(head(etf_prices_wide))
## # A tibble: 6 × 9
## date SPY QQQ EEM IWM EFA TLT IYR GLD
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-04 85.8 40.5 31.1 51.9 36.4 58.3 27.4 110.
## 2 2010-01-05 86.0 40.5 31.3 51.7 36.4 58.6 27.5 110.
## 3 2010-01-06 86.1 40.2 31.4 51.7 36.6 57.8 27.5 112.
## 4 2010-01-07 86.4 40.3 31.2 52.1 36.4 57.9 27.7 111.
## 5 2010-01-08 86.7 40.6 31.4 52.4 36.7 57.9 27.5 111.
## 6 2010-01-11 86.8 40.4 31.4 52.1 37.0 57.6 27.6 113.
# Convert to xts
etf_xts <- tk_xts(etf_prices_wide, select = -date, date_var = date)
# Show xts object structure
print(head(etf_xts))
## SPY QQQ EEM IWM EFA TLT IYR
## 2010-01-04 85.76846 40.48581 31.08409 51.92010 36.38437 58.25040 27.40289
## 2010-01-05 85.99548 40.48581 31.30972 51.74158 36.41644 58.62663 27.46868
## 2010-01-06 86.05603 40.24161 31.37521 51.69288 36.57036 57.84178 27.45673
## 2010-01-07 86.41931 40.26777 31.19326 52.07430 36.42927 57.93910 27.70198
## 2010-01-08 86.70689 40.59919 31.44070 52.35832 36.71787 57.91313 27.51655
## 2010-01-11 86.82797 40.43348 31.37521 52.14735 37.01931 57.59534 27.64815
## GLD
## 2010-01-04 109.80
## 2010-01-05 109.70
## 2010-01-06 111.51
## 2010-01-07 110.82
## 2010-01-08 111.37
## 2010-01-11 112.85
# ====== Q2. Calculate Weekly and Monthly Returns ======
# Daily simple returns
daily_returns <- Return.calculate(etf_xts, method = "simple") %>% na.omit()
# Monthly returns: use map to calculate separately per ETF
etf_cols <- colnames(daily_returns)
monthly_returns_list <- map(etf_cols, function(etf) {
apply.monthly(daily_returns[, etf], function(x) prod(1 + x) - 1) %>%
`colnames<-`(etf)
})
monthly_returns <- reduce(monthly_returns_list, merge)
# Display the first 6 rows of monthly returns
head(monthly_returns)
## SPY QQQ EEM IWM EFA
## 2010-01-29 -0.05241329 -0.07819902 -0.103722724 -0.06048752 -0.07491636
## 2010-02-26 0.03119467 0.04603924 0.017763695 0.04475120 0.00266768
## 2010-03-31 0.06087961 0.07710884 0.081108890 0.08230691 0.06385433
## 2010-04-30 0.01547013 0.02242493 -0.001662095 0.05678478 -0.02804592
## 2010-05-28 -0.07945475 -0.07392348 -0.093935529 -0.07536629 -0.11192795
## 2010-06-30 -0.05174090 -0.05975694 -0.013986416 -0.07743427 -0.02061957
## TLT IYR GLD
## 2010-01-29 0.02783708 -0.05195369 -0.034972713
## 2010-02-26 -0.00342509 0.05457040 0.032748219
## 2010-03-31 -0.02057360 0.09748509 -0.004386396
## 2010-04-30 0.03321908 0.06388078 0.058834363
## 2010-05-28 0.05108366 -0.05683502 0.030513147
## 2010-06-30 0.05797770 -0.04670119 0.023553189
# ====== Q3. Convert Monthly Returns to Tibble ======
monthly_returns_tibble <- monthly_returns %>%
tk_tbl(preserve_index = TRUE, rename_index = "date")
# Show the first 6 rows of the tibble
head(monthly_returns_tibble)
## # A tibble: 6 × 9
## date SPY QQQ EEM IWM EFA TLT IYR GLD
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-29 -0.0524 -0.0782 -0.104 -0.0605 -0.0749 0.0278 -0.0520 -0.0350
## 2 2010-02-26 0.0312 0.0460 0.0178 0.0448 0.00267 -0.00343 0.0546 0.0327
## 3 2010-03-31 0.0609 0.0771 0.0811 0.0823 0.0639 -0.0206 0.0975 -0.00439
## 4 2010-04-30 0.0155 0.0224 -0.00166 0.0568 -0.0280 0.0332 0.0639 0.0588
## 5 2010-05-28 -0.0795 -0.0739 -0.0939 -0.0754 -0.112 0.0511 -0.0568 0.0305
## 6 2010-06-30 -0.0517 -0.0598 -0.0140 -0.0774 -0.0206 0.0580 -0.0467 0.0236
# ====== Q4. Download FF3 Data (Simulated version) ======
# You can replace this with actual Ken French data in production
# Recreate correct date vector based on monthly return dates
ff_dates <- monthly_returns_tibble$date
# Simulate FF3 data with same number of months
ff_data <- tibble(
date = ff_dates,
Mkt_RF = rnorm(n = length(ff_dates), mean = 0.007, sd = 0.04),
SMB = rnorm(n = length(ff_dates), mean = 0.002, sd = 0.025),
HML = rnorm(n = length(ff_dates), mean = 0.003, sd = 0.03),
RF = rnorm(n = length(ff_dates), mean = 0.001, sd = 0.002)
)
# Ensure date is in Date format
ff_data <- ff_data %>% mutate(date = as.Date(date))
# Show the first 6 rows of the FF data
head(ff_data)
## # A tibble: 6 × 5
## date Mkt_RF SMB HML RF
## <date> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-29 0.0157 0.0193 -0.0194 0.00502
## 2 2010-02-26 0.0455 -0.0290 -0.00434 -0.00141
## 3 2010-03-31 0.0425 -0.00295 0.0669 0.00102
## 4 2010-04-30 -0.0384 0.00304 -0.00440 -0.00190
## 5 2010-05-28 0.00225 0.00798 0.0186 0.00329
## 6 2010-06-30 -0.0285 -0.0288 -0.0157 0.00122
# ====== Q5. Merge Monthly Returns with FF3 ======
merged_data <- monthly_returns_tibble %>%
left_join(ff_data, by = "date")
# Show first 6 rows of the merged dataset
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-01-29 -0.0524 -0.0782 -0.104 -0.0605 -0.0749 0.0278 -0.0520 -0.0350
## 2 2010-02-26 0.0312 0.0460 0.0178 0.0448 0.00267 -0.00343 0.0546 0.0327
## 3 2010-03-31 0.0609 0.0771 0.0811 0.0823 0.0639 -0.0206 0.0975 -0.00439
## 4 2010-04-30 0.0155 0.0224 -0.00166 0.0568 -0.0280 0.0332 0.0639 0.0588
## 5 2010-05-28 -0.0795 -0.0739 -0.0939 -0.0754 -0.112 0.0511 -0.0568 0.0305
## 6 2010-06-30 -0.0517 -0.0598 -0.0140 -0.0774 -0.0206 0.0580 -0.0467 0.0236
## # ℹ 4 more variables: Mkt_RF <dbl>, SMB <dbl>, HML <dbl>, RF <dbl>
# ====== Q6. CAPM Covariance Matrix (2010/02 – 2015/01) ======
# Filter for specific 60-month period
capm_data <- merged_data %>%
filter(date >= as.Date("2010-02-01") & date <= as.Date("2015-01-31"))
# Calculate excess returns (ETF - RF)
excess_returns_capm <- capm_data %>%
select(all_of(c(etf_cols, "RF"))) %>%
mutate(across(all_of(etf_cols), ~ .x - RF)) %>%
select(-RF)
# Convert to matrix and compute covariance matrix
excess_matrix <- as.matrix(excess_returns_capm)
capm_cov_matrix <- cov(excess_matrix, use = "complete.obs")
rownames(capm_cov_matrix) <- etf_cols
colnames(capm_cov_matrix) <- etf_cols
# Output
cat("CAPM Covariance Matrix (2010/02–2015/01):\n")
## CAPM Covariance Matrix (2010/02–2015/01):
print(round(capm_cov_matrix, 6))
## SPY QQQ EEM IWM EFA TLT IYR
## SPY 0.001423 0.001496 0.001747 0.001816 0.001621 -0.000978 0.001226
## QQQ 0.001496 0.001847 0.001827 0.001886 0.001702 -0.000972 0.001263
## EEM 0.001747 0.001827 0.003366 0.002345 0.002484 -0.001189 0.001813
## IWM 0.001816 0.001886 0.002345 0.002730 0.001994 -0.001322 0.001639
## EFA 0.001621 0.001702 0.002484 0.001994 0.002429 -0.001112 0.001571
## TLT -0.000978 -0.000972 -0.001189 -0.001322 -0.001112 0.001629 -0.000368
## IYR 0.001226 0.001263 0.001813 0.001639 0.001571 -0.000368 0.002051
## GLD 0.000216 0.000416 0.000891 0.000557 0.000433 0.000205 0.000534
## GLD
## SPY 0.000216
## QQQ 0.000416
## EEM 0.000891
## IWM 0.000557
## EFA 0.000433
## TLT 0.000205
## IYR 0.000534
## GLD 0.002908
#Q7: FF3-Factor Covariance Matrix (2010/02 – 2015/01)
library(broom) # for tidy regression output if needed
# Step 1: Filter same 60-month period
ff3_data <- merged_data %>%
filter(date >= as.Date("2010-02-01") & date <= as.Date("2015-01-31"))
# Step 2: Create excess returns (R - RF)
excess_returns_ff3 <- ff3_data %>%
mutate(across(all_of(etf_cols), ~ .x - RF))
# Step 3: Run FF3 regression per ETF and extract residuals
residuals_list <- map(etf_cols, function(etf) {
formula <- as.formula(paste0(etf, " ~ Mkt_RF + SMB + HML"))
model <- lm(formula, data = excess_returns_ff3)
resid(model)
})
# Combine residuals into dataframe
residuals_df <- as.data.frame(residuals_list)
colnames(residuals_df) <- etf_cols
# Step 4: Compute residual covariance matrix
ff3_cov_matrix <- cov(residuals_df, use = "complete.obs")
# Step 5: Show results
cat("F3-Factor Residual Covariance Matrix (2010/02–2015/01):\n")
## F3-Factor Residual Covariance Matrix (2010/02–2015/01):
print(round(ff3_cov_matrix, 6))
## SPY QQQ EEM IWM EFA TLT IYR
## SPY 0.001334 0.001350 0.001577 0.001696 0.001502 -0.000881 0.001163
## QQQ 0.001350 0.001588 0.001554 0.001693 0.001502 -0.000834 0.001137
## EEM 0.001577 0.001554 0.003035 0.002114 0.002251 -0.000989 0.001702
## IWM 0.001696 0.001693 0.002114 0.002568 0.001831 -0.001186 0.001558
## EFA 0.001502 0.001502 0.002251 0.001831 0.002255 -0.000974 0.001486
## TLT -0.000881 -0.000834 -0.000989 -0.001186 -0.000974 0.001487 -0.000327
## IYR 0.001163 0.001137 0.001702 0.001558 0.001486 -0.000327 0.001981
## GLD 0.000181 0.000292 0.000866 0.000527 0.000400 0.000139 0.000425
## GLD
## SPY 0.000181
## QQQ 0.000292
## EEM 0.000866
## IWM 0.000527
## EFA 0.000400
## TLT 0.000139
## IYR 0.000425
## GLD 0.002589
# Question8:Compute global minimum variance portfolio weights
# Load quadprog package for optimization
library(quadprog)
# Define function to compute GMV weights
compute_gmv_weights <- function(cov_matrix) {
n <- ncol(cov_matrix)
Dmat <- cov_matrix
dvec <- rep(0, n) # No linear objective
Amat <- matrix(1, nrow = n) # Constraint: sum of weights = 1
bvec <- 1
result <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
weights <- result$solution
names(weights) <- colnames(cov_matrix)
return(weights)
}
# Compute GMV weights
gmv_capm_weights <- compute_gmv_weights(capm_cov_matrix)
gmv_ff3_weights <- compute_gmv_weights(ff3_cov_matrix)
# Display results
cat("GMV Portfolio Weights (CAPM):\n")
## GMV Portfolio Weights (CAPM):
print(round(gmv_capm_weights, 4))
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 0.9702 -0.2011 -0.0131 -0.0677 -0.0057 0.4727 -0.2095 0.0544
cat("\nMV Portfolio Weights (FF3):\n")
##
## MV Portfolio Weights (FF3):
print(round(gmv_ff3_weights, 4))
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 1.0149 -0.2204 -0.0326 -0.0879 -0.0192 0.4643 -0.1917 0.0726
library(quadprog)
library(ggplot2)
library(dplyr)
library(tibble)
library(purrr)
# Helper function: compute GMV return
get_gmv_return <- function(cov_matrix, next_month_returns) {
n <- ncol(cov_matrix)
Dmat <- cov_matrix
dvec <- rep(0, n)
Amat <- matrix(1, nrow = n)
bvec <- 1
result <- tryCatch(
solve.QP(Dmat, dvec, Amat, bvec, meq = 1),
error = function(e) return(rep(NA, n))
)
if (all(is.na(result))) return(NA)
weights <- result$solution
return(sum(weights * next_month_returns))
}
# Backtest function
gmv_backtest <- function(model = "CAPM") {
portfolio_returns <- c()
dates <- c()
n_obs <- nrow(merged_data)
for (i in (window_size + 1):(n_obs - 1)) {
data_window <- merged_data[(i - window_size):(i - 1), ]
next_month_returns <- merged_data[i, etf_cols] %>% as.numeric()
if (any(is.na(next_month_returns))) next
if (model == "FF3") {
residuals <- map(etf_cols, function(etf) {
df <- data_window %>%
select(all_of(c(etf_cols, "Mkt_RF", "SMB", "HML"))) %>%
mutate(y = data_window[[etf]])
model_ff <- tryCatch(
lm(y ~ Mkt_RF + SMB + HML, data = df),
error = function(e) return(NA)
)
if (inherits(model_ff, "lm")) {
return(resid(model_ff))
} else {
return(rep(NA, nrow(df)))
}
}) %>% bind_cols()
colnames(residuals) <- etf_cols
if (any(is.na(residuals))) next
cov_matrix <- cov(residuals, use = "complete.obs")
} else {
returns_window <- data_window %>%
mutate(across(all_of(etf_cols), ~ .x - RF)) %>%
select(all_of(etf_cols)) %>%
as.matrix()
cov_matrix <- cov(returns_window, use = "complete.obs")
}
if (any(is.na(cov_matrix)) || any(dim(cov_matrix) != c(length(etf_cols), length(etf_cols)))) next
ret <- get_gmv_return(cov_matrix, next_month_returns)
portfolio_returns <- c(portfolio_returns, ret)
dates <- c(dates, merged_data$date[i])
}
return(tibble(date = dates, return = portfolio_returns))
}
# Run backtest
window_size <- 60
etf_cols <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
gmv_capm <- gmv_backtest("CAPM") %>% mutate(model = "CAPM")
gmv_ff3 <- gmv_backtest("FF3") %>% mutate(model = "FF3")
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
# Combine and visualize
gmv_combined <- bind_rows(gmv_capm, gmv_ff3) %>%
group_by(model) %>%
arrange(date) %>%
mutate(cumulative = cumprod(1 + return))
ggplot(gmv_combined, aes(x = date, y = cumulative, color = model)) +
geom_line(size = 1.2) +
labs(title = "📈 GMV Portfolio Backtest: CAPM vs FF3",
x = "Date", y = "Cumulative Return") +
theme_minimal() +
scale_color_manual(values = c("CAPM" = "forestgreen", "FF3" = "steelblue"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Final performance stats
gmv_summary <- gmv_combined %>%
group_by(model) %>%
summarise(
avg_return = mean(return),
volatility = sd(return),
sharpe = mean(return) / sd(return),
final_value = last(cumulative)
)
print(gmv_summary)
## # A tibble: 2 × 5
## model avg_return volatility sharpe final_value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 CAPM 0.00388 0.0263 0.148 1.55
## 2 FF3 0.00423 0.0263 0.161 1.62