Using loops

Loops are great and (almost) everyone loves then. While they can be a bit more verbose than fancy on-liners, the structure of a loop is very flexible and this can help solve complex problems. Let use it in our problem.

The first step in using loops is the understand the vector that will be used as iterator in the loop. In our problem we are processing each stock, so the number of iterations in the loop is simply the number of stocks in the sample. We can find the unique stocks with the command unique.

# Check unique tickers
unique.tickers <- unique(df.stocks$ticker)

# create a empty vector to store betas
beta.vec <- c()

for (i.ticker in unique.tickers){
  
  # message to prompt
  cat('\nRunning ols for',i.ticker)
  
  # filter the data.frame for stock i.ticker
  df.temp <- df.stocks[df.stocks$ticker==i.ticker, ]
  
  # calculate beta with lm
  my.ols <- lm(data = df.temp, formula = ret.adjusted.prices ~ ret.MktIdx)
  
  # save beta
  my.beta <- coef(my.ols)[2]
  
  # store beta em beta.vec
  beta.vec <- c(beta.vec, my.beta)
}

## 
## Running ols for ZBH
## Running ols for OXY
## Running ols for C
## Running ols for XLNX
## Running ols for VZ
## Running ols for BEN
## Running ols for WY
## Running ols for ROL
## Running ols for VTR
## Running ols for ^GSPC

# print result
print(data.frame(unique.tickers,beta.vec))

##    unique.tickers  beta.vec
## 1             ZBH 0.9130036
## 2             OXY 1.0046153
## 3               C 1.3468402
## 4            XLNX 1.2272195
## 5              VZ 0.5864176
## 6             BEN 1.2555433
## 7              WY 0.9400250
## 8             ROL 0.8057286
## 9             VTR 0.4923402
## 10          ^GSPC 1.0000000

As you can see, the result is a lengthy code, but it works quite well. The final result is a dataframe with the tickers and their betas. Notice that, as expected, the betas are all positive and ^GSPC has a beta equal to 1.

Using function by

Another way of solving the problem is to calculate the betas using one of the functions from the apply family. In this case, we will use function by. Be aware that you can also solve the problem using tapply and lapply. The code, however, will increase in complexity.

The function by works similarly to tapply. The difference is that it is oriented to dataframes. That is, given a grouping variable, the original dataframe is broken into smaller dataframes and each piece is passed to a function. This helps a lot our problem since we need to work with two columns, the vector of returns of the asset and the vector of returns of the market index.

Given the functional form of by, will need to encapsulate a procedure that takes a dataframe as input and returns a coefficient beta, calculated from columns ret and ret.MktIdx. The next code does that.

get.beta <- function(df.temp){
  
  # estimate model
  my.ols <- lm(data=df.temp, formula = ret.adjusted.prices ~ ret.MktIdx)
  
  # isolate beta
  my.beta <- coef(my.ols)[2]
  
  # return beta
  return(my.beta)
}

The previous function accepts a single dataframe called df.temp, uses it to calculate a linear model with function lm and then returns the resulting beta, which is the second coefficient in coef(my.ols). Now, lets use it with function by.

# get betas with by
my.l <- by(data = df.stocks, 
           INDICES = df.stocks$ticker, 
           FUN = get.beta)

# my.l is an objetct of class by. To get only its elements, we can unclass it
betas <- unclass(my.l)

# print result
print(data.frame(betas))

##           betas
## ^GSPC 1.0000000
## BEN   1.2555433
## C     1.3468402
## OXY   1.0046153
## ROL   0.8057286
## VTR   0.4923402
## VZ    0.5864176
## WY    0.9400250
## XLNX  1.2272195
## ZBH   0.9130036

Again, it worked well. Needless to say that the results are identical to the previous case.

Using dplyr

Now, let’s solve our problem using package dplyr. If you are not familiar with the tidyverse and the work of Hadley Wickham, you will be a happier person after reading the rest of this post, trust me.

