tickers <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
# Download adjusted prices from Yahoo Finance
getSymbols(tickers,
src = "yahoo",
from = "2010-01-01",
to = "2025-12-31",
auto.assign = TRUE)## [1] "SPY" "QQQ" "EEM" "IWM" "EFA" "TLT" "IYR" "GLD"
# Extract adjusted prices and merge
prices <- do.call(merge, lapply(tickers, function(t) Ad(get(t))))
colnames(prices) <- tickers
head(prices)## SPY QQQ EEM IWM EFA TLT IYR
## 2010-01-04 84.79637 40.29080 30.35152 51.36655 35.12844 56.13516 26.76812
## 2010-01-05 85.02084 40.29080 30.57181 51.18993 35.15940 56.49771 26.83237
## 2010-01-06 85.08072 40.04777 30.63576 51.14178 35.30801 55.74138 26.82070
## 2010-01-07 85.43986 40.07379 30.45810 51.51910 35.17178 55.83518 27.06027
## 2010-01-08 85.72416 40.40362 30.69972 51.80009 35.45043 55.81017 26.87912
## 2010-01-11 85.84389 40.23870 30.63576 51.59136 35.74147 55.50389 27.00768
## 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
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 2025-12-22 682.9648 618.4302 54.01 253.1297 95.70 86.39350 93.34799 408.23
## 2025-12-23 686.0863 621.3265 54.31 251.6324 96.29 86.53195 93.27818 413.64
## 2025-12-24 688.4997 623.1443 54.42 252.2613 96.41 87.05608 93.95628 411.93
## 2025-12-26 688.4299 623.1043 54.80 250.9736 96.57 86.76930 94.05600 416.74
## 2025-12-29 685.9766 620.0881 54.66 249.4363 96.28 87.09564 94.23550 398.60
## 2025-12-30 685.1389 618.6499 54.88 247.5896 96.44 86.88797 94.44491 398.89
# Convert daily adjusted prices to monthly prices (last trading day of each month)
prices_monthly <- to.monthly(prices, indexAt = "lastof", OHLC = FALSE)
# Discrete (simple) monthly returns: (P_t / P_{t-1}) - 1
returns_monthly <- Return.calculate(prices_monthly, method = "discrete")
returns_monthly <- returns_monthly[-1, ] # Remove first NA row
# Convert to data frame for merging
returns_df <- data.frame(
date = as.yearmon(index(returns_monthly)),
coredata(returns_monthly)
)
head(returns_df)## date SPY QQQ EEM IWM EFA
## 1 Feb 2010 0.03119460 0.04603869 0.017763847 0.04475136 0.002667503
## 2 Mar 2010 0.06087957 0.07710929 0.081109347 0.08230669 0.063854086
## 3 Apr 2010 0.01547016 0.02242545 -0.001662322 0.05678491 -0.028045558
## 4 May 2010 -0.07945456 -0.07392391 -0.093935888 -0.07536668 -0.111927991
## 5 Jun 2010 -0.05174101 -0.05975650 -0.013986079 -0.07743387 -0.020619541
## 6 Jul 2010 0.06830076 0.07258219 0.109324472 0.06730888 0.116104123
## TLT IYR GLD
## 1 -0.003423972 0.05457108 0.032748219
## 2 -0.020573551 0.09748380 -0.004386396
## 3 0.033219036 0.06388182 0.058834363
## 4 0.051082817 -0.05683569 0.030513147
## 5 0.057977934 -0.04670129 0.023553189
## 6 -0.009463760 0.09404773 -0.050871157
# Download FF3 factors using frenchdata package
ff3_raw <- download_french_data("Fama/French 3 Factors")
ff3 <- ff3_raw$subsets$data[[1]] # Monthly data
# Clean and format
ff3 <- ff3 %>%
rename(date_raw = date) %>%
mutate(
date = as.yearmon(as.character(date_raw), "%Y%m"),
`Mkt-RF` = as.numeric(`Mkt-RF`) / 100,
SMB = as.numeric(SMB) / 100,
HML = as.numeric(HML) / 100,
RF = as.numeric(RF) / 100
) %>%
filter(!is.na(date)) %>%
select(date, `Mkt-RF`, SMB, HML, RF)
head(ff3)## # A tibble: 6 × 5
## date `Mkt-RF` SMB HML RF
## <yearmon> <dbl> <dbl> <dbl> <dbl>
## 1 Jul 1926 0.0289 -0.0255 -0.0239 0.0022
## 2 Aug 1926 0.0264 -0.0114 0.0381 0.0025
## 3 Sep 1926 0.0038 -0.0136 0.0005 0.0023
## 4 Oct 1926 -0.0327 -0.0014 0.0082 0.0032
## 5 Nov 1926 0.0254 -0.0011 -0.0061 0.0031
## 6 Dec 1926 0.0262 -0.0007 0.0006 0.0028
## # A tibble: 6 × 5
## date `Mkt-RF` SMB HML RF
## <yearmon> <dbl> <dbl> <dbl> <dbl>
## 1 Sep 2025 0.0339 -0.0185 -0.0105 0.0033
## 2 Oct 2025 0.0196 -0.0055 -0.031 0.0037
## 3 Nov 2025 -0.0013 0.0038 0.0376 0.003
## 4 Dec 2025 -0.0036 -0.0106 0.0242 0.0034
## 5 Jan 2026 0.0102 0.022 0.0372 0.003
## 6 Feb 2026 -0.0117 0.0014 0.0283 0.0028
merged_df <- inner_join(returns_df, ff3, by = "date")
# Compute excess returns for each ETF (R_i - RF)
for (t in tickers) {
merged_df[[paste0(t, "_excess")]] <- merged_df[[t]] - merged_df$RF
}
cat("Merged data dimensions:", dim(merged_df), "\n")## Merged data dimensions: 191 21
cat("Date range:", as.character(min(merged_df$date)),
"to", as.character(max(merged_df$date)), "\n")## Date range: Feb 2010 to Dec 2025
## date SPY QQQ EEM IWM EFA
## 1 Feb 2010 0.03119460 0.04603869 0.017763847 0.04475136 0.002667503
## 2 Mar 2010 0.06087957 0.07710929 0.081109347 0.08230669 0.063854086
## 3 Apr 2010 0.01547016 0.02242545 -0.001662322 0.05678491 -0.028045558
## 4 May 2010 -0.07945456 -0.07392391 -0.093935888 -0.07536668 -0.111927991
## 5 Jun 2010 -0.05174101 -0.05975650 -0.013986079 -0.07743387 -0.020619541
## 6 Jul 2010 0.06830076 0.07258219 0.109324472 0.06730888 0.116104123
## TLT IYR GLD Mkt-RF SMB HML RF
## 1 -0.003423972 0.05457108 0.032748219 0.0339 0.0118 0.0318 0e+00
## 2 -0.020573551 0.09748380 -0.004386396 0.0630 0.0146 0.0219 1e-04
## 3 0.033219036 0.06388182 0.058834363 0.0199 0.0484 0.0296 1e-04
## 4 0.051082817 -0.05683569 0.030513147 -0.0790 0.0013 -0.0248 1e-04
## 5 0.057977934 -0.04670129 0.023553189 -0.0556 -0.0179 -0.0473 1e-04
## 6 -0.009463760 0.09404773 -0.050871157 0.0692 0.0022 -0.0050 1e-04
## SPY_excess QQQ_excess EEM_excess IWM_excess EFA_excess TLT_excess
## 1 0.03119460 0.04603869 0.017763847 0.04475136 0.002667503 -0.003423972
## 2 0.06077957 0.07700929 0.081009347 0.08220669 0.063754086 -0.020673551
## 3 0.01537016 0.02232545 -0.001762322 0.05668491 -0.028145558 0.033119036
## 4 -0.07955456 -0.07402391 -0.094035888 -0.07546668 -0.112027991 0.050982817
## 5 -0.05184101 -0.05985650 -0.014086079 -0.07753387 -0.020719541 0.057877934
## 6 0.06820076 0.07248219 0.109224472 0.06720888 0.116004123 -0.009563760
## IYR_excess GLD_excess
## 1 0.05457108 0.032748219
## 2 0.09738380 -0.004486396
## 3 0.06378182 0.058734363
## 4 -0.05693569 0.030413147
## 5 -0.04680129 0.023453189
## 6 0.09394773 -0.050971157
# Solve for MVP weights using quadprog
# Minimise: w' Sigma w subject to: sum(w) = 1, w >= 0
compute_mvp_weights <- function(cov_matrix) {
n <- ncol(cov_matrix)
Dmat <- 2 * cov_matrix
dvec <- rep(0, n)
# Constraints: (1) sum(w) = 1 (2) w_i >= 0
Amat <- cbind(rep(1, n), diag(n))
bvec <- c(1, rep(0, n))
sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
weights <- sol$solution
names(weights) <- colnames(cov_matrix)
return(weights)
}Under CAPM, the return of asset \(i\) is:
\[R_i - R_f = \alpha_i + \beta_i (R_m - R_f) + \epsilon_i\]
The covariance matrix is estimated as:
\[\Sigma = \hat{\beta}\hat{\beta}' \sigma^2_{mkt} + D\]
where \(D\) is the diagonal matrix of residual variances.
# Filter: 60-month window 2020/03 – 2025/02
window_data <- merged_df %>%
filter(date >= as.yearmon("2020-03") & date <= as.yearmon("2025-02"))
cat("Window observations:", nrow(window_data), "\n")## Window observations: 60
# Fit CAPM for each ETF: excess_return ~ Mkt-RF
capm_betas <- numeric(length(tickers))
capm_alphas <- numeric(length(tickers))
capm_resid_var <- numeric(length(tickers))
for (i in seq_along(tickers)) {
t <- tickers[i]
exc <- window_data[[paste0(t, "_excess")]]
mkt <- window_data$`Mkt-RF`
fit <- lm(exc ~ mkt)
capm_betas[i] <- coef(fit)[2]
capm_alphas[i] <- coef(fit)[1]
capm_resid_var[i] <- var(residuals(fit))
}
names(capm_betas) <- tickers
names(capm_resid_var) <- tickers
# Market variance
sigma2_mkt <- var(window_data$`Mkt-RF`)
# CAPM covariance matrix: Sigma = beta %*% t(beta) * sigma2_mkt + D
capm_cov <- outer(capm_betas, capm_betas) * sigma2_mkt +
diag(capm_resid_var)
colnames(capm_cov) <- rownames(capm_cov) <- tickers
cat("\nCAPM Betas:\n")##
## CAPM Betas:
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 0.9552 1.0634 0.6963 1.1858 0.8243 0.3310 1.0036 0.1746
##
## CAPM Covariance Matrix:
## SPY QQQ EEM IWM EFA TLT IYR GLD
## SPY 0.002623 0.002884 0.001888 0.003216 0.002236 0.000898 0.002722 0.000473
## QQQ 0.002884 0.003791 0.002102 0.003580 0.002489 0.000999 0.003030 0.000527
## EEM 0.001888 0.002102 0.002712 0.002344 0.001629 0.000654 0.001984 0.000345
## IWM 0.003216 0.003580 0.002344 0.004983 0.002775 0.001114 0.003378 0.000588
## EFA 0.002236 0.002489 0.001629 0.002775 0.002618 0.000775 0.002349 0.000409
## TLT 0.000898 0.000999 0.000654 0.001114 0.000775 0.001934 0.000943 0.000164
## IYR 0.002722 0.003030 0.001984 0.003378 0.002349 0.000943 0.003888 0.000497
## GLD 0.000473 0.000527 0.000345 0.000588 0.000409 0.000164 0.000497 0.001733
Under FF3, the return of asset \(i\) is:
\[R_i - R_f = \alpha_i + \beta_{i,mkt}(R_m-R_f) + \beta_{i,SMB} \cdot SMB + \beta_{i,HML} \cdot HML + \epsilon_i\]
\[\Sigma = B \Sigma_F B' + D\]
where \(B\) is the \(n \times 3\) matrix of factor loadings, \(\Sigma_F\) is the \(3 \times 3\) factor covariance matrix, and \(D\) is the diagonal residual variance matrix.
ff3_betas <- matrix(NA, nrow = length(tickers), ncol = 3,
dimnames = list(tickers, c("Mkt-RF","SMB","HML")))
ff3_resid_var <- numeric(length(tickers))
names(ff3_resid_var) <- tickers
for (i in seq_along(tickers)) {
t <- tickers[i]
exc <- window_data[[paste0(t, "_excess")]]
fit <- lm(exc ~ `Mkt-RF` + SMB + HML, data = window_data)
ff3_betas[t, ] <- coef(fit)[2:4]
ff3_resid_var[t] <- var(residuals(fit))
}
# Factor covariance matrix (3x3)
factor_cov <- cov(window_data[, c("Mkt-RF","SMB","HML")])
# FF3 covariance matrix: Sigma = B %*% Sigma_F %*% t(B) + D
ff3_cov <- ff3_betas %*% factor_cov %*% t(ff3_betas) + diag(ff3_resid_var)
colnames(ff3_cov) <- rownames(ff3_cov) <- tickers
cat("FF3 Factor Loadings (Betas):\n")## FF3 Factor Loadings (Betas):
## Mkt-RF SMB HML
## SPY 0.9853 -0.1487 0.0194
## QQQ 1.0813 -0.0890 -0.3994
## EEM 0.6794 0.0834 0.1476
## IWM 1.0058 0.8895 0.2660
## EFA 0.8477 -0.1152 0.2169
## TLT 0.3443 -0.0658 -0.2622
## IYR 0.9953 0.0409 0.2032
## GLD 0.2420 -0.3330 -0.0197
##
## FF3 Covariance Matrix:
## SPY QQQ EEM IWM EFA TLT IYR GLD
## SPY 0.002623 0.002885 0.001881 0.003102 0.002257 0.000899 0.002722 0.000518
## QQQ 0.002885 0.003791 0.001964 0.003243 0.002314 0.001235 0.002848 0.000584
## EEM 0.001881 0.001964 0.002712 0.002511 0.001690 0.000563 0.002054 0.000308
## IWM 0.003102 0.003243 0.002511 0.004983 0.002819 0.000887 0.003546 0.000292
## EFA 0.002257 0.002314 0.001690 0.002819 0.002618 0.000661 0.002438 0.000428
## TLT 0.000899 0.001235 0.000563 0.000887 0.000661 0.001934 0.000824 0.000204
## IYR 0.002722 0.002848 0.002054 0.003546 0.002438 0.000824 0.003888 0.000470
## GLD 0.000518 0.000584 0.000308 0.000292 0.000428 0.000204 0.000470 0.001733
## FF3 MVP Weights:
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 0.0000 0.0000 0.1565 0.0000 0.0821 0.3391 0.0000 0.4223
# Bar plot
barplot(ff3_weights,
main = "FF3 MVP Weights (Window: 2020/03–2025/02)",
ylab = "Weight",
col = "darkorange",
las = 2)
abline(h = 0, lty = 2)weights_compare <- data.frame(
ETF = tickers,
CAPM = round(capm_weights, 4),
FF3 = round(ff3_weights, 4)
)
print(weights_compare)## ETF CAPM FF3
## SPY SPY 0.0000 0.0000
## QQQ QQQ 0.0000 0.0000
## EEM EEM 0.1401 0.1565
## IWM IWM 0.0000 0.0000
## EFA EFA 0.0838 0.0821
## TLT TLT 0.3425 0.3391
## IYR IYR 0.0000 0.0000
## GLD GLD 0.4336 0.4223
# Side-by-side bar chart
barplot(
t(as.matrix(weights_compare[, c("CAPM","FF3")])),
beside = TRUE,
names.arg = tickers,
col = c("steelblue","darkorange"),
legend = c("CAPM","FF3"),
main = "MVP Weights: CAPM vs FF3",
ylab = "Weight",
las = 2
)# Get March 2025 actual returns
march2025 <- merged_df %>%
filter(date == as.yearmon("2025-03"))
if (nrow(march2025) == 0) {
# If not in merged_df yet, fetch directly from prices_monthly
march2025_returns <- as.numeric(
coredata(returns_monthly[as.yearmon(index(returns_monthly)) == as.yearmon("2025-03"), ])
)
names(march2025_returns) <- tickers
} else {
march2025_returns <- as.numeric(march2025[1, tickers])
names(march2025_returns) <- tickers
}
cat("March 2025 ETF Returns:\n")## March 2025 ETF Returns:
## SPY QQQ EEM IWM EFA TLT IYR GLD
## -0.0557 -0.0759 0.0113 -0.0685 0.0018 -0.0120 -0.0234 0.0945
# Realized portfolio returns = w' * r
ret_capm_mar <- sum(capm_weights * march2025_returns)
ret_ff3_mar <- sum(ff3_weights * march2025_returns)
cat("\n--- Realized MVP Returns: March 2025 ---\n")##
## --- Realized MVP Returns: March 2025 ---
## CAPM MVP Return: 0.0386 (3.86%)
## FF3 MVP Return: 0.0377 (3.77%)
For April 2025, the estimation window shifts to 2020/04 – 2025/03 (rolling forward by one month).
# New 60-month window: 2020/04 – 2025/03
window_apr <- merged_df %>%
filter(date >= as.yearmon("2020-04") & date <= as.yearmon("2025-03"))
cat("April window observations:", nrow(window_apr), "\n")## April window observations: 60
# --- CAPM Covariance for April window ---
capm_betas_apr <- numeric(length(tickers))
capm_resid_apr <- numeric(length(tickers))
for (i in seq_along(tickers)) {
t <- tickers[i]
exc <- window_apr[[paste0(t, "_excess")]]
mkt <- window_apr$`Mkt-RF`
fit <- lm(exc ~ mkt)
capm_betas_apr[i] <- coef(fit)[2]
capm_resid_apr[i] <- var(residuals(fit))
}
names(capm_betas_apr) <- tickers
names(capm_resid_apr) <- tickers
sigma2_mkt_apr <- var(window_apr$`Mkt-RF`)
capm_cov_apr <- outer(capm_betas_apr, capm_betas_apr) * sigma2_mkt_apr +
diag(capm_resid_apr)
colnames(capm_cov_apr) <- rownames(capm_cov_apr) <- tickers
capm_weights_apr <- compute_mvp_weights(capm_cov_apr)
# --- FF3 Covariance for April window ---
ff3_betas_apr <- matrix(NA, nrow = length(tickers), ncol = 3,
dimnames = list(tickers, c("Mkt-RF","SMB","HML")))
ff3_resid_apr <- numeric(length(tickers))
names(ff3_resid_apr) <- tickers
for (i in seq_along(tickers)) {
t <- tickers[i]
exc <- window_apr[[paste0(t, "_excess")]]
fit <- lm(exc ~ `Mkt-RF` + SMB + HML, data = window_apr)
ff3_betas_apr[t, ] <- coef(fit)[2:4]
ff3_resid_apr[t] <- var(residuals(fit))
}
factor_cov_apr <- cov(window_apr[, c("Mkt-RF","SMB","HML")])
ff3_cov_apr <- ff3_betas_apr %*% factor_cov_apr %*% t(ff3_betas_apr) +
diag(ff3_resid_apr)
colnames(ff3_cov_apr) <- rownames(ff3_cov_apr) <- tickers
ff3_weights_apr <- compute_mvp_weights(ff3_cov_apr)
cat("\nCAPM MVP Weights (April 2025 window):\n")##
## CAPM MVP Weights (April 2025 window):
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 0.0000 0.0000 0.1847 0.0000 0.1140 0.3046 0.0000 0.3967
##
## FF3 MVP Weights (April 2025 window):
## SPY QQQ EEM IWM EFA TLT IYR GLD
## 0.0000 0.0000 0.1949 0.0000 0.1051 0.3064 0.0000 0.3936
# --- April 2025 Realized Returns ---
april2025 <- merged_df %>%
filter(date == as.yearmon("2025-04"))
if (nrow(april2025) == 0) {
april2025_returns <- as.numeric(
coredata(returns_monthly[as.yearmon(index(returns_monthly)) == as.yearmon("2025-04"), ])
)
names(april2025_returns) <- tickers
} else {
april2025_returns <- as.numeric(april2025[1, tickers])
names(april2025_returns) <- tickers
}
cat("\nApril 2025 ETF Returns:\n")##
## April 2025 ETF Returns:
## SPY QQQ EEM IWM EFA TLT IYR GLD
## -0.0087 0.0140 0.0014 -0.0232 0.0370 -0.0136 -0.0215 0.0542
ret_capm_apr <- sum(capm_weights_apr * april2025_returns)
ret_ff3_apr <- sum(ff3_weights_apr * april2025_returns)
cat("\n--- Realized MVP Returns: April 2025 ---\n")##
## --- Realized MVP Returns: April 2025 ---
## CAPM MVP Return: 0.0218 (2.18%)
## FF3 MVP Return: 0.0213 (2.13%)
summary_tbl <- data.frame(
Model = c("CAPM", "FF3"),
March_2025_Return = c(
paste0(round(ret_capm_mar * 100, 2), "%"),
paste0(round(ret_ff3_mar * 100, 2), "%")
),
April_2025_Return = c(
paste0(round(ret_capm_apr * 100, 2), "%"),
paste0(round(ret_ff3_apr * 100, 2), "%")
)
)
knitr::kable(summary_tbl,
caption = "Realized MVP Portfolio Returns",
align = "lcc")| Model | March_2025_Return | April_2025_Return |
|---|---|---|
| CAPM | 3.86% | 2.18% |
| FF3 | 3.77% | 2.13% |