# Part I: Computer Questions (40%)

## Q1 — Download ETF Daily Data from Yahoo Finance
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(xts)
library(zoo)
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library(quadprog)
library(frenchdata)

# Define the 8 ETF tickers
tickers <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")

# Download adjusted daily 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 only
prices_list <- lapply(tickers, function(t) Ad(get(t)))

# Merge into a single xts object
prices      <- do.call(merge, prices_list)
colnames(prices) <- tickers

# Quick check
cat("Dimensions:", dim(prices), "\n")
## Dimensions: 4023 8
cat("Period:", as.character(start(prices)), "to", as.character(end(prices)), "\n")
## Period: 2010-01-04 to 2025-12-30
head(prices)
##                 SPY      QQQ      EEM      IWM      EFA      TLT      IYR
## 2010-01-04 84.79637 40.29078 30.35151 51.36656 35.12845 56.13522 26.76812
## 2010-01-05 85.02085 40.29078 30.57181 51.18993 35.15940 56.49773 26.83238
## 2010-01-06 85.08070 40.04777 30.63577 51.14177 35.30801 55.74139 26.82069
## 2010-01-07 85.43986 40.07379 30.45810 51.51908 35.17178 55.83514 27.06028
## 2010-01-08 85.72421 40.40364 30.69972 51.80008 35.45044 55.81014 26.87913
## 2010-01-11 85.84387 40.23873 30.63577 51.59138 35.74147 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

Q2 — Monthly Discrete Returns

# Convert daily prices to monthly (last trading day of each month)
prices_monthly <- to.monthly(prices, indexAt = "lastof", OHLC = FALSE)

# Compute discrete returns: (P_t - P_{t-1}) / P_{t-1}
returns_monthly <- na.omit(Return.calculate(prices_monthly, method = "discrete"))

# Convert to data frame
returns_df      <- data.frame(
  Date = as.yearmon(index(returns_monthly)),
  coredata(returns_monthly)
)
colnames(returns_df) <- c("Date", tickers)

cat("Monthly returns dimensions:", dim(returns_df), "\n")
## Monthly returns dimensions: 191 9
head(returns_df)
##       Date         SPY         QQQ          EEM         IWM          EFA
## 1 Feb 2010  0.03119460  0.04603816  0.017763983  0.04475120  0.002667738
## 2 Mar 2010  0.06087975  0.07710959  0.081108354  0.08230661  0.063853603
## 3 Apr 2010  0.01547007  0.02242490 -0.001661876  0.05678456 -0.028045564
## 4 May 2010 -0.07945473 -0.07392331 -0.093935709 -0.07536618 -0.111927733
## 5 Jun 2010 -0.05174083 -0.05975696 -0.013986500 -0.07743393 -0.020619794
## 6 Jul 2010  0.06830105  0.07258271  0.109324939  0.06730918  0.116104341
##            TLT         IYR          GLD
## 1 -0.003424037  0.05457078  0.032748219
## 2 -0.020574209  0.09748448 -0.004386396
## 3  0.033218575  0.06388122  0.058834363
## 4  0.051083319 -0.05683540  0.030513147
## 5  0.057978580 -0.04670136  0.023553189
## 6 -0.009463404  0.09404841 -0.050871157

Q3 — Download Fama-French 3 Factors

# Download FF3 monthly data
ff3_raw     <- download_french_data("Fama/French 3 Factors")
## New names:
## New names:
## • `` -> `...1`
ff3_monthly <- ff3_raw$subsets$data[[1]]

# Rename columns
colnames(ff3_monthly)[colnames(ff3_monthly) == "date"]   <- "Date"
colnames(ff3_monthly)[colnames(ff3_monthly) == "Mkt-RF"] <- "Mkt_RF"

# Convert from percentage to decimal
ff3_monthly$Date   <- as.yearmon(as.character(ff3_monthly$Date), "%Y%m")
ff3_monthly$Mkt_RF <- as.numeric(ff3_monthly$Mkt_RF) / 100
ff3_monthly$SMB    <- as.numeric(ff3_monthly$SMB)    / 100
ff3_monthly$HML    <- as.numeric(ff3_monthly$HML)    / 100
ff3_monthly$RF     <- as.numeric(ff3_monthly$RF)     / 100

# Filter date range
ff3 <- ff3_monthly[ff3_monthly$Date >= as.yearmon("2010-01") &
                   ff3_monthly$Date <= as.yearmon("2025-04"), ]


cat("FF3 dimensions:", dim(ff3), "\n")
## FF3 dimensions: 184 5
head(ff3)
## # A tibble: 6 × 5
##   Date       Mkt_RF     SMB     HML     RF
##   <yearmon>   <dbl>   <dbl>   <dbl>  <dbl>
## 1 Jan 2010  -0.0335  0.0043  0.0033 0     
## 2 Feb 2010   0.0339  0.0118  0.0318 0     
## 3 Mar 2010   0.063   0.0146  0.0219 0.0001
## 4 Apr 2010   0.0199  0.0484  0.0296 0.0001
## 5 May 2010  -0.079   0.0013 -0.0248 0.0001
## 6 Jun 2010  -0.0556 -0.0179 -0.0473 0.0001

