install.packages("tidyverse")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
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
install.packages("dplyr")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.4'
## (as 'lib' is unspecified)
library(dplyr)
data <- read.csv("FamaFrench_mon_69_98_3stocks.csv")
head(data)
## date Mkt.RF SMB HML RF ge ibm mobil CRSP
## 1 196901 -1.20 -0.80 1.57 0.53 -1.1984 -5.9524 -1.4043 -0.6714
## 2 196902 -5.82 -3.90 0.93 0.46 -6.0377 -0.7004 -7.8431 -5.3641
## 3 196903 2.59 -0.28 -0.45 0.46 6.6474 7.0303 21.5130 3.0505
## 4 196904 1.52 -0.85 0.06 0.53 5.9621 4.4586 2.9961 2.0528
## 5 196905 0.02 -0.27 0.74 0.48 -3.5806 -2.5000 2.6667 0.5038
## 6 196906 -7.25 -5.31 -1.15 0.51 -3.8196 5.8777 -12.9870 -6.7388
glimpse(data)
## Rows: 360
## Columns: 9
## $ date <int> 196901, 196902, 196903, 196904, 196905, 196906, 196907, 196908,…
## $ Mkt.RF <dbl> -1.20, -5.82, 2.59, 1.52, 0.02, -7.25, -7.05, 4.65, -2.88, 4.96…
## $ SMB <dbl> -0.80, -3.90, -0.28, -0.85, -0.27, -5.31, -3.27, 0.89, 1.20, 3.…
## $ HML <dbl> 1.57, 0.93, -0.45, 0.06, 0.74, -1.15, 1.36, -3.83, -3.24, -3.19…
## $ RF <dbl> 0.53, 0.46, 0.46, 0.53, 0.48, 0.51, 0.53, 0.50, 0.62, 0.60, 0.5…
## $ ge <dbl> -1.1984, -6.0377, 6.6474, 5.9621, -3.5806, -3.8196, -4.3056, -2…
## $ ibm <dbl> -5.9524, -0.7004, 7.0303, 4.4586, -2.5000, 5.8777, -3.9230, 6.6…
## $ mobil <dbl> -1.4043, -7.8431, 21.5130, 2.9961, 2.6667, -12.9870, -6.0981, 1…
## $ CRSP <dbl> -0.6714, -5.3641, 3.0505, 2.0528, 0.5038, -6.7388, -6.5173, 5.1…
colnames(data)[2]<- 'Mkt_RF'
stock.rets<-data %>% 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]
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,])
bbeta
## [,1]
## ge 1.0580825
## ibm 0.8149949
## mobil 0.8158072
sigeps = crossprod(fit$residuals)/(N-2)
sigeps = diag(diag(sigeps))
sigeps
## [,1] [,2] [,3]
## [1,] 0.001702494 0.000000000 0.000000000
## [2,] 0.000000000 0.003225874 0.000000000
## [3,] 0.000000000 0.000000000 0.002913458
cov_1f = sigF*bbeta%*%t(bbeta)+sigeps
cov_1f
## ge ibm mobil
## ge 0.004070402 0.001823896 0.001825714
## ibm 0.001823896 0.004630742 0.001406268
## mobil 0.001825714 0.001406268 0.004321127
ones = rep(1,N)
X = as.matrix(cbind(ones, stock.rets$Mkt_RF))
data1 = as.matrix(data[,c(6,7,8)]/100)
b_hat = solve(t(X)%*%X)%*%t(X)%*%data1
E_hat = data1 - 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);
cov_1f.1
## ge ibm mobil
## ge 0.004070402 0.001823896 0.001825714
## ibm 0.001823896 0.004630742 0.001406268
## mobil 0.001825714 0.001406268 0.004321127
N <- dim(data)[1]
stock.rets<-data %>% 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,]
bbeta3
## ge ibm mobil
## Mkt_RF 1.134331448 0.8050676 0.9803403
## SMB -0.369412445 -0.3099907 -0.3727822
## HML 0.009630701 -0.2981744 0.3726475
sigeps3 = crossprod(fit3$residuals)/(N-4)
sigeps3 = diag(diag(sigeps3))
cov_3f = t(bbeta3) %*% sigF3 %*% (bbeta3) + sigeps3
cov_3f
## ge ibm mobil
## ge 0.004079122 0.001905590 0.001929618
## ibm 0.001905590 0.004647834 0.001424956
## mobil 0.001929618 0.001424956 0.004335927
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)%*%data1
E_hat.3 = data1 - 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
##
## ge 4.332517e-03 -1.258773e-04 -4.819330e-06
## ibm -1.258773e-04 3.199233e-03 -1.196042e-05
## mobil -4.819330e-06 -1.196042e-05 2.841930e-03
frontier <- function(return, Q) {
#return <- log(tail(assets, -1) / head(assets, -1))
n = ncol(return)
#Q = cov(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("call" = call,
"er" = as.vector(rbase),
"sd" = as.vector(s))
class(efficient.port) <- "portfolio"
efficient.port
}
Q.3f = cov_3f
Q.1f = cov_1f
head(data1)
## ge ibm mobil
## [1,] -0.011984 -0.059524 -0.014043
## [2,] -0.060377 -0.007004 -0.078431
## [3,] 0.066474 0.070303 0.215130
## [4,] 0.059621 0.044586 0.029961
## [5,] -0.035806 -0.025000 0.026667
## [6,] -0.038196 0.058777 -0.129870
Q = cov(data1)
xy.3f = frontier(data1, Q.3f)
xy.1f = frontier(data1, Q.1f)
xy = frontier(data1, Q)
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)
xy.all
## # A tibble: 100 × 6
## s s1 s3 er er1 er3
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.0645 0.0661 0.0660 0.00980 0.00980 0.00980
## 2 0.0641 0.0657 0.0655 0.00985 0.00985 0.00985
## 3 0.0637 0.0652 0.0651 0.00990 0.00990 0.00990
## 4 0.0632 0.0648 0.0647 0.00995 0.00995 0.00995
## 5 0.0628 0.0643 0.0642 0.0100 0.0100 0.0100
## 6 0.0624 0.0639 0.0638 0.0101 0.0101 0.0101
## 7 0.0620 0.0635 0.0634 0.0101 0.0101 0.0101
## 8 0.0616 0.0631 0.0630 0.0102 0.0102 0.0102
## 9 0.0612 0.0626 0.0626 0.0102 0.0102 0.0102
## 10 0.0608 0.0622 0.0622 0.0102 0.0102 0.0102
## # ℹ 90 more rows
class(xy.all)
## [1] "tbl_df" "tbl" "data.frame"
plot(xy$sd, xy$er, type = 'l', col="red", xlim = c(0.03, 0.07))
lines(xy.1f$sd, xy.1f$er, col = "blue")
lines(xy.3f$sd, xy.3f$er, col = "black")
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")
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00