- Import Data
library(fBasics)
library(kableExtra)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyquant)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## ── Attaching core tidyquant packages ─────────────────────── tidyquant 1.0.11 ──
## ✔ PerformanceAnalytics 2.0.8 ✔ TTR 0.24.4
## ✔ quantmod 0.4.27 ✔ xts 0.14.1
## ── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
## ✖ zoo::as.Date() masks base::as.Date()
## ✖ zoo::as.Date.numeric() masks base::as.Date.numeric()
## ✖ dplyr::filter() masks stats::filter()
## ✖ xts::first() masks dplyr::first()
## ✖ dplyr::group_rows() masks kableExtra::group_rows()
## ✖ PerformanceAnalytics::kurtosis() masks fBasics::kurtosis()
## ✖ dplyr::lag() masks stats::lag()
## ✖ xts::last() masks dplyr::last()
## ✖ PerformanceAnalytics::legend() masks graphics::legend()
## ✖ PerformanceAnalytics::skewness() masks fBasics::skewness()
## ✖ quantmod::summary() masks base::summary()
## ✖ TTR::volatility() masks fBasics::volatility()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
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(quantmod)
symbols <- c("SPY","QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
portfolioPrices <- NULL
for (symbol in symbols) {
portfolioPrices <- cbind(portfolioPrices,
getSymbols.yahoo(symbol,
from = '2010-01-01',
to = Sys.Date(),
auto.assign = FALSE)[,6])
}
colnames(portfolioPrices) <- symbols
head(portfolioPrices)
## SPY QQQ EEM IWM EFA TLT IYR
## 2010-01-04 85.76843 40.48581 31.08410 51.92010 36.38437 58.68239 27.40289
## 2010-01-05 85.99551 40.48581 31.30971 51.74158 36.41644 59.06138 27.46869
## 2010-01-06 86.05603 40.24160 31.37521 51.69289 36.57037 58.27078 27.45673
## 2010-01-07 86.41927 40.26779 31.19327 52.07430 36.42926 58.36877 27.70198
## 2010-01-08 86.70687 40.59919 31.44071 52.35833 36.71788 58.34262 27.51655
## 2010-01-11 86.82797 40.43350 31.37521 52.14734 37.01931 58.02246 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
###2. Calculate weekly and monthly returns using discrete returns
#Change data from daily to weekly
weekly_prices <- to.weekly(portfolioPrices, indexAt = 'last', OHLC = FALSE)
#Calculate weekly returns
weekly_returns <- na.omit(Return.calculate(weekly_prices), method = "discrete")
head(weekly_returns)
## SPY QQQ EEM IWM EFA
## 2010-01-15 -0.008117138 -0.015037068 -0.02893499 -0.01301956 -0.003493580
## 2010-01-22 -0.038982690 -0.036859441 -0.05578064 -0.03062165 -0.055740437
## 2010-01-29 -0.016665133 -0.031023614 -0.03357750 -0.02624320 -0.025802973
## 2010-02-05 -0.006797858 0.004440274 -0.02821291 -0.01397459 -0.019054795
## 2010-02-12 0.012938122 0.018147849 0.03333343 0.02952601 0.005245014
## 2010-02-19 0.028693133 0.024451694 0.02445359 0.03343184 0.022995024
## TLT IYR GLD
## 2010-01-15 2.004719e-02 -0.006304674 -0.004579349
## 2010-01-22 1.010117e-02 -0.041784732 -0.033285246
## 2010-01-29 3.368498e-03 -0.008447704 -0.011290465
## 2010-02-05 -5.451749e-05 0.003223355 -0.012080019
## 2010-02-12 -1.945950e-02 -0.007573790 0.022544905
## 2010-02-19 -8.205491e-03 0.050184938 0.022701796
#Change data from daily to monthly
monthly_prices <- to.monthly(portfolioPrices, indexAt = 'last', OHLC = FALSE)
#Calculate monthly returns
monthly_returns <- na.omit(Return.calculate(monthly_prices), method = "discrete")
head(monthly_returns)
## SPY QQQ EEM IWM EFA
## 2010-02-26 0.03119448 0.04603882 0.017764175 0.04475128 0.002667794
## 2010-03-31 0.06087953 0.07710885 0.081108991 0.08230691 0.063854092
## 2010-04-30 0.01547048 0.02242557 -0.001662405 0.05678464 -0.028045927
## 2010-05-28 -0.07945517 -0.07392389 -0.093935698 -0.07536617 -0.111927633
## 2010-06-30 -0.05174046 -0.05975694 -0.013986485 -0.07743433 -0.020619566
## 2010-07-30 0.06830100 0.07258231 0.109324916 0.06730944 0.116103750
## TLT IYR GLD
## 2010-02-26 -0.003423913 0.05457025 0.032748219
## 2010-03-31 -0.020573127 0.09748502 -0.004386396
## 2010-04-30 0.033218546 0.06388104 0.058834363
## 2010-05-28 0.051082630 -0.05683525 0.030513147
## 2010-06-30 0.057979415 -0.04670120 0.023553189
## 2010-07-30 -0.009463741 0.09404779 -0.050871157
###3. Download Fama French 3 factors data and change to digit numbers
# Load necessary libraries
library(readr)
library(dplyr)
# Unzip the correct file (it's a zip, not a CSV)
unzip("/cloud/project/F-F_Research_Data_Factors_daily.CSV", exdir = "/cloud/project")
## Warning in unzip("/cloud/project/F-F_Research_Data_Factors_daily.CSV", exdir =
## "/cloud/project"): error 1 in extracting from zip file
# Check the unzipped files
list.files("/cloud/project")
## [1] "F-F_Research_Data_Factors_daily.CSV"
## [2] "F-F_Research_Data_Factors_weekly.csv"
## [3] "mid term.Rmd"
## [4] "mid-term.Rmd"
## [5] "project.Rproj"
# Read the actual CSV file (you may need to adjust the name if it's different)
ff_data <- read_csv("/cloud/project/F-F_Research_Data_Factors_daily.CSV", skip = 3)
## New names:
## • `` -> `...1`
## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)
## Rows: 25902 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): ...1
## dbl (4): Mkt-RF, SMB, HML, RF
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Check column names to confirm what's there
colnames(ff_data)
## [1] "...1" "Mkt-RF" "SMB" "HML" "RF"
# Drop 'RF' column if it exists, otherwise just show column names and decide
FF3 <- ff_data %>% select(-RF)
# Convert percentages to decimals (skip the first column which is usually Date)
FF3 <- FF3 %>% mutate(across(-1, ~ . / 100))
# Preview result
head(FF3)
## # A tibble: 6 × 4
## ...1 `Mkt-RF` SMB HML
## <chr> <dbl> <dbl> <dbl>
## 1 19260701 0.001 -0.0025 -0.0027
## 2 19260702 0.0045 -0.0033 -0.0006
## 3 19260706 0.0017 0.003 -0.0039
## 4 19260707 0.0009 -0.0058 0.0002
## 5 19260708 0.0021 -0.0038 0.0019
## 6 19260709 -0.0071 0.0043 0.0057
###4. Merge monthly return data in question 2 and 3
library(dplyr)
library(lubridate)
library(xts)
# 1. Ensure FF3 date column is properly formatted
# Rename the first column to 'time' if it's unnamed
names(FF3)[1] <- "time"
# Convert 'time' column to Date (from 'YYYYMMDD' to Date format)
FF3$time <- as.Date(FF3$time, format = "%Y%m%d")
# 2. Filter FF3 between Feb 2010 and Feb 2024
FF3_subset <- FF3 %>%
filter(time >= as.Date("2010-02-01") & time <= as.Date("2024-02-29"))
# 3. Subset the monthly returns
monthly_returns_subset <- monthly_returns["2010-02-26/2024-02-29"]
# 4. Convert xts to tibble with date
monthly_returns_subset_df <- as_tibble(index(monthly_returns_subset)) %>%
rename(time = value) %>%
bind_cols(as_tibble(coredata(monthly_returns_subset)))
# 5. Merge the datasets on 'time'
merged_monthlyreturns <- left_join(
monthly_returns_subset_df,
FF3_subset[, c("time", "Mkt-RF", "SMB", "HML")],
by = "time"
)
# 6. Convert to xts
merged_monthlyreturns_xts <- xts(merged_monthlyreturns[, -1], order.by = merged_monthlyreturns$time)
# 7. Preview
head(merged_monthlyreturns_xts)
## SPY QQQ EEM IWM EFA
## 2010-02-26 0.03119448 0.04603882 0.017764175 0.04475128 0.002667794
## 2010-03-31 0.06087953 0.07710885 0.081108991 0.08230691 0.063854092
## 2010-04-30 0.01547048 0.02242557 -0.001662405 0.05678464 -0.028045927
## 2010-05-28 -0.07945517 -0.07392389 -0.093935698 -0.07536617 -0.111927633
## 2010-06-30 -0.05174046 -0.05975694 -0.013986485 -0.07743433 -0.020619566
## 2010-07-30 0.06830100 0.07258231 0.109324916 0.06730944 0.116103750
## TLT IYR GLD Mkt-RF SMB HML
## 2010-02-26 -0.003423913 0.05457025 0.032748219 0.0013 -0.0057 0.0042
## 2010-03-31 -0.020573127 0.09748502 -0.004386396 -0.0036 -0.0051 0.0020
## 2010-04-30 0.033218546 0.06388104 0.058834363 -0.0172 -0.0097 -0.0088
## 2010-05-28 0.051082630 -0.05683525 0.030513147 -0.0118 0.0010 -0.0068
## 2010-06-30 0.057979415 -0.04670120 0.023553189 -0.0098 0.0005 -0.0035
## 2010-07-30 -0.009463741 0.09404779 -0.050871157 0.0006 0.0005 -0.0021
###5. Based on CAPM model, compute MVP monthly returns based on estimated covariance matrix for the 8-asset portfolio by using past 60-month returns from 2019/03 - 2024/02
monthly_returns_60 <- merged_monthlyreturns_xts["2019-03-29/2024-02-29"]
portfolio_8 <- monthly_returns_60[, 1:8]
rf <- monthly_returns_60[, 9]
#Find covariance matrix using CAMP
cov_ma_sim <- function(returns, rf){
n <- nrow(returns)
X <- as.matrix(cbind(rep(1, n), rf))
Y <- as.matrix(returns)
b_hat <- solve(t(X) %*% X) %*% t(X) %*% Y
E_hat = Y - X %*% b_hat
b_hat <- as.matrix(b_hat[-1, ])
diagD_hat <- diag(t(E_hat) %*% E_hat) / (n - 2)
diag(diagD_hat)
cov_sfm <- as.numeric(var(rf)) * b_hat %*% t(b_hat) + diag(diagD_hat)
return(cov_sfm)
}
cov_ma <- cov_ma_sim(portfolio_8, rf)
head(cov_ma)
## SPY QQQ EEM IWM EFA
## SPY 0.0028372300 0.0004706850 0.0006305482 0.0003710877 0.0005480242
## QQQ 0.0004706850 0.0040445581 0.0007325720 0.0004311304 0.0006366955
## EEM 0.0006305482 0.0007325720 0.0030581809 0.0005775592 0.0008529423
## IWM 0.0003710877 0.0004311304 0.0005775592 0.0048295825 0.0005019702
## EFA 0.0005480242 0.0006366955 0.0008529423 0.0005019702 0.0027448685
## TLT 0.0002242704 0.0002605578 0.0003490535 0.0002054235 0.0003033705
## TLT IYR GLD
## SPY 0.0002242704 0.0003645923 0.0002493140
## QQQ 0.0002605578 0.0004235840 0.0002896534
## EEM 0.0003490535 0.0005674499 0.0003880311
## IWM 0.0002054235 0.0003339533 0.0002283625
## EFA 0.0003033705 0.0004931840 0.0003372469
## TLT 0.0022592366 0.0002018279 0.0001380131
#Find MPV weights
MVP_sim <- function(cov_ma){
x <- ncol(cov_ma)
one_ma <- matrix(rep(1,x), ncol = 1)
numerator <- inv(cov_ma) %*% one_ma
denominator <- t(one_ma) %*% inv(cov_ma) %*% one_ma
weights <- numerator/as.vector(denominator)
colnames(weights) <- "Weight"
return(weights)
}
capm_model<- MVP_sim(cov_ma)
capm_model
## Weight
## SPY 0.12480450
## QQQ 0.07071843
## EEM 0.05312760
## IWM 0.07410709
## EFA 0.09037234
## TLT 0.21477027
## IYR 0.10345745
## GLD 0.26864232
###6. Based on FF 3-factor model, compute MVP monthly returns covariance matrix for the 8-asset portfolio by using past 60-month returns from 2019/03 - 2024/02
#Find covariance matrix using Fama-French 3-factor
N <- nrow(portfolio_8)
ones <- rep(1, N)
portfolio_ma <- as.matrix(portfolio_8)
sigF3 <- as.matrix(var(cbind(monthly_returns_60$`Mkt-RF`,
monthly_returns_60$SMB,
monthly_returns_60$HML)))
X.3 <- cbind(ones, monthly_returns_60$`Mkt-RF`, monthly_returns_60$SMB, monthly_returns_60$HML)
b_hat.3 <- solve(t(X.3) %*% (X.3)) %*% t(X.3) %*% portfolio_ma
E_hat.3 <- portfolio_ma - X.3 %*% b_hat.3
b_hat.3 <- as.matrix(b_hat.3[-1, ])
diagD_hat.3 <- diag(t(E_hat.3) %*% E_hat.3)/(N-4)
cov_3f.3 = t(b_hat.3) %*% sigF3 %*% b_hat.3 + diag(diagD_hat.3)
cov_3f.3
## SPY QQQ EEM IWM EFA
## SPY 0.0028934017 0.0011043494 0.0010716588 0.0011131009 0.0010179418
## QQQ 0.0011043494 0.0041318726 0.0012182462 0.0012460517 0.0011567139
## EEM 0.0010716588 0.0012182462 0.0031138867 0.0011619540 0.0012142948
## IWM 0.0011131009 0.0012460517 0.0011619540 0.0049366304 0.0011093705
## EFA 0.0010179418 0.0011567139 0.0012142948 0.0011093705 0.0027957711
## TLT 0.0004713182 0.0005320650 0.0005432111 0.0005355972 0.0005056403
## IYR 0.0009496624 0.0010662937 0.0010279004 0.0011175376 0.0009721463
## GLD 0.0003564138 0.0004046295 0.0004782744 0.0003872371 0.0004243401
## TLT IYR GLD
## SPY 0.0004713182 0.0009496624 0.0003564138
## QQQ 0.0005320650 0.0010662937 0.0004046295
## EEM 0.0005432111 0.0010279004 0.0004782744
## IWM 0.0005355972 0.0011175376 0.0003872371
## EFA 0.0005056403 0.0009721463 0.0004243401
## TLT 0.0023296179 0.0004619018 0.0001902717
## IYR 0.0004619018 0.0036816141 0.0003491122
## GLD 0.0001902717 0.0003491122 0.0018368439
#Find MVP weights
ff3_model <- MVP_sim(cov_3f.3)
ff3_model
## Weight
## SPY 0.10562485
## QQQ 0.04081132
## EEM 0.04987760
## IWM 0.03384879
## EFA 0.08699072
## TLT 0.25130397
## IYR 0.07972748
## GLD 0.35181528
###7. You can invest in the 8-portfolio in 2024/03 based on the optimal weights of MVP from question 5 and 6. What are the realized portfolio returns in the March of 2024 using the weights from question 5 and 6?
march_return <- monthly_returns['2024-03-28']
capm_return <- march_return %*% capm_model
capm_return <- as.data.frame(capm_return)
colnames(capm_return) <- "Return"
capm_return <- capm_return %>% mutate(Return = scales::percent(Return, accuracy = 0.001))
capm_return %>% kbl(format = "html", caption = "Portfolio Return in March using CAMP Model")
Portfolio Return in March using CAMP Model
|
Return
|
|
3.895%
|
ff3_return <- march_return %*% ff3_model
ff3_return <- as.data.frame(ff3_return)
colnames(ff3_return) <- "Return"
ff3_return <- ff3_return %>% mutate(Return = scales::percent(Return, accuracy = 0.001))
ff3_return %>% kbl(format = "html", caption = "Portfolio Return in March using Fama-French Model")
Portfolio Return in March using Fama-French Model
|
Return
|
|
4.340%
|