Q4 — Merge Monthly Returns with Fama-French Factors

# Merge ETF returns and FF3 factors by Date
merged_df <- merge(returns_df, ff3, by = "Date")

cat("Merged dataset dimensions:", dim(merged_df), "\n")
## Merged dataset dimensions: 183 13
head(merged_df)
##       Date         SPY         QQQ          EEM         IWM          EFA
## 1 Feb 2010  0.03119460  0.04603816  0.017763983  0.04475120  0.002667738
## 2 Mar 2010  0.06087975  0.07710959  0.081108354  0.08230661  0.063853603
## 3 Apr 2010  0.01547007  0.02242490 -0.001661876  0.05678456 -0.028045564
## 4 May 2010 -0.07945473 -0.07392331 -0.093935709 -0.07536618 -0.111927733
## 5 Jun 2010 -0.05174083 -0.05975696 -0.013986500 -0.07743393 -0.020619794
## 6 Jul 2010  0.06830105  0.07258271  0.109324939  0.06730918  0.116104341
##            TLT         IYR          GLD  Mkt_RF     SMB     HML    RF
## 1 -0.003424037  0.05457078  0.032748219  0.0339  0.0118  0.0318 0e+00
## 2 -0.020574209  0.09748448 -0.004386396  0.0630  0.0146  0.0219 1e-04
## 3  0.033218575  0.06388122  0.058834363  0.0199  0.0484  0.0296 1e-04
## 4  0.051083319 -0.05683540  0.030513147 -0.0790  0.0013 -0.0248 1e-04
## 5  0.057978580 -0.04670136  0.023553189 -0.0556 -0.0179 -0.0473 1e-04
## 6 -0.009463404  0.09404841 -0.050871157  0.0692  0.0022 -0.0050 1e-04
# Estimation window: 60 months ending Feb 2020
window_data <- merged_df[merged_df$Date >= as.yearmon("2015-03") &
                         merged_df$Date <= as.yearmon("2020-02"), ]

# ETF returns matrix
etf_returns <- as.matrix(window_data[, tickers])
rf_vec      <- window_data$RF
mkt         <- window_data$Mkt_RF
n           <- length(tickers)

# Excess returns: R_i - RF
excess_ret  <- etf_returns - rf_vec

# CAPM residuals for each ETF
resid_capm  <- matrix(NA, nrow = nrow(window_data), ncol = n)
colnames(resid_capm) <- tickers

for (i in 1:n) {
  fit             <- lm(excess_ret[, i] ~ mkt)
  resid_capm[, i] <- residuals(fit)
}

# Covariance matrix
cov_capm <- cov(resid_capm)

# MVP optimization
Amat <- cbind(rep(1, n), diag(n))
bvec <- c(1, rep(0, n))

weights_capm        <- solve.QP(2 * cov_capm, rep(0, n), Amat, bvec, meq = 1)$solution
names(weights_capm) <- tickers

cat("=== MVP Weights (CAPM) ===\n")
## === MVP Weights (CAPM) ===
print(round(weights_capm, 4))
##    SPY    QQQ    EEM    IWM    EFA    TLT    IYR    GLD 
## 0.8149 0.0294 0.0007 0.1453 0.0000 0.0000 0.0000 0.0097
cat("Sum of weights:", round(sum(weights_capm), 6), "\n")
## Sum of weights: 1
var_capm <- t(weights_capm) %*% cov_capm %*% weights_capm
cat("Monthly portfolio volatility (CAPM):", round(sqrt(var_capm) * 100, 4), "%\n")
## Monthly portfolio volatility (CAPM): 0.177 %

Q6 — MVP via Fama-French 3-Factor Model (same 60-month window)

smb <- window_data$SMB
hml <- window_data$HML

# Compute FF3 residuals for each ETF
resid_ff3 <- matrix(NA, nrow = nrow(window_data), ncol = n)
colnames(resid_ff3) <- tickers

for (i in 1:n) {
  fit             <- lm(excess_ret[, i] ~ mkt + smb + hml)
  resid_ff3[, i] <- residuals(fit)
}

# Covariance matrix from FF3 residuals
cov_ff3 <- cov(resid_ff3)

# MVP optimization
weights_ff3        <- solve.QP(2 * cov_ff3, rep(0, n), Amat, bvec, meq = 1)$solution
names(weights_ff3) <- tickers

