Today, we will return to the Fama French (FF) model of asset returns and use it as a proxy for fitting and evaluating multiple linear models. In a previous post, we reviewed how to run the FF three-factor model on the returns of a portfolio. That is, we ran one model on one set of returns. Today, we will run multiple models on multiple streams of returns, which will allow us to compare those models and hopefully build a code scaffolding that can be used when we wish to explore other factor models.
Let’s get to it!
We will start by importing daily prices and calculating the returns of our five usual ETFs: SPY, EFA, IJS, EEM and AGG. I covered the logic for this task in a previous post and the code is below.
library(tidyverse)
library(broom)
library(tidyquant)
library(timetk)
symbols % reduce(merge) %__%
`colnames%
gather(asset, returns, -date) %__% group_by(asset) %__% mutate(returns = (log(returns) - log(lag(returns)))) %__% na.omit()
asset_returns_long %__% head()
# A tibble: 6 x 3
# Groups: asset [1]
date asset returns
1 2013-01-02 SPY 0.0253 2 2013-01-03 SPY -0.00226
3 2013-01-04 SPY 0.00438
4 2013-01-07 SPY -0.00274
5 2013-01-08 SPY -0.00288
6 2013-01-09 SPY 0.00254
We now have the returns of our five ETFs saved in a tidy tibble called asset_returns_long. Normally we would combine these into one portfolio, but we will leave them as individual assets today so we can explore how to model the returns of multiple assets saved in one data object.
We are going to model those individual assets on the Fama French set of five factors to see if those factors explain our asset returns. Those five factors are the market returns (similar to CAPM), firm size (small versus big), firm value (high versus low book-to-market), firm profitability (high versus low operating profits), and firm investment (high versus low total asset growth). To learn more about the theory behind using these factors, see this article.
Our next step is to import our factor data, and luckily FF make them available on their website. I presented the code showing how to import this data in a previous post on my blog so I won’t go through the steps again, but the code is below.
factors_data_address % select(-RF)
Global_5_Factors %__% head()
# A tibble: 6 x 6
date MKT SMB HML RMW CMA
1 1990-07-02 0.00700 -0.000600 -0.0031 0.0022 0.0004 2 1990-07-03 0.0018 0.0008 -0.0017 0.0007 0.0004 3 1990-07-04 0.0063 -0.0019 -0.0016 -0.0007 0.000300
4 1990-07-05 -0.0074 0.0031 0.0017 -0.0013 0.000600
5 1990-07-06 0.002 -0.0015 0.0002 0.002 -0.000300
6 1990-07-09 0.00240 0.0018 0.0001 0.0004 -0.00240 Next, we mash asset_returns_long and Global_5_Factors into one data object, using left_join(..., by = &quot;date&quot;) because each object has a column called date.
data_joined_tidy %
left_join(Global_5_Factors, by = &quot;date&quot;) %__% na.omit()
data_joined_tidy %__% head(5)
# A tibble: 5 x 8
# Groups: asset [1]
date asset returns MKT SMB HML RMW CMA
1 2013-01-02 SPY 0.0253 0.0199 -0.0043 0.0028 -0.0028 -0.0023 2 2013-01-03 SPY -0.00226 -0.0021 0.00120 0.000600 0.0008 0.0013 3 2013-01-04 SPY 0.00438 0.0064 0.0011 0.0056 -0.0043 0.0036 4 2013-01-07 SPY -0.00274 -0.0014 0.00580 0 -0.0015 0.0001 5 2013-01-08 SPY -0.00288 -0.0027 0.0018 -0.00120 -0.0002 0.00120
We now have a tibble called data_joined_tidy with the asset names in the asset column, asset returns in the returns column, and our five FF factors in the MKT, SMB, HML, RMW and CMA columns. We want to explore whether there is a linear relationship between the returns of our five assets and any/all of the five FF factors. Specifically, we will run several linear regressions, save the results, examine the results, and then quickly visualize the results.
The broom and purrr packages will do a lot of the heavy lifting for us eventually, but let’s start with a simple example: regress the return of one ETF on one of the FF factors. We will use do() for this, which I believe has been declared not a best practice but it’s so irresistibly simple to plunk into the pipes.
data_joined_tidy %__% filter(asset == &quot;SPY&quot;) %__% do(model = lm(returns ~ MKT, data = .))
Source: local data frame [1 x 2]
Groups: # A tibble: 1 x 2
asset model * 1 SPY That worked, but our model is stored as a nested S3 object. Let’s use tidy(), glance(), and augment() to view the results.
data_joined_tidy %__% filter(asset == &quot;SPY&quot;) %__% do(model = lm(returns ~ MKT, data = .)) %__% tidy(model)
# A tibble: 2 x 6
# Groups: asset [1]
asset term estimate std.error statistic p.value
1 SPY (Intercept) 0.000126 0.000100 1.26 0.208
2 SPY MKT 1.02 0.0154 66.0 0 data_joined_tidy %__% filter(asset == &quot;SPY&quot;) %__% do(model = lm(returns ~ MKT, data = .)) %__% glance(model)
# A tibble: 1 x 12
# Groups: asset [1]
asset r.squared adj.r.squared sigma statistic p.value df logLik
1 SPY 0.776 0.776 0.00354 4356. 0 2 5319.
# ... with 4 more variables: AIC , BIC , deviance ,
# df.residual data_joined_tidy %__% filter(asset == &quot;SPY&quot;) %__% do(model = lm(returns ~ MKT, data = .)) %__% augment(model) %__% head(5)
# A tibble: 5 x 10
# Groups: asset [1]
asset returns MKT .fitted .se.fit .resid .hat .sigma .cooksd
1 SPY 0.0253 0.0199 0.0204 3.17e-4 4.90e-3 7.99e-3 0.00354 7.77e-3
2 SPY -0.00226 -0.0021 -0.00201 1.07e-4 -2.47e-4 9.17e-4 0.00354 2.24e-6
3 SPY 0.00438 0.0064 0.00665 1.36e-4 -2.27e-3 1.47e-3 0.00354 3.02e-4
4 SPY -0.00274 -0.0014 -0.00130 1.04e-4 -1.44e-3 8.59e-4 0.00354 7.07e-5
5 SPY -0.00288 -0.0027 -0.00263 1.11e-4 -2.56e-4 9.82e-4 0.00354 2.56e-6
# ... with 1 more variable: .std.resid We can quickly pop our augmented results into ggplot() and inspect our residuals versus our fitted values. The important takeaway here is that our augmented results are in a data frame, so we can use all of our ggplot() code flows to create visualizations.
data_joined_tidy %__% filter(asset == &quot;SPY&quot;) %__% do(model = lm(returns ~ MKT, data = .)) %__% augment(model) %__% ggplot(aes(y = .resid, x = .fitted)) +
geom_point(color = &quot;cornflowerblue&quot;)
Alright, we have run a simple linear regression and seen how tidy(), glance(), and augment() clean up the model results. We could, of course, keep repeating this process for any combination of the FF factors and any of our ETFs, but let’s look at a more efficient approach for fitting multiple models on all of our ETFs.
First, let’s create and save three models as functions. That will allow us to pass them to the map functions from purrr. We will save a one-factor model as a function called model_one, a three-factor model as a function called model_two and a five-factor model as a function called model_three. Note that each function takes a data frame as an argument.
model_one % group_by(asset) %__% nest() %__%
mutate(one_factor_model = map(data, model_one),
three_factor_model = map(data, model_two),
five_factor_model = map(data, model_three))
# A tibble: 5 x 5
asset data one_factor_model three_factor_mod… five_factor_mod…
1 SPY % nest() %__%
mutate(one_factor_model = map(data, model_one),
tidied_one = map(one_factor_model, tidy)) %__% unnest(tidied_one)
# A tibble: 10 x 6
asset term estimate std.error statistic p.value
1 SPY (Intercept) 0.000126 0.000100 1.26 2.08e- 1
2 SPY MKT 1.02 0.0154 66.0 0. 3 EFA (Intercept) -0.000288 0.0000972 -2.97 3.07e- 3
4 EFA MKT 1.29 0.0150 86.0 0. 5 IJS (Intercept) 0.0000634 0.000171 0.371 7.11e- 1
6 IJS MKT 1.13 0.0264 42.9 3.25e-248
7 EEM (Intercept) -0.000491 0.000203 -2.42 1.57e- 2
8 EEM MKT 1.40 0.0313 44.8 4.94e-263
9 AGG (Intercept) 0.0000888 0.0000572 1.55 1.21e- 1
10 AGG MKT -0.0282 0.00883 -3.19 1.46e- 3
Let’s use glance() and augment() as well.
data_joined_tidy %__% group_by(asset) %__% nest() %__%
mutate(one_factor_model = map(data, model_one),
tidied_one = map(one_factor_model, tidy),
glanced_one = map(one_factor_model, glance),
augmented_one = map(one_factor_model, augment))
# A tibble: 5 x 6
asset data one_factor_model tidied_one glanced_one augmented_one 1 SPY