We use quantmod to download adjusted daily prices for
the 8 ETFs from Yahoo Finance.
library(quantmod)
library(tidyverse)
library(lubridate)
library(PerformanceAnalytics)
library(quadprog)
library(knitr)
library(kableExtra)
tickers <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
# Download adjusted daily prices
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 close prices and merge into one xts object
prices_daily <- do.call(merge, lapply(tickers, function(tk) Ad(get(tk))))
colnames(prices_daily) <- tickers
head(prices_daily)## SPY QQQ EEM IWM EFA TLT IYR
## 2010-01-04 84.79637 40.29079 30.35151 51.36655 35.12844 56.13515 26.76810
## 2010-01-05 85.02084 40.29079 30.57180 51.18993 35.15940 56.49769 26.83238
## 2010-01-06 85.08067 40.04774 30.63576 51.14176 35.30801 55.74141 26.82071
## 2010-01-07 85.43984 40.07379 30.45811 51.51911 35.17179 55.83516 27.06028
## 2010-01-08 85.72419 40.40363 30.69972 51.80010 35.45044 55.81015 26.87913
## 2010-01-11 85.84389 40.23870 30.63576 51.59137 35.74148 55.50392 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.76929 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.89Convert daily adjusted prices to monthly prices, then compute discrete (simple) returns: \[R_t = \frac{P_t - P_{t-1}}{P_{t-1}}\]
# Convert to monthly (last trading day of each month)
prices_monthly <- to.monthly(prices_daily, indexAt = "lastof", OHLC = FALSE)
# Discrete (simple) monthly returns
returns_monthly <- na.omit(Return.calculate(prices_monthly, method = "discrete"))
head(returns_monthly)## SPY QQQ EEM IWM EFA
## 2010-02-28 0.03119509 0.04603870 0.017763843 0.04475129 0.002667738
## 2010-03-31 0.06087956 0.07710940 0.081108848 0.08230693 0.063854079
## 2010-04-30 0.01547007 0.02242564 -0.001661812 0.05678441 -0.028045885
## 2010-05-31 -0.07945448 -0.07392408 -0.093936056 -0.07536598 -0.111927677
## 2010-06-30 -0.05174110 -0.05975717 -0.013986641 -0.07743434 -0.020619793
## 2010-07-31 0.06830067 0.07258326 0.109325026 0.06730912 0.116104269
## TLT IYR GLD
## 2010-02-28 -0.003424896 0.05456963 0.032748219
## 2010-03-31 -0.020573687 0.09748511 -0.004386396
## 2010-04-30 0.033217760 0.06388116 0.058834363
## 2010-05-31 0.051084471 -0.05683504 0.030513147
## 2010-06-30 0.057978433 -0.04670160 0.023553189
## 2010-07-31 -0.009463991 0.09404820 -0.050871157## SPY QQQ EEM IWM EFA
## 2025-07-31 0.023031543 0.024236854 0.006633454 0.0166829141 -0.020919628
## 2025-08-31 0.020519538 0.009539755 0.026770994 0.0719267543 0.045246876
## 2025-09-30 0.035620413 0.053762154 0.070998795 0.0317911543 0.020660248
## 2025-10-31 0.023837414 0.047803817 0.035580494 0.0176475324 0.011995297
## 2025-11-30 0.001949944 -0.015610297 -0.017721479 0.0102343663 0.007408238
## 2025-12-31 0.008267742 0.001579351 0.024767431 0.0004491716 0.031490000
## TLT IYR GLD
## 2025-07-31 -0.0113966251 0.001160733 -0.006134551
## 2025-08-31 0.0001270414 0.029089345 0.049874625
## 2025-09-30 0.0359097587 0.000554634 0.117584158
## 2025-10-31 0.0138109048 -0.024927912 0.035586671
## 2025-11-30 0.0027232220 0.023663708 0.053678176
## 2025-12-31 -0.0187683309 -0.013658495 0.028385092Download from Ken French’s data library and rescale from percentage to decimal.
library(frenchdata) # install.packages("frenchdata") if needed
# Download FF3 monthly factors
ff3_raw <- download_french_data("Fama/French 3 Factors")
ff3_monthly <- ff3_raw$subsets$data[[1]] # Monthly data subset
# Rename and convert to decimal
ff3_monthly <- ff3_monthly %>%
rename(
Date = date,
Mkt_RF = `Mkt-RF`,
SMB = SMB,
HML = HML,
RF = RF
) %>%
mutate(
Date = ymd(paste0(Date, "01")), # "YYYYMM" -> Date
Mkt_RF = Mkt_RF / 100,
SMB = SMB / 100,
HML = HML / 100,
RF = RF / 100
) %>%
filter(Date >= as.Date("2010-01-01"),
Date <= as.Date("2025-12-31"))
head(ff3_monthly)## # A tibble: 6 × 5
## Date Mkt_RF SMB HML RF
## <date> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-01 -0.0335 0.0043 0.0033 0
## 2 2010-02-01 0.0339 0.0118 0.0318 0
## 3 2010-03-01 0.063 0.0146 0.0219 0.0001
## 4 2010-04-01 0.0199 0.0484 0.0296 0.0001
## 5 2010-05-01 -0.079 0.0013 -0.0248 0.0001
## 6 2010-06-01 -0.0556 -0.0179 -0.0473 0.0001# Convert xts returns to data frame with Date column
returns_df <- as.data.frame(returns_monthly) %>%
rownames_to_column("Date") %>%
mutate(Date = as.Date(as.yearmon(Date))) # match end-of-month format
# Merge on month-year
merged_df <- returns_df %>%
left_join(ff3_monthly, by = "Date") %>%
drop_na()
cat("Merged dataset dimensions:", nrow(merged_df), "rows x", ncol(merged_df), "cols\n")## Merged dataset dimensions: 191 rows x 13 cols## Date SPY QQQ EEM IWM EFA
## 1 2010-02-01 0.03119509 0.04603870 0.017763843 0.04475129 0.002667738
## 2 2010-03-01 0.06087956 0.07710940 0.081108848 0.08230693 0.063854079
## 3 2010-04-01 0.01547007 0.02242564 -0.001661812 0.05678441 -0.028045885
## 4 2010-05-01 -0.07945448 -0.07392408 -0.093936056 -0.07536598 -0.111927677
## 5 2010-06-01 -0.05174110 -0.05975717 -0.013986641 -0.07743434 -0.020619793
## 6 2010-07-01 0.06830067 0.07258326 0.109325026 0.06730912 0.116104269
## TLT IYR GLD Mkt_RF SMB HML RF
## 1 -0.003424896 0.05456963 0.032748219 0.0339 0.0118 0.0318 0e+00
## 2 -0.020573687 0.09748511 -0.004386396 0.0630 0.0146 0.0219 1e-04
## 3 0.033217760 0.06388116 0.058834363 0.0199 0.0484 0.0296 1e-04
## 4 0.051084471 -0.05683504 0.030513147 -0.0790 0.0013 -0.0248 1e-04
## 5 0.057978433 -0.04670160 0.023553189 -0.0556 -0.0179 -0.0473 1e-04
## 6 -0.009463991 0.09404820 -0.050871157 0.0692 0.0022 -0.0050 1e-04Method: Estimate betas from CAPM, build the factor-model covariance matrix, then solve for the Minimum Variance Portfolio (MVP) weights.
\[\Sigma_{CAPM} = \beta \sigma_m^2 \beta^\top + D\]
where \(D\) is the diagonal matrix of residual variances.
# --- Helper: solve MVP weights given covariance matrix ---
solve_mvp <- function(sigma) {
n <- ncol(sigma)
Dmat <- 2 * sigma
dvec <- rep(0, n)
Amat <- cbind(rep(1, n), # equality: weights sum to 1
diag(n)) # inequality: weights >= 0 (long-only)
bvec <- c(1, rep(0, n))
meq <- 1L
sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = meq)
w <- sol$solution
names(w) <- colnames(sigma)
w
}
# --- Window: 2020/03 – 2025/02 (60 months) ---
window_start <- as.Date("2020-03-01")
window_end <- as.Date("2025-02-28")
window_data <- merged_df %>%
filter(Date >= window_start, Date <= window_end)
cat("Window observations:", nrow(window_data), "\n")## Window observations: 60# Excess returns for each ETF
etf_excess <- window_data %>%
select(all_of(tickers)) %>%
sweep(1, window_data$RF, "-") # R_i - RF
mkt_excess <- window_data$Mkt_RF # already excess
# --- OLS betas (CAPM: one factor) ---
betas_capm <- sapply(tickers, function(tk) {
fit <- lm(etf_excess[[tk]] ~ mkt_excess)
coef(fit)[2]
})
# Residual variances
resid_var_capm <- sapply(tickers, function(tk) {
fit <- lm(etf_excess[[tk]] ~ mkt_excess)
var(residuals(fit))
})
# Market variance
var_mkt <- var(mkt_excess)
# CAPM covariance matrix
sigma_capm <- outer(betas_capm, betas_capm) * var_mkt +
diag(resid_var_capm)
colnames(sigma_capm) <- rownames(sigma_capm) <- tickers
# MVP weights
w_mvp_capm <- solve_mvp(sigma_capm)
kable(data.frame(Ticker = tickers,
Beta = round(betas_capm, 4),
Weight = round(w_mvp_capm, 4)),
caption = "CAPM – Betas and MVP Weights") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Ticker | Beta | Weight | |
|---|---|---|---|
| SPY.mkt_excess | SPY | 0.9552 | 0.0000 |
| QQQ.mkt_excess | QQQ | 1.0634 | 0.0000 |
| EEM.mkt_excess | EEM | 0.6963 | 0.1401 |
| IWM.mkt_excess | IWM | 1.1858 | 0.0000 |
| EFA.mkt_excess | EFA | 0.8243 | 0.0838 |
| TLT.mkt_excess | TLT | 0.3310 | 0.3425 |
| IYR.mkt_excess | IYR | 1.0036 | 0.0000 |
| GLD.mkt_excess | GLD | 0.1746 | 0.4336 |
\[\Sigma_{FF3} = B \Sigma_F B^\top + D\]
where \(B\) is the \((8 \times 3)\) factor-loading matrix and \(\Sigma_F\) is the \((3 \times 3)\) factor covariance matrix.
factors_ff3 <- window_data[, c("Mkt_RF", "SMB", "HML")]
# --- OLS factor loadings for each ETF ---
B_ff3 <- sapply(tickers, function(tk) {
fit <- lm(etf_excess[[tk]] ~ Mkt_RF + SMB + HML, data = window_data)
coef(fit)[2:4] # Mkt_RF, SMB, HML loadings
})
B_ff3 <- t(B_ff3) # 8 x 3 matrix
colnames(B_ff3) <- c("Mkt_RF", "SMB", "HML")
# Factor covariance matrix (3x3)
Sigma_F <- cov(factors_ff3)
# Residual variances
resid_var_ff3 <- sapply(tickers, function(tk) {
fit <- lm(etf_excess[[tk]] ~ Mkt_RF + SMB + HML, data = window_data)
var(residuals(fit))
})
# FF3 covariance matrix
sigma_ff3 <- B_ff3 %*% Sigma_F %*% t(B_ff3) + diag(resid_var_ff3)
colnames(sigma_ff3) <- rownames(sigma_ff3) <- tickers
# MVP weights
w_mvp_ff3 <- solve_mvp(sigma_ff3)
kable(data.frame(Ticker = tickers,
b_MktRF = round(B_ff3[, "Mkt_RF"], 4),
b_SMB = round(B_ff3[, "SMB"], 4),
b_HML = round(B_ff3[, "HML"], 4),
Weight_FF3 = round(w_mvp_ff3, 4)),
caption = "FF3 – Factor Loadings and MVP Weights") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Ticker | b_MktRF | b_SMB | b_HML | Weight_FF3 | |
|---|---|---|---|---|---|
| SPY | SPY | 0.9853 | -0.1487 | 0.0194 | 0.0000 |
| QQQ | QQQ | 1.0813 | -0.0890 | -0.3994 | 0.0000 |
| EEM | EEM | 0.6794 | 0.0834 | 0.1476 | 0.1565 |
| IWM | IWM | 1.0058 | 0.8895 | 0.2660 | 0.0000 |
| EFA | EFA | 0.8477 | -0.1152 | 0.2169 | 0.0821 |
| TLT | TLT | 0.3443 | -0.0658 | -0.2622 | 0.3391 |
| IYR | IYR | 0.9953 | 0.0409 | 0.2032 | 0.0000 |
| GLD | GLD | 0.2420 | -0.3330 | -0.0197 | 0.4223 |
Apply the MVP weights (estimated on 2020/03–2025/02) to the actual returns in March 2025.
# March 2025 returns
mar2025 <- merged_df %>%
filter(year(Date) == 2025, month(Date) == 3) %>%
select(all_of(tickers))
if (nrow(mar2025) == 0) {
# Fallback: compute from daily prices if monthly merge doesn't include Mar 2025
prices_mar <- prices_monthly["2025-03"]
prices_feb <- prices_monthly["2025-02"]
mar2025_ret <- as.numeric((prices_mar - prices_feb) / prices_feb)
names(mar2025_ret) <- tickers
mar2025 <- as.data.frame(t(mar2025_ret))
}
ret_mar2025 <- as.numeric(mar2025[1, tickers])
names(ret_mar2025) <- tickers
rp_capm_mar <- sum(w_mvp_capm * ret_mar2025)
rp_ff3_mar <- sum(w_mvp_ff3 * ret_mar2025)
cat(sprintf("Individual ETF returns in March 2025:\n"))## Individual ETF returns in March 2025:## SPY QQQ EEM IWM EFA TLT IYR GLD
## -0.0557 -0.0759 0.0113 -0.0685 0.0018 -0.0120 -0.0234 0.0945cat(sprintf("\nRealized MVP Return (CAPM weights) – March 2025: %.4f (%.2f%%)\n",
rp_capm_mar, rp_capm_mar * 100))##
## Realized MVP Return (CAPM weights) – March 2025: 0.0386 (3.86%)cat(sprintf("Realized MVP Return (FF3 weights) – March 2025: %.4f (%.2f%%)\n",
rp_ff3_mar, rp_ff3_mar * 100))## Realized MVP Return (FF3 weights) – March 2025: 0.0377 (3.77%)kable(data.frame(
Model = c("CAPM", "FF3"),
Return_Mar = c(rp_capm_mar, rp_ff3_mar),
Return_Pct = c(rp_capm_mar * 100, rp_ff3_mar * 100)
), digits = 4,
caption = "Realized MVP Portfolio Returns – March 2025") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Model | Return_Mar | Return_Pct |
|---|---|---|
| CAPM | 0.0386 | 3.8576 |
| FF3 | 0.0377 | 3.7730 |
Re-estimate MVP weights using the window 2020/04–2025/03, then apply to April 2025 returns.
# --- New 60-month window: 2020/04 – 2025/03 ---
window2_start <- as.Date("2020-04-01")
window2_end <- as.Date("2025-03-31")
window2_data <- merged_df %>%
filter(Date >= window2_start, Date <= window2_end)
cat("Window 2 observations:", nrow(window2_data), "\n")## Window 2 observations: 60etf_excess2 <- window2_data %>%
select(all_of(tickers)) %>%
sweep(1, window2_data$RF, "-")
mkt_excess2 <- window2_data$Mkt_RF
# --- CAPM MVP (rolling) ---
betas_capm2 <- sapply(tickers, function(tk) {
fit <- lm(etf_excess2[[tk]] ~ mkt_excess2)
coef(fit)[2]
})
resid_var_capm2 <- sapply(tickers, function(tk) {
fit <- lm(etf_excess2[[tk]] ~ mkt_excess2)
var(residuals(fit))
})
var_mkt2 <- var(mkt_excess2)
sigma_capm2 <- outer(betas_capm2, betas_capm2) * var_mkt2 + diag(resid_var_capm2)
colnames(sigma_capm2) <- rownames(sigma_capm2) <- tickers
w_mvp_capm2 <- solve_mvp(sigma_capm2)
# --- FF3 MVP (rolling) ---
B_ff3_2 <- sapply(tickers, function(tk) {
fit <- lm(etf_excess2[[tk]] ~ Mkt_RF + SMB + HML, data = window2_data)
coef(fit)[2:4]
})
B_ff3_2 <- t(B_ff3_2)
colnames(B_ff3_2) <- c("Mkt_RF", "SMB", "HML")
Sigma_F2 <- cov(window2_data[, c("Mkt_RF", "SMB", "HML")])
resid_var_ff32 <- sapply(tickers, function(tk) {
fit <- lm(etf_excess2[[tk]] ~ Mkt_RF + SMB + HML, data = window2_data)
var(residuals(fit))
})
sigma_ff3_2 <- B_ff3_2 %*% Sigma_F2 %*% t(B_ff3_2) + diag(resid_var_ff32)
colnames(sigma_ff3_2) <- rownames(sigma_ff3_2) <- tickers
w_mvp_ff3_2 <- solve_mvp(sigma_ff3_2)
# --- April 2025 realized returns ---
apr2025 <- merged_df %>%
filter(year(Date) == 2025, month(Date) == 4) %>%
select(all_of(tickers))
if (nrow(apr2025) == 0) {
prices_apr <- prices_monthly["2025-04"]
prices_mar2 <- prices_monthly["2025-03"]
apr2025_ret <- as.numeric((prices_apr - prices_mar2) / prices_mar2)
names(apr2025_ret) <- tickers
apr2025 <- as.data.frame(t(apr2025_ret))
}
ret_apr2025 <- as.numeric(apr2025[1, tickers])
names(ret_apr2025) <- tickers
rp_capm_apr <- sum(w_mvp_capm2 * ret_apr2025)
rp_ff3_apr <- sum(w_mvp_ff3_2 * ret_apr2025)
cat(sprintf("Individual ETF returns in April 2025:\n"))## Individual ETF returns in April 2025:## SPY QQQ EEM IWM EFA TLT IYR GLD
## -0.0087 0.0140 0.0014 -0.0232 0.0370 -0.0136 -0.0215 0.0542cat(sprintf("\nRealized MVP Return (CAPM weights) – April 2025: %.4f (%.2f%%)\n",
rp_capm_apr, rp_capm_apr * 100))##
## Realized MVP Return (CAPM weights) – April 2025: 0.0218 (2.18%)cat(sprintf("Realized MVP Return (FF3 weights) – April 2025: %.4f (%.2f%%)\n",
rp_ff3_apr, rp_ff3_apr * 100))## Realized MVP Return (FF3 weights) – April 2025: 0.0213 (2.13%)kable(data.frame(
Model = c("CAPM", "FF3"),
Return_Apr = c(rp_capm_apr, rp_ff3_apr),
Return_Pct = c(rp_capm_apr * 100, rp_ff3_apr * 100)
), digits = 4,
caption = "Realized MVP Portfolio Returns – April 2025 (Rolling Window)") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Model | Return_Apr | Return_Pct |
|---|---|---|
| CAPM | 0.0218 | 2.1839 |
| FF3 | 0.0213 | 2.1333 |
weights_summary <- data.frame(
Ticker = tickers,
CAPM_w_Q5 = round(w_mvp_capm, 4),
FF3_w_Q6 = round(w_mvp_ff3, 4),
CAPM_w_Q8roll = round(w_mvp_capm2, 4),
FF3_w_Q8roll = round(w_mvp_ff3_2, 4)
)
kable(weights_summary,
caption = "MVP Weights: Original Window vs Rolling Window") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Ticker | CAPM_w_Q5 | FF3_w_Q6 | CAPM_w_Q8roll | FF3_w_Q8roll | |
|---|---|---|---|---|---|
| SPY | SPY | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| QQQ | QQQ | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| EEM | EEM | 0.1401 | 0.1565 | 0.1847 | 0.1949 |
| IWM | IWM | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| EFA | EFA | 0.0838 | 0.0821 | 0.1140 | 0.1051 |
| TLT | TLT | 0.3425 | 0.3391 | 0.3046 | 0.3064 |
| IYR | IYR | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| GLD | GLD | 0.4336 | 0.4223 | 0.3967 | 0.3936 |
returns_summary <- data.frame(
Question = c("Q7 – March 2025", "Q8 – April 2025"),
CAPM_Ret = c(rp_capm_mar, rp_capm_apr),
FF3_Ret = c(rp_ff3_mar, rp_ff3_apr)
)
kable(returns_summary, digits = 4,
caption = "Realized MVP Portfolio Returns Summary") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Question | CAPM_Ret | FF3_Ret |
|---|---|---|
| Q7 – March 2025 | 0.0386 | 0.0377 |
| Q8 – April 2025 | 0.0218 | 0.0213 |
End of Part I. Part II textbook problems (Chapter 5–8) to be answered below.