Package dplyr is one of my favorites and most used packages. It allows for the representation of data processing procedures in a simpler and more intuitive way. It really helps to tackle computational problems if you can fit it within a flexible structure. This is what, in my opinion, dplyr does best. It combines clever functions with dataframes in the long (tidy) format.

Have a look in the next set of code.

library(dplyr)

beta.tab <- df.stocks %>% 
  group_by(ticker) %>% # group by column ticker
  do(ols.model = lm(data = ., formula = ret.adjusted.prices ~ret.MktIdx)) %>%   # estimate model
  mutate(beta = coef(ols.model)[2]) # get coefficients

print(beta.tab)

## Source: local data frame [10 x 3]
## Groups: <by row>
## 
## # A tibble: 10 x 3
##    ticker ols.model  beta
##    <chr>  <list>    <dbl>
##  1 ^GSPC  <lm>      1    
##  2 BEN    <lm>      1.26 
##  3 C      <lm>      1.35 
##  4 OXY    <lm>      1.00 
##  5 ROL    <lm>      0.806
##  6 VTR    <lm>      0.492
##  7 VZ     <lm>      0.586
##  8 WY     <lm>      0.940
##  9 XLNX   <lm>      1.23 
## 10 ZBH    <lm>      0.913

After loading dplyr, we use the pipeline operator %>% to streamline all calculations. This means that we don’t need to keep a copy of intermediate calculations. Also, It looks pretty, don’t you agree?

The line beta.tab <- df.stocks %>% passes the dataframe df.stocks for the next line, group_by(ticker) %>%, which will group the dataframe according the the values of column ticker and pass the result for the next step. The line do(ols.model = lm(data = ., formula = ret ~ret.MktIdx)) estimates the model by passing a temporary dataframe and saves it in a column called ols.model. Notice that the model is a S3 object and not a single value. The dataframe alternative tibble is flexible with its content. The final line, mutate(beta = coef(ols.model)[2]) retrieves the beta from each element of the column ols.model.

What I really like about dplyr is that it makes it easy to extend the original code. As an example, if I wanted to use a second grouping variable, I can just add it in the second line as group_by(ticker, newgroupingvariable). This becomes handy if you need, lets say, to estimated the model in different time periods.

As an example, let’s assume that I want to split the sample for each stock in half and see if the betas change significantly from time period to the other. This robustness check is a very common procedure in scientific research. First, let’s build a new column in df.stocks that sets the time periods as Sample 1 and Sample 2. We can use tapply for that;

set.sample <- function(ref.dates){
  my.n <- length(ref.dates) # get number of dates
  my.n.half <- floor(my.n/2) # get aproximate half of observations
  
  # create grouping variable
  samples.vec <- c(rep('Sample 1', my.n.half ), rep('Sample 2', my.n-my.n.half))
  
  # return
  return(samples.vec)
}

# build group
my.l <- tapply(X = df.stocks$ref.date, 
               INDEX = df.stocks$ticker,
               FUN = set.sample )

# unsort it
my.l <- my.l[my.tickers]

# save it in dataframe
df.stocks$my.sample <- unlist(my.l)

We proceed by calling the same functions as before, but using an additional grouping variable.

beta.tab <- df.stocks %>% 
  group_by(ticker,my.sample) %>% # group by column ticker
  do(ols.model = lm(data = ., formula = ret.adjusted.prices ~ret.MktIdx)) %>%   # estimate model
  mutate(beta = coef(ols.model)[2]) # get coefficients

print(beta.tab)

