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.76845 40.48582 31.08410 51.92011 36.38437 58.25041 27.40289
## 2010-01-05 85.99549 40.48582 31.30970 51.74158 36.41643 58.62660 27.46869
## 2010-01-06 86.05602 40.24161 31.37521 51.69289 36.57036 57.84180 27.45673
## 2010-01-07 86.41930 40.26778 31.19325 52.07430 36.42926 57.93911 27.70199
## 2010-01-08 86.70685 40.59920 31.44071 52.35833 36.71788 57.91312 27.51654
## 2010-01-11 86.82796 40.43348 31.37521 52.14734 37.01931 57.59533 27.64816
## 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.05241330 -0.07819965 -0.103723073 -0.06048780 -0.074916242
## 2010-02-26 0.03119467 0.04603905 0.017764178 0.04475137 0.002667793
## 2010-03-31 0.06087998 0.07710927 0.081108397 0.08230654 0.063853738
## 2010-04-30 0.01546961 0.02242539 -0.001661659 0.05678522 -0.028045505
## 2010-05-28 -0.07945460 -0.07392399 -0.093935755 -0.07536642 -0.111928143
## 2010-06-30 -0.05174100 -0.05975676 -0.013986416 -0.07743413 -0.020619632
## TLT IYR GLD
## 2010-01-29 0.02783616 -0.05195355 -0.034972713
## 2010-02-26 -0.00342369 0.05457076 0.032748219
## 2010-03-31 -0.02057352 0.09748406 -0.004386396
## 2010-04-30 0.03321736 0.06388169 0.058834363
## 2010-05-28 0.05108440 -0.05683541 0.030513147
## 2010-06-30 0.05797801 -0.04670138 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.00342 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.0125 -0.000197 0.0139 -0.0000503
## 2 2010-02-26 0.0123 -0.0108 -0.00863 0.000862
## 3 2010-03-31 0.0446 0.0391 0.0170 0.00188
## 4 2010-04-30 -0.00770 0.0143 -0.00317 -0.00169
## 5 2010-05-28 -0.0206 -0.0188 0.0149 0.000977
## 6 2010-06-30 -0.0354 -0.0457 -0.0558 0.00131
# ====== 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.00342 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.001389 0.001454 0.001719 0.001776 0.001593 -0.000982 0.001210
## QQQ 0.001454 0.001797 0.001791 0.001839 0.001666 -0.000984 0.001239
## EEM 0.001719 0.001791 0.003344 0.002311 0.002462 -0.001187 0.001803
## IWM 0.001776 0.001839 0.002311 0.002684 0.001960 -0.001332 0.001618
## EFA 0.001593 0.001666 0.002462 0.001960 0.002407 -0.001111 0.001562
## TLT -0.000982 -0.000984 -0.001187 -0.001332 -0.001111 0.001654 -0.000354
## IYR 0.001210 0.001239 0.001803 0.001618 0.001562 -0.000354 0.002053
## GLD 0.000215 0.000407 0.000896 0.000550 0.000437 0.000233 0.000550
## GLD
## SPY 0.000215
## QQQ 0.000407
## EEM 0.000896
## IWM 0.000550
## EFA 0.000437
## TLT 0.000233
## IYR 0.000550
## GLD 0.002939
#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.001261 0.001328 0.001682 0.001590 0.001527 -0.000878 0.001107
## QQQ 0.001328 0.001663 0.001770 0.001647 0.001610 -0.000871 0.001155
## EEM 0.001682 0.001770 0.003296 0.002269 0.002432 -0.001182 0.001722
## IWM 0.001590 0.001647 0.002269 0.002407 0.001874 -0.001174 0.001480
## EFA 0.001527 0.001610 0.002432 0.001874 0.002360 -0.001064 0.001499
## TLT -0.000878 -0.000871 -0.001182 -0.001174 -0.001064 0.001554 -0.000305
## IYR 0.001107 0.001155 0.001722 0.001480 0.001499 -0.000305 0.001895
## GLD 0.000207 0.000399 0.000936 0.000548 0.000411 0.000270 0.000629
## GLD
## SPY 0.000207
## QQQ 0.000399
## EEM 0.000936
## IWM 0.000548
## EFA 0.000411
## TLT 0.000270
## IYR 0.000629
## GLD 0.002753
# 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.9468 -0.1690 -0.0104 -0.0527 -0.0117 0.4763 -0.2239 0.0447
cat("\nMV Portfolio Weights (FF3):\n")
##
## MV Portfolio Weights (FF3):
print(round(gmv_ff3_weights, 4))
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 0.9861 -0.1847 -0.0321 -0.0536 -0.0131 0.4598 -0.2179 0.0555
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.00409 0.0262 0.156 1.60
## 2 FF3 0.00419 0.0260 0.161 1.62