cat("=== MVP Weights (FF3) ===\n")
## === MVP Weights (FF3) ===
print(round(weights_ff3, 4))
##    SPY    QQQ    EEM    IWM    EFA    TLT    IYR    GLD 
## 0.7033 0.0343 0.0045 0.2452 0.0000 0.0000 0.0000 0.0127
cat("Sum of weights:", round(sum(weights_ff3), 6), "\n")
## Sum of weights: 1
var_ff3 <- t(weights_ff3) %*% cov_ff3 %*% weights_ff3
cat("Monthly portfolio volatility (FF3):", round(sqrt(var_ff3) * 100, 4), "%\n")
## Monthly portfolio volatility (FF3): 0.1694 %

Q7 — Realized Portfolio Returns in March 2025

# Actual ETF returns in March 2025
mar2025_row <- merged_df[merged_df$Date == as.yearmon("2025-03"), ]
ret_mar2025 <- as.numeric(mar2025_row[, tickers])
names(ret_mar2025) <- tickers

cat("ETF Returns - March 2025:\n")
## ETF Returns - March 2025:
print(round(ret_mar2025, 4))
##     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 MVP return
ret_mvp_capm_mar <- sum(weights_capm * ret_mar2025)
ret_mvp_ff3_mar  <- sum(weights_ff3  * ret_mar2025)

cat("\nRealized MVP Return - CAPM (March 2025):", round(ret_mvp_capm_mar * 100, 4), "%\n")
## 
## Realized MVP Return - CAPM (March 2025): -5.6676 %
cat("Realized MVP Return - FF3  (March 2025):", round(ret_mvp_ff3_mar  * 100, 4), "%\n")
## Realized MVP Return - FF3  (March 2025): -5.7349 %

Q8 — Realized MVP Returns in April 2025 (rolling 60-month window)

# New estimation window: 2020/04 - 2025/03
window_apr <- merged_df[merged_df$Date >= as.yearmon("2020-04") & merged_df$Date <= as.yearmon("2025-03"), ]

etf_ret_apr <- as.matrix(window_apr[, tickers])
rf_apr      <- window_apr$RF
mkt_apr     <- window_apr$Mkt_RF
smb_apr     <- window_apr$SMB
hml_apr     <- window_apr$HML
excess_apr  <- etf_ret_apr - rf_apr

# CAPM residuals
resid_capm_apr <- matrix(NA, nrow = nrow(window_apr), ncol = n)
for (i in 1:n) {
  fit                  <- lm(excess_apr[, i] ~ mkt_apr)
  resid_capm_apr[, i] <- residuals(fit)
}
cov_capm_apr <- cov(resid_capm_apr)

# FF3 residuals
resid_ff3_apr <- matrix(NA, nrow = nrow(window_apr), ncol = n)
for (i in 1:n) {
  fit                 <- lm(excess_apr[, i] ~ mkt_apr + smb_apr + hml_apr)
  resid_ff3_apr[, i] <- residuals(fit)
}
cov_ff3_apr <- cov(resid_ff3_apr)

# New MVP weights
w_capm_apr        <- solve.QP(2*cov_capm_apr, rep(0,n), Amat, bvec, meq=1)$solution
w_ff3_apr         <- solve.QP(2*cov_ff3_apr,  rep(0,n), Amat, bvec, meq=1)$solution
names(w_capm_apr) <- tickers
names(w_ff3_apr)  <- tickers

cat("=== MVP Weights CAPM - April 2025 ===\n")
## === MVP Weights CAPM - April 2025 ===
print(round(w_capm_apr, 4))
##    SPY    QQQ    EEM    IWM    EFA    TLT    IYR    GLD 
## 0.7497 0.1019 0.0060 0.1410 0.0000 0.0000 0.0000 0.0014
cat("\n=== MVP Weights FF3 - April 2025 ===\n")
## 
## === MVP Weights FF3 - April 2025 ===
print(round(w_ff3_apr, 4))
##    SPY    QQQ    EEM    IWM    EFA    TLT    IYR    GLD 
## 0.6436 0.0840 0.0088 0.2636 0.0000 0.0000 0.0000 0.0000
# Actual ETF returns in April 2025
apr2025_row <- merged_df[merged_df$Date == as.yearmon("2025-04"), ]
ret_apr2025 <- as.numeric(apr2025_row[, tickers])
names(ret_apr2025) <- tickers

cat("\nETF Returns - April 2025:\n")
## 
## ETF Returns - April 2025:
print(round(ret_apr2025, 4))
##     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_mvp_capm_apr <- sum(w_capm_apr * ret_apr2025)
ret_mvp_ff3_apr  <- sum(w_ff3_apr  * ret_apr2025)

cat("\nRealized MVP Return - CAPM (April 2025):", round(ret_mvp_capm_apr * 100, 4), "%\n")
## 
## Realized MVP Return - CAPM (April 2025): -0.8266 %
cat("Realized MVP Return - FF3  (April 2025):", round(ret_mvp_ff3_apr  * 100, 4), "%\n")
## Realized MVP Return - FF3  (April 2025): -1.0512 %