## Source: local data frame [20 x 4]
## Groups: <by row>
## 
## # A tibble: 20 x 4
##    ticker my.sample ols.model  beta
##    <chr>  <chr>     <list>    <dbl>
##  1 ^GSPC  Sample 1  <lm>      1.00 
##  2 ^GSPC  Sample 2  <lm>      1.   
##  3 BEN    Sample 1  <lm>      1.40 
##  4 BEN    Sample 2  <lm>      1.08 
##  5 C      Sample 1  <lm>      1.54 
##  6 C      Sample 2  <lm>      1.11 
##  7 OXY    Sample 1  <lm>      1.10 
##  8 OXY    Sample 2  <lm>      0.882
##  9 ROL    Sample 1  <lm>      0.781
## 10 ROL    Sample 2  <lm>      0.836
## 11 VTR    Sample 1  <lm>      0.653
## 12 VTR    Sample 2  <lm>      0.295
## 13 VZ     Sample 1  <lm>      0.669
## 14 VZ     Sample 2  <lm>      0.485
## 15 WY     Sample 1  <lm>      1.08 
## 16 WY     Sample 2  <lm>      0.772
## 17 XLNX   Sample 1  <lm>      1.07 
## 18 XLNX   Sample 2  <lm>      1.42 
## 19 ZBH    Sample 1  <lm>      0.888
## 20 ZBH    Sample 2  <lm>      0.944

As we can see, the output now shows the beta for all combinations between ticker and my.sample. It seems that the betas tend to be higher for Sample 2, meaning that the overall systematic risk in the market has increased over time, at least for the majority of the ten chosen stocks. Given the small sample of stocks, It might be interesting to test for this property in a larger dataset.

Back to the model, if you want more information about it, you can just write new lines in the last call to %>%. Let’s say, for example, that you want to get the value of alpha and the corresponding t-statistic of both coefficients. We can use the following code for that:

library(dplyr)

beta.tab <- df.stocks %>% 
  group_by(ticker) %>% # group by column ticker
  do(ols.model = lm(data = ., formula = ret.adjusted.prices ~ret.MktIdx)) %>%   # estimate model
  mutate(beta = coef(ols.model)[2],
         beta.tstat = summary(ols.model)[[4]][2,3],
         alpha = coef(ols.model)[1],
         alpha.tstat = summary(ols.model)[[4]][1,3]) # get coefficients

## Warning in summary.lm(ols.model): essentially perfect fit: summary may be
## unreliable

## Warning in summary.lm(ols.model): essentially perfect fit: summary may be
## unreliable

print(beta.tab)

## Source: local data frame [10 x 6]
## Groups: <by row>
## 
## # A tibble: 10 x 6
##    ticker ols.model  beta beta.tstat     alpha alpha.tstat
##    <chr>  <list>    <dbl>      <dbl>     <dbl>       <dbl>
##  1 ^GSPC  <lm>      1        3.49e17  0              0    
##  2 BEN    <lm>      1.26     2.99e 1 -0.000708      -1.96 
##  3 C      <lm>      1.35     3.46e 1 -0.000184      -0.548
##  4 OXY    <lm>      1.00     2.19e 1 -0.000217      -0.551
##  5 ROL    <lm>      0.806    2.11e 1  0.000879       2.67 
##  6 VTR    <lm>      0.492    8.94e 0 -0.000144      -0.304
##  7 VZ     <lm>      0.586    1.65e 1  0.000287       0.935
##  8 WY     <lm>      0.940    2.24e 1 -0.000463      -1.28 
##  9 XLNX   <lm>      1.23     2.46e 1  0.000623       1.45 
## 10 ZBH    <lm>      0.913    2.17e 1 -0.000164      -0.451

In the previous code, I added line beta.tstat = summary(ols.model)[[4]][2,3] that returns the t-statistic of the beta coefficient. The location of this parameter is found by investigating the elements of an object of type lm. After calling summary, the t-statistic is available in the fourth element of the lm object, which is a matrix with several information from the estimation. The t-statistic for the alpha parameter is found in a similar way.

As you can see, the syntax of dplyr make it easy to extend the model and quickly try new things. It is possible to do the same using other R functions and a loop, but using dplyr is really handy.