In-class Exercise: Factor Models

library(tidyverse)

Data Processing

retdata = read_csv("C:/Users/Admin/Downloads/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'; 

2 approaches to estimate covariance matrix

Method 1: by “lm” function

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]
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

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); 
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

Using FF 3 factor model to compute covariance matrix

Method 1: by “lm” function

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,]
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

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) 
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

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

Create frontier function to plot efficient frontier

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
}

Use different covariance matrix to plot efficient frontier: Q

Q.3f = cov_3f
Q.1f = cov_1f
head(retdata1)

##             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(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) 

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(xx, yy)

Using plot()

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")

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")