Factor Modeling in R

[This article was first published on Posts on Matthew Smith R Shenanigans, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.












The most popular models for analysing the risk of portfolios are factor models, since stocks have a tendency to move together. The principal component of securities often explains a large share of it’s variance. Since we are mostly concerned with multiple assets which form a portfolio we need to account for this. Some questions might be why do stocks with low Price to Book ratios outperform stocks with higher Price to Book ratios? Here the “Price” part of the ratio is simply the share price (per share) and the “Book” part of the ratio is the “Shareholders Equity” / “Shares Outstanding” which are items found on a companies Balance Sheet.

Another important reason Factor models are popular are dimensionality. Assume \(N\) assets which subsequently has \(N\) variances and \(N(N-1)/2\) correlations. If we have a factor model with \(F\) factors, it has \(N\) idiosyncratic variances, \(F\) factor variances and \(NF\) betas. As long as:

\[2F < N(N-1)/(N+1)\]

We have fewer parameters with factor models.

The Capital Asset Pricing Model Sharpe (1964) is the simplest factor model consisting of a single factor.

\[r_i = a_i + b_if + e_i\]

The mean-variance parameters are:

\[\overline{r}_i = a_{i} + b_{i}\overline{y}\]
\[\sigma_{i}^2 = b_{i}^2\sigma_{\epsilon_{i}}^2\]

\[\sigma_{ij} = b_{i}b_{j}\sigma_{y}^2\]

The \(Cov(r_{i}, y) = Cov(a_{i} + b_{i}y + e_{i}, y) = b_{i}\sigma_{y}^2\) and therefore \(b_{i}\) is;

\[b_{i} = \dfrac{Cov(r_{i}, y)}{\sigma_{y}^2}\]

This is alternatively expressed as;

\[b_{i} = Cov(r_{i}, y) / Var(y)\]

Where \(b\) is the slope of the regression line which minimises the squared distance of returns and \(a\) is the intercept or \(alpha\) which are expressed on a line \(a + bf\).

I will go through these calculations using Base R functions, however first we need some data and for this I like the tidy functions which require a few packages.

I download daily price data from Yahoo Finance using quantmod or tidyquant’s wrapper to the quantmod package. The difference here being that quantmod collects data and stores it as an xts object and tidyquant collects the data and stores it as a tibble, from here we are able to more easily use functions from the tidyverse way of handling data, we convert the data back to an xts object using the tk_xts function from the timetk package.

The firms are Apple, Harley Davidson, Nvidia, Microsoft, AMD, Ford, General Electric, 3M, Intel and Amazon.

library(quantmod)
library(ggplot2)
library(tidyquant)
library(tidyr)
library(dplyr)
library(timetk)
library(tibble)
library(rvest)
library(corrplot) 
library(stringr)

start_date <- "2017-01-01"
end_date <- "2019-01-01"

symbols <- c("AAPL", "HOG", "NVDA", "MSFT", "AMD", "F", "GE", "MMM", "INTL", "AMZN")

portfolio_prices <- tq_get(
  symbols,
  from = start_date,
  to = end_date
  ) %>%
  select(symbol, date, adjusted) %>% 
  pivot_wider(names_from = symbol, values_from = adjusted) %>% 
  tk_xts(date_var = date)

The data looks like the following, I removed the Open, High, Low, Close and Volume data and kept just the Adjusted prices, where each asset is it’s own column, the data has been converted into a time series object or an xts object where the date is stored as the index (or rownames) and is not visable in the table below.

AAPL HOG NVDA MSFT AMD F GE MMM INTL AMZN
110.9539 53.59669 101.0116 59.49657 11.43 10.43042 28.13062 165.3643 40.41 753.67
110.8297 54.18777 103.3683 59.23036 11.43 10.91093 28.13950 165.6151 40.86 757.18
111.3933 54.24233 100.7443 59.23036 11.24 10.57954 27.97971 165.0485 38.77 780.45
112.6351 53.74219 102.0910 59.74376 11.32 10.57126 28.05961 165.5315 38.71 795.99
113.6668 52.87832 106.2301 59.55362 11.49 10.46355 27.92646 164.6399 38.33 796.92
113.7815 53.19659 105.4280 59.53460 11.44 10.64582 27.84656 163.9991 38.46 795.90

I also collect the S&P500 data using the same method:

SPY <- tq_get(
  "^GSPC",
  from = start_date,
  to = end_date
) %>% 
  select(date, adjusted) %>% 
  tk_xts(date_var = date)

Which looks similar to the individual asset prices.

adjusted
2257.83
2270.75
2269.00
2276.98
2268.90
2268.90

We can plot the data using the chartSeries function.

chartSeries(SPY)

I compute the daily log returns for all the assets in the portfolio_returns and the SPY returns

portfolio_returns <- diff(log(portfolio_prices), na.pad = FALSE)
SPY_returns <- diff(log(SPY), na.pad = FALSE)

The returns data for each asset looks like:

AAPL HOG NVDA MSFT AMD F GE MMM INTL AMZN
-0.001119731 0.010967925 0.023062871 -0.0044844490 0.000000000 0.0450387901 0.0003155849 0.001515279 0.011074335 0.004646413
0.005072384 0.001006252 -0.025713094 0.0000000000 -0.016762633 -0.0308429974 -0.0056944999 -0.003426713 -0.052504862 0.030269695
0.011086527 -0.009263191 0.013278832 0.0086305184 0.007092228 -0.0007833275 0.0028513212 0.002921873 -0.001548813 0.019715919
0.009117836 -0.016204986 0.039742771 -0.0031877844 0.014906019 -0.0102405428 -0.0047565496 -0.005400948 -0.009865007 0.001167666
0.001008053 0.006000910 -0.007578959 -0.0003193095 -0.004361106 0.0172690614 -0.0028650791 -0.003899774 0.003385783 -0.001280696
0.005358706 -0.008411454 -0.012380417 0.0090613598 -0.021202208 -0.0141067872 0.0031828276 0.007391130 0.009059257 0.003912422

We can use the autoplot to plot xts or time series data with ggplot functionality, (I only plot the first 3 assets):

autoplot(portfolio_returns[,1:3]) +
  xlab("Date") +
  ggtitle("Asset Log Returns") +
  theme_tq()

Technical R:

From here I will talk a little more technically but will keep it as simple as possible.

Next I compute the calculations described at the begining of the post. We first need to define a few things:

  • The number of assets we have in the portfolio, denoted previously as \(N\)
  • The number of days we apply our model which is usually denoted as \(T\).
N_portfolio <- ncol(portfolio_prices)
days <- nrow(portfolio_prices)

(a) We can compute \(\beta\) as the following:

beta <- cov(portfolio_returns, SPY_returns) / as.numeric(var(SPY_returns))

Which returns:

adjusted
AAPL 1.2722563
HOG 0.8957108
NVDA 1.9681667
MSFT 1.4221581
AMD 1.8544079
F 0.8626349

However, in order to understand it better we can break it down each of the computations:

The \(cov\) covariance matrix part of the \(beta\) formula looks like:

adjusted
AAPL 0.0000852
HOG 0.0000600
NVDA 0.0001318
MSFT 0.0000953
AMD 0.0001242
F 0.0000578
GE 0.0000548
MMM 0.0000714
INTL 0.0000576
AMZN 0.0001056

We can compute the covariances using base R as follows:

\[Cov(x, y) = \sum_{i=1}^N \dfrac{((x_i - \overline{x}) (y_{i} - \overline{y}))}{(n-1)}\]
Where \(x_{i}\) is our asset and \(y_{i}\) is the SPY500. In R we can compute this as:

x = portfolio_returns[,1] # which takes the first Asset in the portfolio_returns data
x_bar = colMeans(x)

y = SPY_returns[,1]
y_bar = colMeans(y)

n = length(x)
sum(((x - x_bar)*(y - y_bar))) / (n-1)
## [1] 0.00008521578

Which corresponds to the first result, AAPL in the cov(portfolio_returns, SPY_returns) calculation. We can change the \(1\) in x = portfolio_returns[,1] to \(2\), \(3\), \(4\)\(N\) to compute the covariances for each of the assets (HOG, NVDA and MSFT, respectively), or we can wrap the above into a for loop to compute for all assets. Note: I only changed the \(1\) to an \(i\) from the above equation in the for loop, everything else is constant.

for(i in 1:ncol(portfolio_returns)){
  x = portfolio_returns[,i]
  x_bar = colMeans(x)
  y = SPY_returns[,1]
  y_bar = colMeans(y)
  cov_x_y = sum(((x - x_bar)*(y - y_bar))) / (n-1)
  print(cov_x_y)
}
## [1] 0.00008521578
## [1] 0.00005999475
## [1] 0.0001318279
## [1] 0.00009525621
## [1] 0.0001242083
## [1] 0.00005777932
## [1] 0.00005477266
## [1] 0.00007136762
## [1] 0.0000575575
## [1] 0.0001056044

The variance of the SPY returns is.

x
0.000067

Which is calculated as:

\[Var = \sum(\dfrac{((y - \overline{y})^2)} {(n-1)})\]
In base R we can simply calculate it as:

n = length(SPY_returns)
y = SPY_returns[,1]
y_bar = colMeans(y)

sum(((y - y_bar)**2) / (n-1))
## [1] 0.00006698004

Putting all this together, we can compute \(beta\). Recall, that \(\beta = \dfrac{Cov(r_{i}, y)}{Var(y)}\) where here \(r_{i}\) is each of our assets in our portfolio and \(y\) is the market returns or SPY500 returns.

To compute \(beta\) for each of our assets using base R we can wrap the above code into a function:

myBetaComputation <- function(ASSET){
  x = ASSET                   # which Asset in the portfolio_returns data we want to compute beta for.
  x_bar = colMeans(x)
  
  y = SPY_returns[,1]
  y_bar = colMeans(y)
  
  cov_x_y <- sum(((x - x_bar)*(y - y_bar))) / (n-1)
  
  
  n = length(SPY_returns)
  y = SPY_returns[,1]
  y_bar = colMeans(y)
  
  var_y <- sum(((y - y_bar)**2) / (n-1))
  
  return(cov_x_y / var_y)
}

We can apply this function to individual assets in our data and then all of them:

myBetaComputation(portfolio_returns[,1])        # Computes for the first asset
## [1] 1.272256
myBetaComputation(portfolio_returns[,2])        # Computes for the second asset
## [1] 0.8957108
sapply(portfolio_returns, myBetaComputation)    # Computes for all assets
##      AAPL       HOG      NVDA      MSFT       AMD         F        GE       MMM 
## 1.2722563 0.8957108 1.9681667 1.4221581 1.8544079 0.8626349 0.8177459 1.0655057 
##      INTL      AMZN 
## 0.8593232 1.5766544

The interpretation here is that a value of 1 indicates that the asset moves in perfect correlation with the market, values \(>\) 1 indicate that an asset moves more than when the market moves or moves with more volatility when the market moves and values \(<\) 1 indicates that an asset moves less than the market.

Typically Tech stocks have a higher market \(\beta\) whereas mature non-tech stocks have a lower market \(\beta\). That is, AAPL, NVDA, MSFT and AMD all have betas above 1 whereas HOG and Ford (F) have \(\beta\)’s less than 1.

(b) We can compute \(\alpha\) as the following:

Now that we have \(\beta\) we can estimate \(\alpha\) for each of our assets in our portfolio.

Recall, that alpha is defined as;

\[\alpha = \overline{r_{i}} - \beta_{i}*\overline{f}\]

Where \(\overline{r_{i}}\) is the average returns of our individual assets in our portfolio, \(\beta_{i}\) is the beta of each asset over our sample period and \(\overline{f}\) is the average returns of SPY. Therefore, we can compute \(\alpha\) as;

alpha <- colMeans(portfolio_returns) - beta*colMeans(SPY_returns)

The alpha or \(\alpha\) for our assests, looks like:

adjusted
AAPL 0.0004068
HOG -0.0011516
NVDA 0.0001393
MSFT 0.0007480
AMD 0.0005696
F -0.0009272

There are other ways to calculate \(\beta\) and \(\alpha\) but for a single factor model we can use the CAPM.beta and CAPM.alpha from the PerformanceAnalytics package.

library(PerformanceAnalytics)
CAPM.beta(Ra = portfolio_returns, Rb = SPY_returns, Rf = 0)
##                    AAPL       HOG     NVDA     MSFT      AMD         F
## Beta: adjusted 1.272256 0.8957108 1.968167 1.422158 1.854408 0.8626349
##                       GE      MMM      INTL     AMZN
## Beta: adjusted 0.8177459 1.065506 0.8593232 1.576654
CAPM.alpha(Ra = portfolio_returns, Rb = SPY_returns, Rf = 0)
##                        AAPL          HOG         NVDA         MSFT          AMD
## Alpha: adjusted 0.000406813 -0.001151611 0.0001393021 0.0007480208 0.0005695635
##                             F           GE           MMM          INTL
## Alpha: adjusted -0.0009272497 -0.002875334 0.00001141014 -0.0003782044
##                        AMZN
## Alpha: adjusted 0.001047154

Which is much simpler than what we just did. The Ra is the asset returns, Rb is the market returns and Rf is the risk-free rate of return.

Upon investigating the CAPM.beta functions from the PerformanceAnalytics package I noticed they had functions for CAPM.beta.bull and CAPM.beta.bear, so I wanted to see how these would look like plotted for each asset.

rbind(
  CAPM.beta.bull(Ra = portfolio_returns, Rb = SPY_returns, Rf = 0),
  CAPM.beta(Ra = portfolio_returns, Rb = SPY_returns, Rf = 0),
  CAPM.beta.bear(Ra = portfolio_returns, Rb = SPY_returns, Rf = 0)
) %>% 
  as_tibble(rownames = NA) %>% 
  rownames_to_column("betas") %>% 
  pivot_longer(2:ncol(.)) %>% 
  ggplot(aes(x = name, y = value, color = betas)) +
  geom_point(size = 5, shape = 19) +
  scale_color_brewer(type='seq', palette='Reds') +
  theme_tq() +
  ggtitle("Asset Betas:", subtitle = "Bear, Beta, Bull")

Interesting…

Visualising the covariance matrix of the returns

In order to visualise the the covariance matrix some further formulas are needed.

Ultimately we want to compute the following;

\[\hat{\Sigma} = var(y)\hat{\beta} \hat{\beta}^T + \hat{\psi}\]
Where;
\[\hat{\psi} = diag(\hat{\sigma}_{1}^2, ... \hat{\sigma}_{i}^2, ... \hat{\sigma}_{N}^2)\]

We use our \(\beta\) and \(\alpha\) results from before as well as \(N\). I created a function which takes in the \(asset\) and computes the residules and Sigma values. What we are calculating here is the following;

err
\[\hat{\epsilon_{i}} = x_{i} - \alpha_{i} - \beta_{i}*y\]

Where \(i = 1.,..., N\)

Sigma
\[\hat{\sigma}_{i}^2 = \dfrac{1} {N-2}\hat{\epsilon}_{i}^N\epsilon_{i}\]

The code in base R for the above equations is:

beta <- cov(portfolio_returns, SPY_returns) / as.numeric(var(SPY_returns))
alpha <- colMeans(portfolio_returns) - beta*colMeans(SPY_returns)
N = length(SPY_returns)

mySigmaFunction <- function(ASSET){
  err = portfolio_returns[, ASSET] - alpha[ASSET ,] - beta[ASSET ,] * SPY_returns
  Sigma = (1 / (N - 2)) %*% t(err) %*% err
  return(Sigma)
}

mySigmaFunction("AAPL")                               # For Apple
##              AAPL
## [1,] 0.0001178595
mySigmaFunction("F")                                  # For Ford
##                 F
## [1,] 0.0001599711
sapply(colnames(portfolio_returns), mySigmaFunction)  # For all assets in the portfolio
##          AAPL           HOG          NVDA          MSFT           AMD 
## 0.00011785947 0.00021461371 0.00056144410 0.00006590872 0.00126760753 
##             F            GE           MMM          INTL          AMZN 
## 0.00015997110 0.00035027461 0.00007165307 0.00019184188 0.00017665040

Now that we have the \(\sigma\) values. we want to create a matrix with the Sigma values down the digaonal.

Sigma2 <- sapply(colnames(portfolio_returns), mySigmaFunction)
diag_Sigma2 <- diag(Sigma2)
colnames(diag_Sigma2) = colnames(portfolio_returns)
rownames(diag_Sigma2) = colnames(portfolio_returns)

Which creates our diagonal matrix as follows:

AAPL HOG NVDA MSFT AMD F GE MMM INTL AMZN
AAPL 0.0001179 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.0000000 0.0000000 0.0000000 0.0000000
HOG 0.0000000 0.0002146 0.0000000 0.0000000 0.0000000 0.00000 0.0000000 0.0000000 0.0000000 0.0000000
NVDA 0.0000000 0.0000000 0.0005614 0.0000000 0.0000000 0.00000 0.0000000 0.0000000 0.0000000 0.0000000
MSFT 0.0000000 0.0000000 0.0000000 0.0000659 0.0000000 0.00000 0.0000000 0.0000000 0.0000000 0.0000000
AMD 0.0000000 0.0000000 0.0000000 0.0000000 0.0012676 0.00000 0.0000000 0.0000000 0.0000000 0.0000000
F 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00016 0.0000000 0.0000000 0.0000000 0.0000000
GE 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.0003503 0.0000000 0.0000000 0.0000000
MMM 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.0000000 0.0000717 0.0000000 0.0000000
INTL 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.0000000 0.0000000 0.0001918 0.0000000
AMZN 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.0000000 0.0000000 0.0000000 0.0001767

Now that we have our diagonal matrix \(Diag(\psi)\) we can compute the covariance matrix of the returns using: \[\Sigma = var(y)\beta \beta^T + \psi\]

Then use the cov2cor function from the stats package.

Sigma <- as.numeric(var(SPY_returns)) * beta %*% t(beta) + diag_Sigma2
corrplot(cov2cor(Sigma), main = "Covariance matrix of Portfolio Assets")

The correlation table looks like:

AAPL HOG NVDA MSFT AMD F GE MMM INTL AMZN
AAPL 1.0000000 0.3097538 0.3891498 0.5677312 0.2714307 0.3373737 0.2330694 0.4966757 0.3133839 0.4821634
HOG 0.3097538 1.0000000 0.2515805 0.3670311 0.1754765 0.2181079 0.1506764 0.3210946 0.2025988 0.3117126
NVDA 0.3891498 0.2515805 1.0000000 0.4611084 0.2204547 0.2740132 0.1892978 0.4033975 0.2545288 0.3916107
MSFT 0.5677312 0.3670311 0.4611084 1.0000000 0.3216216 0.3997582 0.2761668 0.5885171 0.3713324 0.5713213
AMD 0.2714307 0.1754765 0.2204547 0.3216216 1.0000000 0.1911233 0.1320346 0.2813684 0.1775330 0.2731471
F 0.3373737 0.2181079 0.2740132 0.3997582 0.1911233 1.0000000 0.1641118 0.3497257 0.2206640 0.3395071
GE 0.2330694 0.1506764 0.1892978 0.2761668 0.1320346 0.1641118 1.0000000 0.2416026 0.1524423 0.2345432
MMM 0.4966757 0.3210946 0.4033975 0.5885171 0.2813684 0.3497257 0.2416026 1.0000000 0.3248576 0.4998164
INTL 0.3133839 0.2025988 0.2545288 0.3713324 0.1775330 0.2206640 0.1524423 0.3248576 1.0000000 0.3153656
AMZN 0.4821634 0.3117126 0.3916107 0.5713213 0.2731471 0.3395071 0.2345432 0.4998164 0.3153656 1.0000000

ETF analysis and a randomly sampled portfolio

Analysing my randomly selected portfolio and a series of U.S. funds

Since factor modelling is all about risk and portfolio analysis I thought it would be interesting to compare the performance of a number of U.S. Exchnage Traded Funds (ETFs) and a randomly selected portfolio of assets from the SP500.

In order to construct the randomly created portfolio I first scraped the wikipedia page which contains the list of SP500 firms along with their ticker symbol, I filter out all Class A, B and C shares since a few firms have multiple asset classes and I don’t want to sample two of the same asset.

Just as before I collect the data and put it into time series format.

start_date <- "2016-01-01"
end_date <- "2017-12-31"

set.seed(4746)
url <- "https://en.wikipedia.org/wiki/List_of_S%26P_500_companies"

symbols <- url %>%
  read_html() %>%
  html_nodes(xpath = '//*[@id="constituents"]') %>% 
  html_table() %>% 
  .[[1]] %>% 
  filter(!str_detect(Security, "Class A|Class B|Class C")) %>%     # Removes firms with Class A shares
  sample_n(20) %>% 
  pull(Symbol)

asset_prices <- tq_get(
  symbols,
  from = start_date,
  to = end_date
) %>%
  select(symbol, date, adjusted) %>% 
  pivot_wider(names_from = symbol, values_from = adjusted) %>% 
  tk_xts(date_var = date)

Which looks like:

MMM MNST TJX TMUS NWL INTC EIX XRX GPN HES TIF IR GRMN DRI GLW ZTS CCI BK
132.8015 48.11333 33.27088 38.95 38.54331 30.45350 51.77293 24.64223 62.47044 44.85172 68.47929 51.15904 30.85206 55.28202 16.24429 45.97239 73.96610 36.90321
133.3804 48.69000 33.79228 40.22 38.52557 30.31014 51.79048 24.59438 63.43876 44.74931 68.27712 50.61163 30.96944 56.28635 16.23521 46.69207 75.48198 36.76471
130.6940 48.76667 32.92327 40.05 37.46084 29.63818 51.44798 24.28336 63.63840 41.95655 67.00899 49.53537 30.12268 55.86901 15.87241 46.70180 74.84963 35.84144
127.5101 48.40000 32.87630 40.51 35.78388 28.52719 51.19327 23.61348 61.25259 40.60672 65.17110 48.06017 29.25916 55.47833 15.50961 45.28187 72.98727 34.93664
127.0759 48.08333 31.43894 39.88 34.94097 28.23153 51.18450 22.94359 60.07465 40.56948 62.72673 47.60556 28.12736 55.06099 15.55497 44.62054 71.96515 34.29958
127.0488 48.68000 31.84760 39.68 34.91435 28.72431 51.99249 22.48903 58.41755 38.81003 62.34076 47.76327 27.83393 55.01660 15.70009 43.35622 72.49354 34.54886

Next, I collect the U.S. ETFs. I scraped a website (which I cannot remember the name) for the symbols of the ETFs.

myETFs <- c("SPY", "IVV", "VTI", "VOO", "QQQ", "VEA", "IEFA", "AGG", "VWO",
            "EFA","IEMG","VTV", "IJH", "IWF","BND", "IJR", "IWM", "VUG",
            "GLD", "IWD", "VIG", "VNQ", "USMV", "LQD", "VO", "VYM", "EEM",
            "VB", "VCSH", "XLF", "VCIT", "VEU", "XLK", "ITOT", "IVW", "BNDX",
            "VGT", "DIA", "BSV", "SHV", "IWB", "IWR", "TIP", "SCHF", "MBB", "SDY",
            "MDY", "SCHX", "IEF", "HYG", "DVY", "XLV", "SHY", "IXUS", "TLT", "IVE",
            "PFF", "IAU", "VXUS", "RSP", "SCHB", "VV", "GOVT", "EMB", "MUB", "QUAL",
            "XLY", "VBR", "EWJ", "XLP", "VGK", "SPLV", "MINT", "BIV", "IGSB", "EFAV",
            "VT", "GDX", "XLU", "IWS", "XLI", "SCHD", "IWP", "ACWI", "VMBS", "XLE", "JNK",
            "VOE", "FLOT", "IWV", "JPST", "SCZ", "IEI", "IWN", "DGRO", "VBK", "IGIB", "IWO")
            

etf_prices <- tq_get(
  myETFs,
  from = start_date,
  to = end_date
) %>% 
  select(symbol, date, adjusted) %>% 
  pivot_wider(names_from = symbol, values_from = adjusted) %>% 
  tk_xts(date_var = date)

The data looks like:

SPY IVV VTI VOO QQQ VEA IEFA AGG VWO EFA IEMG VTV IJH IWF BND IJR IWM VUG GLD IWD VIG VNQ USMV LQD VO VYM EEM VB VCSH XLF VCIT VEU XLK ITOT IVW BNDX VGT DIA BSV SHV IWB IWR TIP SCHF MBB SDY MDY SCHX IEF HYG DVY XLV SHY IXUS TLT IVE PFF IAU VXUS RSP SCHB VV GOVT EMB MUB QUAL XLY VBR EWJ XLP VGK SPLV MINT BIV IGSB EFAV VT GDX XLU IWS XLI SCHD IWP ACWI VMBS XLE JNK VOE FLOT IWV JPST SCZ IEI IWN DGRO VBK IGIB IWO
186.8362 187.3943 95.69553 171.0293 105.7754 32.31758 48.12894 97.70714 28.82161 51.85117 35.29821 73.21329 129.8436 93.32362 72.46418 51.13569 104.60918 99.84885 102.89 88.57166 71.16028 67.15640 38.30936 100.1451 112.0916 58.70417 29.25530 103.10962 71.86420 13.61419 73.99710 38.23042 39.89517 42.61963 107.9444 48.90886 101.77430 156.9451 74.18727 105.2681 104.11566 37.17856 101.8868 24.86581 97.72220 63.99211 238.7614 44.54472 98.23306 65.49518 65.93277 66.68895 80.17432 44.14714 110.3894 80.04479 31.46391 10.38 39.78462 71.26929 45.05346 85.64318 23.57018 87.65514 101.0075 59.58526 72.85894 90.46735 45.07748 45.11800 43.21453 34.96312 93.40020 75.13937 48.21395 57.09138 51.90483 13.86359 38.17657 62.81808 48.48287 34.12284 87.50332 50.81642 48.00686 53.82844 80.80476 78.58987 46.91737 110.6951 NA 45.01355 115.1113 84.09489 23.20199 115.5700 47.86354 131.6090
187.1522 187.7839 95.90974 171.3355 105.5919 32.28184 48.04810 97.75242 28.88509 51.77047 35.38096 73.44967 129.9475 93.46677 72.55397 51.29261 104.83739 99.86796 103.18 88.74566 71.39217 68.45230 38.56908 100.1978 112.2051 58.99791 29.32071 103.16643 71.84602 13.66662 74.01466 38.14082 39.79135 42.69417 108.1532 48.88116 101.28743 157.0367 74.17794 105.2777 104.32082 37.26314 101.8220 24.83831 97.79476 64.25666 238.8375 44.63781 98.20529 65.61782 66.33069 67.00954 80.12685 44.10166 109.9440 80.26437 31.49616 10.40 39.71283 71.43853 45.09076 85.81048 23.54203 87.88744 101.1264 59.74447 72.76409 90.67188 45.64425 45.40763 42.95955 35.24865 93.40020 75.22964 48.19091 57.19847 51.91398 13.79471 38.45059 62.96579 48.61266 34.31977 87.66763 50.86277 48.09793 54.03376 80.97237 78.76543 46.90808 110.9662 NA 44.89528 115.1582 84.47743 23.28449 115.6185 47.92159 131.6866
184.7914 185.3255 94.62438 169.1641 104.5776 31.69214 47.20373 98.12343 28.36815 50.91853 34.74654 72.40415 128.1251 92.48373 72.87716 50.70773 103.23986 98.76987 104.67 87.39035 70.52956 68.26473 38.32790 100.5933 110.1342 58.22349 28.76008 101.62282 71.87331 13.45691 74.25190 37.48678 39.30056 42.08852 106.9669 49.06577 99.91275 154.7919 74.26179 105.2586 102.94081 36.64525 102.1459 24.37987 98.11230 63.51594 235.4678 44.06064 98.79820 65.61782 65.65863 66.46268 80.17432 43.33769 111.4255 79.06572 31.42361 10.57 39.03973 70.13145 44.53121 84.68593 23.63587 87.87912 101.5744 59.03268 72.05276 89.25890 44.85076 45.25377 42.24736 34.93548 93.40947 75.53665 48.22778 56.64514 51.14435 14.02101 38.37987 61.87635 47.86177 33.90801 86.25658 50.09352 48.20723 51.95383 80.90055 77.24093 46.90808 109.4701 NA 44.14021 115.4954 83.48843 22.96364 113.8152 47.99304 129.3380
180.3579 180.8912 92.30511 165.0441 101.3029 31.09350 46.30547 98.11440 27.45217 49.86035 33.69835 70.71315 124.8108 90.12631 72.88611 49.40485 100.48224 96.22995 106.15 85.32989 69.23102 66.93474 37.70643 100.6724 107.6568 56.96841 27.87242 98.88598 71.90975 13.07825 74.29585 36.68042 38.13966 41.07753 104.2529 48.99192 96.69561 151.1452 74.37366 105.2681 100.43254 35.80653 102.0626 23.93059 98.23930 62.23734 229.4516 43.00870 99.01125 65.19263 64.49137 65.11437 80.20277 42.43728 111.6255 77.26316 31.22207 10.72 38.23201 68.54225 43.44009 82.62268 23.68278 87.46435 101.7207 57.62780 70.57320 87.00002 44.17064 44.71072 41.46484 34.35522 93.40020 75.63598 48.24622 55.87762 49.95324 14.64089 38.12353 60.51917 46.56396 33.20084 84.22703 48.92575 48.27098 50.68622 80.46945 75.47622 46.91737 106.8519 NA 43.46700 115.7578 81.26778 22.53279 110.6160 48.06894 125.7083
178.3782 178.8873 91.25257 163.2903 100.4722 30.70930 45.71261 98.33157 27.20731 49.25951 33.37653 69.87674 123.1961 89.22913 72.97586 48.58223 98.75159 95.33240 105.68 84.23099 68.51684 66.03955 37.32612 100.7515 106.2667 56.32751 27.57342 97.51283 71.91882 12.87435 74.32222 36.24140 37.83765 40.59300 103.2470 48.98268 95.86506 149.5692 74.45753 105.2777 99.38823 35.36956 102.1459 23.63719 98.34815 61.63771 226.3864 42.52462 99.26135 65.02093 64.04923 64.13377 80.25970 41.90978 112.1254 76.28411 31.31075 10.66 37.78327 67.70534 42.94580 81.73975 23.71093 87.38972 101.7664 56.97219 69.81445 85.77296 43.22602 44.36678 41.10435 34.02365 93.39094 75.83461 48.27848 55.38676 49.36685 14.28668 38.10586 59.79903 46.09118 32.84278 83.17358 48.36041 48.30741 50.03457 80.18208 74.31207 46.88947 105.6269 NA 42.88478 115.9920 79.95219 22.25777 108.9970 48.12699 123.5344
178.5548 179.0729 91.19671 163.3088 100.7813 30.84332 45.87430 98.05104 27.08941 49.46577 33.33975 69.83128 122.7901 89.27684 72.67066 48.58699 98.32368 95.42790 104.74 84.19437 68.66522 66.41467 37.38177 100.4263 105.7561 56.38982 27.56408 97.03933 71.98254 12.89183 74.19040 36.31308 38.08304 40.58370 103.5317 48.92731 96.33284 150.1098 74.44821 105.2681 99.32295 35.24034 101.6092 23.71054 98.25743 61.62007 225.7010 42.51531 98.94638 64.92279 64.19955 63.37006 80.28816 41.99163 110.8984 76.18345 31.14146 10.57 37.86405 67.48906 42.91782 81.75836 23.67340 87.24041 101.6019 57.05649 70.36455 85.55912 43.52830 44.78313 41.22744 34.09733 93.41876 75.65400 48.26467 55.53847 49.39434 13.69632 38.30916 59.57745 46.11899 32.90544 82.78702 48.41602 48.23457 48.96336 80.06235 74.02564 46.90808 105.5428 NA 42.95756 115.8515 79.92420 22.32194 108.0954 48.04214 122.4960

I also collect the SPY500.

SPY_prices <- tq_get(
  "^GSPC",
  from = start_date,
  to = end_date
) %>% 
  select(date, adjusted) %>% 
  tk_xts(date_var = date)

Which looks like:

adjusted
2012.66
2016.71
1990.26
1943.09
1922.03
1923.67

Now, we have a series of 3 data sets, the daily prices of the adjusted close for the randomly selected assets from the SPY500 wikipedia page, the ETFs and the SPY500. Next I calculate the daily returns.

asset_returns <- diff(log(asset_prices), na.pad = FALSE)
etf_returns <- diff(log(etf_prices), na.pad = FALSE)
SPY_returns <- diff(log(SPY_prices), na.pad = FALSE)

The asset_returns looks like:

MMM MNST TJX TMUS NWL INTC EIX XRX GPN HES TIF IR GRMN DRI GLW ZTS CCI BK
0.0043496940 0.011914276 0.015549833 0.032085611 -0.0004605231 -0.004718686 0.0003390193 -0.001943595 0.015381467 -0.0022860685 -0.002956558 -0.010757820 0.003797292 0.0180042454 -0.0005589373 0.0155331932 0.020287018 -0.003760040
-0.0203468502 0.001573356 -0.026052579 -0.004235761 -0.0280260856 -0.022418948 -0.0066352645 -0.012726538 0.003142188 -0.0644414487 -0.018747883 -0.021494331 -0.027722605 -0.0074420504 -0.0225996450 0.0002084719 -0.008412745 -0.025433678
-0.0246630740 -0.007547151 -0.001427791 0.011420159 -0.0457985304 -0.038205692 -0.0049630049 -0.027973817 -0.038210941 -0.0327010158 -0.027810684 -0.030233211 -0.029085489 -0.0070174571 -0.0231226086 -0.0308759534 -0.025196097 -0.025568744
-0.0034106468 -0.006564266 -0.044704812 -0.015673838 -0.0238374478 -0.010418300 -0.0001713849 -0.028779059 -0.019418177 -0.0009174858 -0.038228556 -0.009504269 -0.039450060 -0.0075510143 0.0029199836 -0.0147124623 -0.014103137 -0.018403085
-0.0002136353 0.012332679 0.012914805 -0.005027688 -0.0007621181 0.017304330 0.0156624749 -0.020011095 -0.027971612 -0.0443372159 -0.006172111 0.003307394 -0.010486917 -0.0008065403 0.0092862449 -0.0287440451 0.007315467 0.007241532
0.0028436862 0.008590726 0.016384277 0.009531073 0.0156331210 0.019154091 0.0033726404 0.013734922 0.000000000 -0.0514319659 0.008367407 0.016185405 -0.001507319 0.0088375907 0.0017316113 0.0087102953 0.004411381 0.007985230

However, I want to assume that I own all of these assets in a single portfolio. So I average over the rows and join the data to the ETFs and call it all_returns.

portfolio_returns <- xts(rowMeans(asset_returns), index(asset_returns))
colnames(portfolio_returns) <- c("myPortfolio")
all_returns <- cbind(portfolio_returns, etf_returns)

The all_returns data looks like the following, where we can now see myPortfolio has been joined with the ETF dataset.

myPortfolio SPY IVV VTI VOO QQQ VEA IEFA AGG VWO EFA IEMG VTV IJH IWF BND IJR IWM VUG GLD IWD VIG VNQ USMV LQD VO VYM EEM VB VCSH XLF VCIT VEU XLK ITOT IVW BNDX VGT DIA BSV SHV IWB IWR TIP SCHF MBB SDY MDY SCHX IEF HYG DVY XLV SHY IXUS TLT IVE PFF IAU VXUS RSP SCHB VV GOVT EMB MUB QUAL XLY VBR EWJ XLP VGK SPLV MINT BIV IGSB EFAV VT GDX XLU IWS XLI SCHD IWP ACWI VMBS XLE JNK VOE FLOT IWV JPST SCZ IEI IWN DGRO VBK IGIB IWO
0.006099968 0.0016898660 0.002077056 0.0022360355 0.0017890863 -0.001736693 -0.001106388 -0.001681088 0.00046332859 0.002200231 -0.001557725 0.0023415111 0.0032234753 0.0007997348 0.0015327343 0.0012383279 0.003064097 0.002179106 0.000191391 0.002814589 0.0019626291 0.003253374 0.019113014 0.006756640 0.0005267369 0.001011874 0.004991256 0.0022334070 0.0005508056 -0.0002531205 0.003843879 0.0002372109 -0.002346225 -0.002605712 0.0017475020 0.001932397 -0.0005666427 -0.004795300 0.0005836581 -0.0001257842 0.00009081159 0.0019685044 0.002272383 -0.0006361631 -0.001106790 0.0007422784 0.0041255033 0.0003189161 0.0020875396 -0.0002827452 0.001870799 0.0060170694 0.004795681 -0.0005922602 -0.001030700 -0.004042725 0.002739571 0.0010243955 0.001924928 -0.001806146 0.002371882 0.0008275628 0.0019515134 -0.001194977 0.00264674400 0.0011761117 0.002668490 -0.001302665 0.002258164 0.012494986 0.006398807 -0.005917852 0.008133639 0.00000000000 0.001200660 -0.0004779842 0.001874008 0.0001763650 -0.00498036 0.0071519848 0.0023486170 0.0026733893 0.005754619 0.001875951 0.0009116517 0.0018953483 0.003807047 0.0020722096 0.002231295 -0.0001980912 0.002446725 NA -0.002630911 0.0004069319 0.0045386533 0.003549552 0.0004194414 0.001211963 0.0005894971
-0.017322911 -0.0126945384 -0.013178202 -0.0134924019 -0.0127544774 -0.009652265 -0.018436329 -0.017729720 0.00378821023 -0.018058472 -0.016592869 -0.0180938802 -0.0143367503 -0.0141237737 -0.0105732344 0.0044445716 -0.011468264 -0.015355426 -0.011056244 0.014337489 -0.0153896498 -0.012156174 -0.002743903 -0.006272724 0.0039388941 -0.018628687 -0.013212984 -0.0193056059 -0.0150753935 0.0003798793 -0.015464153 0.0032002389 -0.017296894 -0.012410734 -0.0142874994 -0.011028538 0.0037696378 -0.013665012 -0.0143976824 0.0011297112 -0.00018126091 -0.0133167477 -0.016720820 0.0031757481 -0.018629431 0.0032416931 -0.0115945095 -0.0142095765 -0.0130144362 0.0060192314 0.000000000 -0.0101836091 -0.008194411 0.0005922602 -0.017474795 0.013385186 -0.015046456 -0.0023059217 0.016213994 -0.017094464 -0.018466064 -0.0124870774 -0.0131916717 0.003978224 -0.00009471654 0.0044202968 -0.011985546 -0.009823880 -0.015706040 -0.017537150 -0.003394215 -0.016716991 -0.008924441 0.00009924540 0.004072733 0.0007647688 -0.009721009 -0.0149360906 0.01627195 -0.0018408588 -0.0174535242 -0.0155668847 -0.012070245 -0.016226286 -0.0152395434 0.0022698272 -0.039253655 -0.0008873999 -0.019544672 0.0000000000 -0.013574412 NA -0.016961509 0.0029245001 -0.0117764103 -0.013875422 -0.0157191560 0.001489930 -0.0179955058
-0.022721172 -0.0242840940 -0.024218220 -0.0248157103 -0.0246568769 -0.031814151 -0.019069875 -0.019212776 -0.00009198023 -0.032821816 -0.021000906 -0.0306311536 -0.0236320930 -0.0262082119 -0.0258207050 0.0001227881 -0.026029826 -0.027074043 -0.026052010 0.014040682 -0.0238600026 -0.018582970 -0.019675179 -0.016347691 0.0007861946 -0.022751878 -0.021792064 -0.0313506787 -0.0273006374 0.0005068469 -0.028542194 0.0005917693 -0.021745295 -0.029984083 -0.0243136677 -0.025700013 -0.0015062361 -0.032729302 -0.0238405527 0.0015052940 0.00009044932 -0.0246679043 -0.023153370 -0.0008155098 -0.018599946 0.0012935879 -0.0203358057 -0.0258820165 -0.0241643775 0.0021541043 -0.006500818 -0.0179376149 -0.020495249 0.0003548886 -0.020995449 0.001793079 -0.023062135 -0.0064343040 0.014091356 -0.020906723 -0.022920923 -0.0248076812 -0.0246652173 0.001982686 -0.00473089487 0.0014396023 -0.024085978 -0.020748172 -0.025632797 -0.015280173 -0.012072667 -0.018696026 -0.016748990 -0.00009924540 0.001314087 0.0003822792 -0.013642270 -0.0235646381 0.04326116 -0.0067016093 -0.0221778810 -0.0274902986 -0.021075943 -0.023810522 -0.0235878969 0.0013215006 -0.024701342 -0.0053430008 -0.023111858 0.0001980912 -0.024207994 NA -0.015369037 0.0022689596 -0.0269583342 -0.018940649 -0.0285117041 0.001580334 -0.0284649716
-0.016273779 -0.0110371371 -0.011139358 -0.0114683030 -0.0106831096 -0.008234532 -0.012433352 -0.012885863 0.00221091926 -0.008959566 -0.012123531 -0.0095957668 -0.0118987200 -0.0130211706 -0.0100044748 0.0012306568 -0.016790854 -0.017373511 -0.009370855 -0.004437547 -0.0129618969 -0.010369474 -0.013464265 -0.010137208 0.0007854976 -0.012995981 -0.011313767 -0.0107854401 -0.0139835513 0.0001261363 -0.015713048 0.0003548290 -0.012040822 -0.007949941 -0.0118655695 -0.009695182 -0.0001885182 -0.008626436 -0.0104818256 0.0011270489 0.00009081159 -0.0104525822 -0.012278941 0.0008155098 -0.012336194 0.0011074564 -0.0096812509 -0.0134486432 -0.0113193358 0.0025227807 -0.002637331 -0.0068795036 -0.015174172 0.0007095865 -0.012508199 0.004468584 -0.012752597 0.0028360806 -0.005612737 -0.011806658 -0.012285245 -0.0114437766 -0.0107437626 0.001188048 -0.00085371766 0.0004493159 -0.011441888 -0.010809474 -0.014204613 -0.021617824 -0.007722146 -0.008731859 -0.009698102 -0.00009916959 0.002622702 0.0006684507 -0.008823294 -0.0118083221 -0.02449108 -0.0004634958 -0.0119706508 -0.0102051133 -0.010843211 -0.012586070 -0.0116223810 0.0007544131 -0.012939813 -0.0035775482 -0.015544198 -0.0005950313 -0.011530313 NA -0.013485184 0.0020213357 -0.0163207962 -0.012280109 -0.0147442744 0.001206870 -0.0174444534
-0.003287114 0.0009894013 0.001036607 -0.0006123565 0.0001135276 0.003072098 0.004354785 0.003530857 -0.00285692501 -0.004342847 0.004178409 -0.0011027594 -0.0006508431 -0.0033014245 0.0005345031 -0.0041909610 0.000097994 -0.004342530 0.001001152 -0.008934590 -0.0004349583 0.002163316 0.005664191 0.001489776 -0.0032326845 -0.004816595 0.001105438 -0.0003388982 -0.0048675476 0.0008856206 0.001356817 -0.0017752166 0.001975840 0.006464163 -0.0002292283 0.002753215 -0.0011312023 0.004867711 0.0036078899 -0.0001251127 -0.00009081159 -0.0006570038 -0.003659859 -0.0052679258 0.003098229 -0.0009228731 -0.0002861645 -0.0030322035 -0.0002189089 -0.0031781230 -0.001510439 0.0023442262 -0.011979547 0.0003544986 0.001951243 -0.011003744 -0.001320438 -0.0054213750 -0.008478619 0.002135701 -0.003199632 -0.0006516613 0.0002276234 -0.001584279 -0.00170998028 -0.0016183255 0.001478523 0.007848505 -0.002496171 0.006968673 0.009340380 0.002990074 0.002163358 0.00029781110 -0.002384484 -0.0002860690 0.002735430 0.0005566762 -0.04220025 0.0053208513 -0.0037123602 0.0006030788 0.001905938 -0.004658536 0.0011492677 -0.0015090848 -0.021642045 -0.0014944173 -0.003861882 0.0003969400 -0.000796942 NA 0.001695691 -0.0012123094 -0.0003501705 0.002878845 -0.0083062241 -0.001764683 -0.0084416475
0.005140787 0.0080360599 0.007689359 0.0068197139 0.0082617000 0.011531038 0.004047425 0.004688442 0.00221249905 0.003675782 0.004702448 0.0008273783 0.0061003752 0.0038377668 0.0097876198 0.0025911534 0.002638930 0.002993541 0.009163331 -0.005072985 0.0057480147 0.008340049 -0.006697424 0.005197513 0.0022727943 0.004549844 0.005352840 0.0020317774 0.0028262158 0.0000000000 0.007652486 0.0005922597 0.004185365 0.011825830 0.0070922542 0.010122713 0.0013197206 0.011331708 0.0068127609 0.0001251127 0.00009081159 0.0077616649 0.004124727 0.0027283407 0.004244698 0.0015682489 0.0052809133 0.0041247242 0.0078515507 0.0044839921 0.001007182 0.0009635456 0.012714693 0.0001185908 0.003891078 0.014321948 0.006106584 -0.0007773691 -0.004741593 0.002367465 0.004726248 0.0064978364 0.0075872456 0.003165548 -0.00467027771 0.0003596236 0.010125999 0.010992048 0.001411921 -0.005221970 0.005441886 0.008282949 0.005119095 -0.00019864151 0.001431193 -0.0003823817 0.003208727 0.0055494219 -0.02252194 -0.0039301169 0.0004647994 0.0080082458 0.005967129 0.008138341 0.0064871928 0.0001887472 0.002367193 0.0023899644 0.002866649 0.0001980912 0.007238678 NA 0.000000000 0.0018584037 -0.0026886552 0.008587052 0.0051881359 0.001579281 0.0088344163

Next, I compute (just as before) the \(\beta\) and \(\alpha\) of the portfolios. This time using just the CAPM.beta and CAPM.alpha functions in the PerformanceAnalytics package.

library(PerformanceAnalytics)
betas <- CAPM.beta(Ra = all_returns, Rb = SPY_returns, Rf = 0)
alphas <- CAPM.alpha(Ra = all_returns, Rb = SPY_returns, Rf = 0)

beta_alphas <- cbind(t(betas), t(alphas))
colnames(beta_alphas) <- c("beta", "alpha")

I rank the ETF’s by their \(\alpha\) value. We can see the GDX has the highest \(\alpha\) and a negative \(\beta\) which makes sense since its a Goldminers ETF whihc tracks the Arca Gold Miners Index (GDMNTR) which tracks the performance of companies involved in the gold mining industry. As far as I know there are very few Gold mining companies listed in the SPY500. However, this is a nice way to present and rank ETFs by their \(\alpha\) values and see their corresponding \(\beta\) values.

beta_alphas %>% 
  as_tibble(rownames = NA) %>% 
  rownames_to_column("ETFs") %>% 
  arrange(desc(alpha)) %>% 
  head(10) %>% # I only take the top 10 observations since there are many ETF's in the data
  pivot_longer(cols = c("beta", "alpha")) %>% 
  ggplot(aes(x = reorder(ETFs, if_else(name == "alpha", value, -Inf), FUN = max), y = value)) +
  geom_bar(color = "black", fill = "darkred", stat = "identity") +
  scale_color_brewer(type='seq', palette='Reds') + 
  facet_wrap(~name, scales = "free") +
  theme_tq() +
  coord_flip() +
  ggtitle("Highest ranked  alpha stocks", subtitle = "With corresponding beta values") +
  ylab("Alpha / Beta value") +
  xlab("ETF")

We can also rank the \(alpha\) and \(beta\) values by taking their ratio of \(\dfrac{\alpha}{\beta}\) and plot the results.

beta_alphas %>% 
  as_tibble(rownames = NA) %>% 
  rownames_to_column("ETFs") %>% 
  mutate(alpha_divided_beta = alpha/beta) %>% 
  arrange(desc(alpha_divided_beta )) %>% 
  head(10) %>% 
  ggplot(aes(x = ETFs, y = alpha_divided_beta )) +
  geom_bar(color = "black", fill = "darkred", stat = "identity") +
  scale_color_brewer(type='seq', palette='Reds') + 
  theme_tq() +
  coord_flip() +
  ggtitle("Best Performers: Alpha over Beta ranked ETFs") +
  ylab("Alpha / Beta value") +
  xlab("ETF")

Changing a few lines of code and we can plot and see the worst performers with negative \(\alpha\) to \(\beta\) ratios.

beta_alphas %>% 
  as_tibble(rownames = NA) %>% 
  rownames_to_column("ETFs") %>% 
  mutate(alpha_divided_beta = alpha/beta) %>% 
  arrange(alpha_divided_beta) %>% 
  head(10) %>% 
  ggplot(aes(x = ETFs, y = alpha_divided_beta )) +
  geom_bar(color = "black", fill = "darkblue", stat = "identity") +
  scale_color_brewer(type='seq', palette='Reds') + 
  theme_tq() +
  coord_flip() +
  ggtitle("Worst Performers: Alpha over Beta ranked ETFs") +
  ylab("Alpha / Beta value") +
  xlab("ETF")

We can take our ETF \(\beta\) and \(\alpha\) values before and compute \(\sigma\).

beta <- cov(all_returns, SPY_returns) / as.numeric(var(SPY_returns))
alpha <- colMeans(all_returns) - beta*colMeans(SPY_returns)
N = length(SPY_returns)

mySigmaFunction <- function(ASSET){
  err = all_returns[, ASSET] - alpha[ASSET ,] - beta[ASSET ,] * SPY_returns
  Sigma = (1 / (N - 2)) %*% t(err) %*% err
  return(Sigma)
}

Sigma2 <- sapply(colnames(all_returns), mySigmaFunction)
diag_Sigma2 <- diag(Sigma2)
colnames(diag_Sigma2) = colnames(all_returns)
rownames(diag_Sigma2) = colnames(all_returns)

Sigma <- as.numeric(var(SPY_returns)) * beta %*% t(beta) + diag_Sigma2

Finally, I take a random sample of the ETF’s (since there are too many to analyse) and plot the correlations between the ETF’s.

set.seed(1234)
Sigma_subsample <- sample(colnames(Sigma), size = 10)
Sigma_sampled <- Sigma[ncol = Sigma_subsample, nrow = Sigma_subsample, drop = FALSE]
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(cov2cor(Sigma_sampled), method = "color", col = col(200),
         main = "Covariance matrix of ETFs",
         type = "upper", order = "hclust", number.cex = .7,
         addCoef.col = "black",
         tl.col = "black", tl.srt = 90,
         diag = FALSE)

Sharpe Ratio, CAPM and Fama-French Factor Analysis

Using simple plots still doesn’t give us enough information about an ETF, Portfolio or Asset. The Sharpe ratio Sharpe (1966) is a better measure. The Sharpe ratio is a reward-to-variability ratio which allows us to compare the performance of a portfolio compared to a risk-free asset, after adjusting for it’s risk. It takes the differenece between the portfolio’s return and the risk-free return, this is then divided by the standard deviation (which measures the portfolios volatility).

The Sharpe ratio tells us what the additional unit of return we can expect per unit of risk increase. The Sharpe ratio is defined as:

\[S = \dfrac{\overline{d}}{\sigma_{d}}\]

Where \[\widetilde{d} = \widetilde{R}_{F} - \widetilde{R}_{B}\]

The single asset model with just the market factor looks like:

\[ x_t = \alpha + \beta y_{t} + \epsilon_{t}\]

Taking the expected value of \(x\) and \(y\) at time \(t\) we have:

\[E[x_{t}] = \alpha + \beta E[y_{t}]\]
With variance:

\[var[x_{t}] = \beta^2var[y_t] + \sigma^2\]

Thus, the Sharpe ratio becomes:

\[\dfrac{E[x_{t}]}{\sqrt{var[x_{t}]}}\]

Recall that \(\sigma = \sqrt{var[x_t]}\), we just plug in the values for \(E[x_{t}]\) and \(var[x_t]\) defined previously and the equation becomes:

\[\dfrac{\alpha + \beta E[y_{t}]}{\sqrt{\beta^2 var[y_t] + \sigma^2}}\]
Which after some algebra we obtain:

\[\dfrac{\dfrac{\alpha}{\beta} + E[y_{t}]}{\sqrt{var[y_{t}] + \dfrac{\sigma^{2}}{\beta^{2}}}}\]

Finally we end up with:

\[\dfrac{\dfrac{\alpha}{\beta} + E[y_t]}{\sqrt{(var[y_t])}}\]
We can represent the Sharpe ratio as \(\dfrac{\overline{x}}{\sqrt(Diag(\sigma^{2}))}\) where \(\overline{x}\) is the average value of \(x\) over the historic period from \(t=1\) to \(T\) - simply computed as \(\overline{x} = \dfrac{1}{T}\sum_{t=1}^{T} D_{t}\). In R we can simply take colMeans(all_returns) / sqrt(diag(var(all_returns))).

Machine Learning and Clustering

We can cluster our ETF’s based on their \(\beta\), \(\alpha\) and \(sharpe\) values. Why would we want to do this? Perhaps with more variables than the 3 here the model can cluster these firms in a higher dimensional space and thus we can select our ETF’s based on the clusters and use it as a portfolio diversification tool, that is one cluster might contain riskier ETF’s whereas another might contain value stocks or growth stocks. With such few variables here it doesn’t quite work but with more variables we can better classifiy ETFs.

betas <- CAPM.beta(Ra = all_returns, Rb = SPY_returns, Rf = 0)
alphas <- CAPM.alpha(Ra = all_returns, Rb = SPY_returns, Rf = 0)
sharpe <- colMeans(all_returns) / sqrt(diag(var(all_returns)))

beta_alphas_sharpe <- cbind(t(betas), t(alphas), sharpe)
colnames(beta_alphas_sharpe) <- c("beta", "alpha", "sharpe")
beta_alphas_sharpe <- na.omit(beta_alphas_sharpe)
library(factoextra)
kmeans_model <- kmeans(beta_alphas_sharpe, centers = 5, iter.max = 10, nstart = 1)
fviz_cluster(kmeans_model, data = beta_alphas_sharpe)

The beta_alphas_sharpe data looks like:

beta alpha sharpe
myPortfolio 1.0977199 0.0000345 0.0834437
SPY 0.9990763 0.0000799 0.0988232
IVV 0.9923880 0.0000872 0.1000217
VTI 1.0282825 0.0000718 0.0969068
VOO 0.9978598 0.0000848 0.0995858
QQQ 1.1472461 0.0000931 0.0870931

We can plot the \(\beta\), \(\alpha\), and \(sharpe\) values in a 3D plot and colour them according to the clusters from the kmeans model. It also gives me the opportunity to use the threejs package which we can interact with.

#library(rgl)
#options(rgl.printRglwidget = TRUE)
#plot3d(beta_alphas_sharpe, col = kmeans_model$cluster, size=2, type='s')

library(threejs)
scatterplot3js(beta_alphas_sharpe, color = factor(kmeans_model$cluster), size = 0.5, bg = "black")

Factor Modelling Fama and French

Finally, I analyse the performance of the various ETF’s. The CAPM formula tries to explain the performance of a portfolio through a sigle factor, the market as a whole. CAPM is defined as follows:

\[R_{a} = R_{rf} + \beta_{a}(R_{m} - R_{rf})\]

We can take the model further by adding Factors to the model. Thus, the 3 Factor model looks like:

\[R_{a} = R_{rf} + \beta_{a}(R_{m} - R_{rf}) + \beta_{b}(SMB) + \beta_{c} HML\]

Which takes the market factor from the CAPM and adds 2 new factors, the \(SMB\) and \(HML\) or Small-minus-Big and High-minus-Low. The \(SMB\) tries to caputre the size effect, small market cap companies should outperform the market over the long run. The \(HML\) aims to capture the value verses growth effect by sorting assests based on high book-to-market ratios (defined at the begining of the post). Firms which typically have a high book-to-market ratio are value stocks and low book-to-market ratios are growth stocks, the literature also shows that value stocks outperform growth stocks over the long term horizon.

The Fama French factor model allows us to analyse fund managers performance, since, if an asset managers portfolio can be explained by the Fama French factors then that fund manager has not added any value through portfolio selection and skill and subseuently obtained no alpha.

More formally, our equation becomes:

\[X^T = a^T + B F^T + E^T\]

Where \(F\) is a matrix of factors and \(B\) is a scalar of \(\beta\)’s such that \(B = [{b_{1}, b_{2}, b_{3}]}\), in which \(b_{1}\) is the market \(\beta\), \(b_{2}\) is the \(SMB\) \(\beta\) and \(b_{3}\) is the \(HML\) \(\beta\). More formally extending to include more factors up to \(N\) we have; \(B = [\beta_{1}, ... , \beta_{N}]^T\)

We want to fit a line which minimises the squared errors such that; \(minimise, \beta\) \(|X^T - a^T-BF^T|_{F}^2\)

The solution becomes \([a, B] = X^T \widetilde{F}(\widetilde{F}^T\widetilde{F})^{-1}\) which we can solve in R using the following:

(For completness I re-download the data)

  1. I download the data as before randomly selecting from Wikipedia and convert the daily prices to daily returns - (I set a seed such that we collect the same data using set.seed).
library(quantmod)
library(tidyquant)
library(timetk)
library(rvest)

start_date <- "2016-01-01"
end_date <- "2019-11-15"


set.seed(4746)
url <- "https://en.wikipedia.org/wiki/List_of_S%26P_500_companies"

symbols <- url %>%
  read_html() %>%
  html_nodes(xpath = '//*[@id="constituents"]') %>% 
  html_table() %>% 
  .[[1]] %>% 
  filter(!str_detect(Security, "Class A|Class B|Class C")) %>%     # Removes firms with Class A shares
  sample_n(20) %>% 
  pull(Symbol)


asset_returns <- tq_get(
  symbols,
  from = start_date,
  to = end_date
) %>% 
  group_by(symbol) %>% 
  tq_transmute(
    select = adjusted,
    mutate_fun = periodReturn,
    period = "daily",
    type = "arithmetic"
  ) %>% 
  select(symbol, date, daily.returns) %>% 
  pivot_wider(names_from = symbol, values_from = daily.returns) %>% 
  .[-1, ]
  1. I download the ETFs just as before and convert to daily returns.
myETFs <- c("SPY", "IVV", "VTI", "VOO", "QQQ", "VEA", "IEFA", "AGG", "VWO",
            "EFA","IEMG","VTV", "IJH", "IWF","BND", "IJR", "IWM", "VUG", 
            "GLD", "IWD", "VIG", "VNQ", "USMV", "LQD", "VO", "VYM", "EEM",
            "VB", "VCSH", "XLF", "VCIT", "VEU", "XLK", "ITOT", "IVW", "BNDX",
            "VGT", "DIA", "BSV", "SHV", "IWB", "IWR", "TIP", "SCHF", "MBB", "SDY",
            "MDY", "SCHX", "IEF", "HYG", "DVY", "XLV", "SHY", "IXUS", "TLT", "IVE",
            "PFF", "IAU", "VXUS", "RSP", "SCHB", "VV", "GOVT", "EMB", "MUB", "QUAL",
            "XLY", "VBR", "EWJ", "XLP", "VGK", "SPLV", "MINT", "BIV", "IGSB", "EFAV",
            "VT", "GDX", "XLU", "IWS", "XLI", "SCHD", "IWP", "ACWI", "VMBS", "XLE", "JNK",
            "VOE", "FLOT", "IWV", "JPST", "SCZ", "IEI", "IWN", "DGRO", "VBK", "IGIB", "IWO")


etf_returns <- tq_get(
  myETFs,
  from = start_date,
  to = end_date
) %>% 
  group_by(symbol) %>% 
  tq_transmute(
    select = adjusted,
    mutate_fun = periodReturn,
    period = "daily",
    type = "arithmetic"
  ) %>% 
  select(symbol, date, daily.returns) %>% 
  pivot_wider(names_from = symbol, values_from = daily.returns) %>% 
  .[-1, ]
  1. I take the average daily return of the randomly selected stocks and join the data with the ETFs and set the data as a time series object. I also download the daily Fama French 3 factors from Kenneth French’s website and clean the data up a little.
portfolio_and_etfs <- asset_returns %>% 
  mutate(myPortfolio = rowMeans(select(., -date), na.rm = TRUE)) %>% 
  select(date, myPortfolio) %>% 
  inner_join(etf_returns, by = c("date")) %>% 
  tk_xts(date_var = date)

# Daily fama french
library(glue)
library(readr)
temp <- tempfile()
base <- "http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/"
factor <- "F-F_Research_Data_Factors_daily"
format<-"_CSV.zip"
full_url <-glue(base,factor,format,sep ="")
download.file(full_url,temp,quiet = TRUE)

North_America_3_Factors <- read_csv(unz(temp, "F-F_Research_Data_Factors_daily.CSV"), skip = 4) %>%
  rename(date = X1) %>% 
  mutate(date = ymd(date)) %>% 
  drop_na(date) %>%
  rename(Mkt_Rf_3 = `Mkt-RF`,
         SMB_3 = SMB,
         HML_3 = HML,
         RF_3 = RF) %>% 
  select(-RF_3) %>% 
  mutate_at(vars(c("Mkt_Rf_3", "SMB_3", "HML_3")), funs(./100)) %>% 
  tk_xts(date_var = date)

Finally we can compute the solution \([a, B] = X^T \widetilde{F}(\widetilde{F}^T\widetilde{F})^{-1}\)

three_factors <- North_America_3_Factors[index(North_America_3_Factors) %in% c(min(index(portfolio_and_etfs)):max(index(North_America_3_Factors)))] # this ensures that the dates from the both time series match
portfolio_and_etfs_new <- portfolio_and_etfs[index(portfolio_and_etfs) %in% c(min(index(three_factors))):max(index(three_factors))] # There is probably a much nicer way to go about doing this but this works...
FF_3Factors <- cbind(alpha = 1, three_factors)
Gamma <- t(solve(t(FF_3Factors) %*% FF_3Factors, t(FF_3Factors) %*% portfolio_and_etfs_new))

Which we can see as:

alpha Mkt_Rf_3 SMB_3 HML_3
myPortfolio 0.0001697 0.9844870 0.0155642 0.0629606
SPY 0.0000391 0.9784199 -0.1331473 -0.0085326
IVV 0.0000423 0.9771912 -0.1341942 -0.0113980
VTI 0.0000440 0.9727841 -0.0164012 -0.0019078
VOO 0.0000430 0.9769683 -0.1356460 -0.0099922
QQQ 0.0000440 1.1645777 -0.1440829 -0.4955734

Instead, we can estimate these using linear regression models. For my random portfolio myPortfolio we can use the lm function to fit a linear model and then use the tidy function from the broom package to tidy the output up a little:

library(broom)
regdata <- cbind(FF_3Factors, portfolio_and_etfs_new)
#sapply(regdata, function(x) sum(is.na(x)))   # One column was causing problems since it consisted of just NA values. I remove it below.
regdata$JPST <- NULL      # This ETF failed to download or does not exist on the Yahoo Finance site

myPortfolioReg <- lm(myPortfolio ~ Mkt_Rf_3 + SMB_3 + HML_3, data = regdata, na.action = na.exclude)
tidy(myPortfolioReg)
## # A tibble: 4 x 5
##   term        estimate std.error statistic p.value
##                          
## 1 (Intercept) 0.000170  0.000112     1.51  0.131  
## 2 Mkt_Rf_3    0.984     0.0136      72.6   0      
## 3 SMB_3       0.0156    0.0226       0.690 0.490  
## 4 HML_3       0.0630    0.0200       3.15  0.00167

We can apply this to all our ETFs in the data using the apply command and applying our own custom lm function.

regressionLists <- apply(regdata, 2, function(y) lm(y ~ Mkt_Rf_3 + SMB_3 + HML_3, data = regdata, na.action = na.exclude))

We can also apply the tidy command to individual ETFs and then use the stars.pval to make the data even more tidy awesome.

library(gtools)
tidy(regressionLists$myPortfolio) %>% 
  mutate(p.value = stars.pval(p.value))
## # A tibble: 4 x 5
##   term        estimate std.error statistic p.value
##                          
## 1 (Intercept) 0.000170  0.000112     1.51  " "    
## 2 Mkt_Rf_3    0.984     0.0136      72.6   ***    
## 3 SMB_3       0.0156    0.0226       0.690 " "    
## 4 HML_3       0.0630    0.0200       3.15  **
tidy(regressionLists$VTV)%>% 
  mutate(p.value = stars.pval(p.value))
## # A tibble: 4 x 5
##   term          estimate std.error statistic p.value
##                            
## 1 (Intercept)  0.0000514 0.0000446      1.15 " "    
## 2 Mkt_Rf_3     0.918     0.00539      170.   ***    
## 3 SMB_3       -0.131     0.00897      -14.7  ***    
## 4 HML_3        0.265     0.00794       33.4  ***
tidy(regressionLists$LQD) %>% 
  mutate(p.value = stars.pval(p.value))
## # A tibble: 4 x 5
##   term         estimate std.error statistic p.value
##                           
## 1 (Intercept)  0.000241 0.0000882      2.74 **     
## 2 Mkt_Rf_3    -0.0119   0.0106        -1.12 " "    
## 3 SMB_3       -0.0325   0.0177        -1.84 .      
## 4 HML_3       -0.157    0.0157       -10.0  ***

Finally, we can apply this same method to all our ETFs using the lapply function to tidy the data and the map function to mutate or convert the p-value to stars. Then take a random sample of 5 ETFs regressions.

library(purrr)
myTidyRegressions <- lapply(regressionLists[4:length(regressionLists)], tidy)
myTidiedRegressions <- map(myTidyRegressions, ~mutate(., p.value = stars.pval(p.value)))
library(rlist)
list.sample(myTidiedRegressions, 5)
## $DGRO
## # A tibble: 4 x 5
##   term         estimate std.error statistic p.value
##                           
## 1 (Intercept)  0.000154 0.0000507      3.04 **     
## 2 Mkt_Rf_3     0.890    0.00612      146.   ***    
## 3 SMB_3       -0.145    0.0102       -14.2  ***    
## 4 HML_3        0.121    0.00902       13.4  ***    
## 
## $EFAV
## # A tibble: 4 x 5
##   term           estimate std.error statistic p.value
##                             
## 1 (Intercept) -0.00000220  0.000142   -0.0155 " "    
## 2 Mkt_Rf_3     0.572       0.0171     33.5    ***    
## 3 SMB_3       -0.112       0.0285     -3.94   ***    
## 4 HML_3       -0.0172      0.0252     -0.684  " "    
## 
## $XLE
## # A tibble: 4 x 5
##   term         estimate std.error statistic p.value
##                           
## 1 (Intercept) -0.000264  0.000270    -0.978 " "    
## 2 Mkt_Rf_3     1.09      0.0326      33.5   ***    
## 3 SMB_3        0.0920    0.0543       1.69  .      
## 4 HML_3        0.716     0.0481      14.9   ***    
## 
## $IVW
## # A tibble: 4 x 5
##   term          estimate std.error statistic p.value
##                            
## 1 (Intercept)  0.0000434 0.0000407      1.07 " "    
## 2 Mkt_Rf_3     1.01      0.00491      206.   ***    
## 3 SMB_3       -0.180     0.00817      -22.0  ***    
## 4 HML_3       -0.312     0.00724      -43.1  ***    
## 
## $VEA
## # A tibble: 4 x 5
##   term         estimate std.error statistic p.value
##                           
## 1 (Intercept) -0.000113  0.000154    -0.730 " "    
## 2 Mkt_Rf_3     0.852     0.0186      45.8   ***    
## 3 SMB_3       -0.0850    0.0310      -2.74  **     
## 4 HML_3        0.150     0.0274       5.45  ***

A few notes here: We should be modelling the excess returns on the ETFs and not just the ETF returns. This is simple enough by replacing for example myPortfolio in the lm regressions with \((myPortfolio - RF_3)\) where \(RF_3\) is the risk free rate which comes with the Fama and French data. We should also probably use \(NeweyWest\) adjusted standard errors by running coeftest(myPortfolioReg, NeweyWest(myPortfolioReg, lag = N, prewhite = FALSE)). Literature suggests that lags of \(4(N/100)^{2/9}\) should be used, where \(N\) is the number of observations.
We could rank the ETFs based on their alphas as before and go long on the high alphas and short on the low alphas. Run our hedge portfolio through a Fama French regression here and see if we are able to obtain better performances.

Finally, here I only used the 3 Factor model. There are more factors from the literature we could use. On the Kenneth French website we can collect data on \(Market\), \(SMB\), \(HML\), \(RMW\), \(CMA\) and \(MOM\). Where the \(RMW\) factor is a profitability factor, the \(CMA\) is an investment factor and the \(MOM\) is a momentum factor.

Below I pull in the daily Fama and French 5 Factor model and plot them.

library(glue)
library(timetk)
library(plotly)
library(pipeR)

temp <- tempfile()
base <- "http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/"
factor <- "North_America_5_Factors"
format<-"_CSV.zip"
full_url <-glue(base,factor,format,sep ="")
download.file(full_url,temp,quiet = TRUE)

North_America_Factors <- read_csv(unz(temp, "North_America_5_Factors.csv"), skip = 6) %>%
  rename(date = X1) %>%
  mutate_at(vars(-date), as.numeric) %>%
  mutate(date = rollback(ymd(parse_date_time(date, "%Y%m") + months(1)))) %>%
  drop_na(date) %>%
  rename(Mkt_Rf = `Mkt-RF`) %>% 
  select(-RF) %>% 
  mutate_at(vars(c("Mkt_Rf", "SMB", "HML", "RMW", "CMA")), funs(./100)) %>% 
  tk_xts(date_var = date)

ax <- list(
  zeroline = TRUE,
  showline = TRUE,
  mirror = "ticks",
  gridcolor = toRGB('white'),
  gridwidth = 2,
  zerolinecolor = toRGB('white'),
  zerolinewidth = 4,
  linecolor = toRGB('white'),
  linewidth = 6,
  color = 'white'
)


North_America_Factors %>>% 
  (cumprod(1+.)) %>>%
  ROC(36, type = "discrete") %>>%
  na.omit() %>>% 
  (
    plot_ly(
      x = colnames(.),
      y = as.Date(index(.)),
      z = data.matrix(.),
      type = "surface",
      colors = c("darkblue", "yellow", "darkred")
    )
  ) %>%
  plotly::layout(
    xaxis = ax, yaxis = ax,
    paper_bgcolor ='rgb(0,0,0)',
    #plot_bgcolor ='rgb(0,0,0)',
    title = "Fama French Factors",
    scene = list(
      xaxis = c(ax, title = "Factors", tickangle = -45),
      yaxis = c(ax, title = "Date", tickangle = -45),
      zaxis = c(ax, title = "Fama French Factor")
    ),
    titlefont = list(
      family = "Agency FB",
      size = 50,
      color = 'white')
  )

Sharpe, William F. 1964. “Capital Asset Prices: A Theory of Market Equilibrium Under Conditions of Risk.” The Journal of Finance 19 (3): 425–42.

———. 1966. “Mutual Fund Performance.” The Journal of Business 39 (1): 119–38.

To leave a comment for the author, please follow the link and comment on their blog: Posts on Matthew Smith R Shenanigans.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)