library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(reshape2)
##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smiths
# Read and inspect data
retdata = read_csv('FamaFrench_mon_69_98_3stocks.csv')
## Rows: 360 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (9): date, Mkt-RF, SMB, HML, RF, ge, ibm, mobil, CRSP
##
## ℹ 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.
head(retdata)
## # A tibble: 6 × 9
## date `Mkt-RF` SMB HML RF ge ibm mobil CRSP
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 196901 -1.2 -0.8 1.57 0.53 -1.20 -5.95 -1.40 -0.671
## 2 196902 -5.82 -3.9 0.93 0.46 -6.04 -0.700 -7.84 -5.36
## 3 196903 2.59 -0.28 -0.45 0.46 6.65 7.03 21.5 3.05
## 4 196904 1.52 -0.85 0.06 0.53 5.96 4.46 3.00 2.05
## 5 196905 0.02 -0.27 0.74 0.48 -3.58 -2.5 2.67 0.504
## 6 196906 -7.25 -5.31 -1.15 0.51 -3.82 5.88 -13.0 -6.74
glimpse(retdata)
## Rows: 360
## Columns: 9
## $ date <dbl> 196901, 196902, 196903, 196904, 196905, 196906, 196907, 19690…
## $ `Mkt-RF` <dbl> -1.20, -5.82, 2.59, 1.52, 0.02, -7.25, -7.05, 4.65, -2.88, 4.…
## $ SMB <dbl> -0.80, -3.90, -0.28, -0.85, -0.27, -5.31, -3.27, 0.89, 1.20, …
## $ HML <dbl> 1.57, 0.93, -0.45, 0.06, 0.74, -1.15, 1.36, -3.83, -3.24, -3.…
## $ RF <dbl> 0.53, 0.46, 0.46, 0.53, 0.48, 0.51, 0.53, 0.50, 0.62, 0.60, 0…
## $ ge <dbl> -1.1984, -6.0377, 6.6474, 5.9621, -3.5806, -3.8196, -4.3056, …
## $ ibm <dbl> -5.9524, -0.7004, 7.0303, 4.4586, -2.5000, 5.8777, -3.9230, 6…
## $ mobil <dbl> -1.4043, -7.8431, 21.5130, 2.9961, 2.6667, -12.9870, -6.0981,…
## $ CRSP <dbl> -0.6714, -5.3641, 3.0505, 2.0528, 0.5038, -6.7388, -6.5173, 5…
colnames(retdata)[2] <- 'Mkt_RF' # Replace 'Mkt-RF' with 'Mkt_RF'
# Single index model to compute covariance matrix
stock.rets <- retdata %>% select(c(2, 6, 7, 8)) / 100
glimpse(stock.rets)
## Rows: 360
## Columns: 4
## $ Mkt_RF <dbl> -0.0120, -0.0582, 0.0259, 0.0152, 0.0002, -0.0725, -0.0705, 0.0…
## $ ge <dbl> -0.011984, -0.060377, 0.066474, 0.059621, -0.035806, -0.038196,…
## $ ibm <dbl> -0.059524, -0.007004, 0.070303, 0.044586, -0.025000, 0.058777, …
## $ mobil <dbl> -0.014043, -0.078431, 0.215130, 0.029961, 0.026667, -0.129870, …
N <- dim(stock.rets)[1]
# Method 1: by "lm" function
fit = lm(formula = cbind(ge, ibm, mobil) ~ Mkt_RF, data = stock.rets)
sigF = as.numeric(var(stock.rets$Mkt_RF))
bbeta = as.matrix(fit$coefficients)
bbeta = as.matrix(bbeta[-1,])
sigeps = crossprod(fit$residuals) / (N - 2)
sigeps = diag(diag(sigeps))
cov_1f = sigF * bbeta %*% t(bbeta) + sigeps
# Method 2: by formula "inv(X'X)*X'Y"
ones = rep(1, N)
X = as.matrix(cbind(ones, stock.rets$Mkt_RF))
retdata1 = as.matrix(retdata[, c(6, 7, 8)] / 100)
b_hat = solve(t(X) %*% X) %*% t(X) %*% retdata1
E_hat = retdata1 - X %*% b_hat
b_hat = as.matrix(b_hat[-1,])
diagD_hat = diag(t(E_hat) %*% E_hat) / (N - 2)
cov_1f.1 = as.numeric(var(stock.rets$Mkt_RF)) * b_hat %*% t(b_hat) + diag(diagD_hat)
# Using FF 3 factor model to compute covariance matrix
N <- dim(retdata)[1]
stock.rets <- retdata %>% select(c(2, 3, 4, 6, 7, 8)) / 100
fit3 = lm(formula = cbind(ge, ibm, mobil) ~ Mkt_RF + SMB + HML, data = stock.rets)
sigF3 = as.matrix(var(cbind(stock.rets$Mkt_RF, stock.rets$SMB, stock.rets$HML)))
bbeta3 = as.matrix(fit3$coefficients)
bbeta3 = bbeta3[-1,]
sigeps3 = crossprod(fit3$residuals) / (N - 4)
sigeps3 = diag(diag(sigeps3))
cov_3f = t(bbeta3) %*% sigF3 %*% (bbeta3) + sigeps3
# Method 2: by formula "inv(X'X)*X'Y"
X.3 = cbind(ones, stock.rets$Mkt_RF, stock.rets$SMB, stock.rets$HML)
b_hat.3 = solve(t(X.3) %*% (X.3)) %*% t(X.3) %*% retdata1
E_hat.3 = retdata1 - 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)
# Create frontier function to plot efficient frontier
frontier <- function(return, Q) {
n = ncol(return)
Ax <- rbind(2 * cov(return), colMeans(return), rep(1, n))
Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2)))
r <- colMeans(return)
rbase <- seq(min(r), max(r), length = 100)
s <- sapply(rbase, function(x) {
b0 <- c(rep(0, ncol(return)), x, 1)
y <- head(solve(Ax, b0), n)
sqrt(y %*% Q %*% y)
})
efficient.port <- list("er" = as.vector(rbase),
"sd" = as.vector(s))
class(efficient.port) <- "portfolio"
efficient.port
}
# Use different covariance matrix to plot efficient frontier: Q
Q.3f = cov_3f
Q.1f = cov_1f
Q = cov(retdata1)
# Draw overlay frontiers on the same graph
xy.3f = frontier(retdata1, Q.3f)
xy.1f = frontier(retdata1, Q.1f)
xy = frontier(retdata1, Q)
# Convert to tibble and rename column names
xx <- cbind(xy$sd, xy.1f$sd, xy.3f$sd) %>%
as_tibble() %>%
rename(s = V1, s1 = V2, s3 = V3)
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
yy <- cbind(xy$er, xy.1f$er, xy.3f$er) %>%
as_tibble() %>%
rename(er = V1, er1 = V2, er3 = V3)
xy.all <- bind_cols(xx, yy)
# Using ggplot
ggplot(data = xy.all, aes(x = s, y = er)) +
geom_point(color = "red", size = 0.3) +
annotate(geom="text", x=0.045, y=0.0125, label="historical covariance",
color="red", size=4) +
geom_point(aes(x = s1, y = er1), color = "blue", size = 0.3, shape = 2) +
annotate(geom="text", x=0.058, y=0.013, label="capm covariance",
color="blue", size=4) +
geom_point(aes(x = s3, y = er3), color = "black", size = 0.3, shape = 4) +
annotate(geom="text", x=0.045, y=0.0145, label="FF3F covariance",
color="black", size=4) +
labs(title="Efficient Frontiers", x="sd", y="ret")
