# R: Fast and efficient way of running a multivariate regression across a (really) large panel (First pass of Fama MacBeth)

I am attempting to run a rolling multivariate regression (14 explanatory variables) across a panel of 5000 stocks:

• For each of the 5000 stocks, I run 284 regressions (by rolling over my sample period).
• In summary: 1,420,000 regressions in total are ran for the panel.

To achieve this, I make use a nested "for loop": loop over securities and over time. Coefficients are exported to a csv file.

As expected, the issue is that the entire procedure takes a HUGE amount of time to complete. Would there be an efficient way of handling this? (As I realize that the "apply" function is more efficient than a "for loop", please keep in mind that given the huge processing time, the time gain from the alternative use of the "apply" function would still be minimal).

Here is a snapshot of the code:

sec = ncol(ret.zoo)
num.factors = ncol(data)
rows = nrow(ret.zoo) - 60 + 1
col.names <- c("gvkey", "date", "intercept", colnames(data))
write.table(as.data.frame(t(col.names)), file = paste(path, "betas.csv", sep = ""),  row.names = FALSE, col.names = FALSE, sep = ",")

for(i in 1:sec) {
beta = data.frame(matrix(nc = num.factors + 3, nr = rows))
df = merge(ret.zoo[,i], data)
names(df) <- c("return", names(data))

for(j in 1:rows) {
#Checks if number of observations >=30. If so, regression is ran. Otherwise, it is not.
no.na = ret.zoo[j:(j+59),i][which(!is.na(coredata(ret.zoo[j:(j+59),i])))]
if(length(no.na) >= 30) {
beta[j,1] = substr(colnames(ret.zoo)[i],2,7)
beta[j,2] = as.character(index(df[(j+59),])) ### Date
beta[j,3:(num.factors+3)] = coef(lm(return ~., data = as.data.frame(df[j:(j+59),]), na.action = na.omit))
}
}
write.table(beta, file = paste(path, "betas.csv", sep = ""), append = T, sep = ",", row.names = FALSE, col.names = FALSE)
rm(beta)
}


Note that:

• sec: number of stocks (securities). Each security has a time series of returns.
• rows: number of time periods (over which we roll the regression)
• beta: matrix of coefficients of all regressions for each security. It is cleared every time for each sec.

## MODEL:

Here is the regression model for each security i at time t :

R(i,t) = a(i,t) + b1(i,t)f1(t) + b2(i,t)f2(t) + .... + bn(i,t)fn(t) + e(i,t)

where b are the regression coefficients, f the factors, and e the residuals.

Note that i is in [1:5000], the number of factors n is 14, and time t is in [1:343] (343 months).

• For each security i, we run this regression over rolling periods of 60 months (hence the j:j+59 in R code).

• Each rolling regression is ran only if the non-NA number of observations of the rolling window for the dependent variable is >= 30 (While the independent variables cannot be NA, the dependent variables (here stock returns) can take NA values, if the stock drops from the index).

• We then obtain 284 = 343 - 60 + 1 beta coefficients for each factor f for each security i. These are stored in the "beta" dataframe (the "beta" dataframe has nr = 284, and ncol = 14+3 (14 factors, intercept, date, and identifier).

So, in summary, we conduct 284 regressions per security, and we have a total of 5000 securities. That makes 1,420,000 regressions in total.

For some perspective, running this script takes about 50min to successfully complete.

Thank you,

• well for starters, theres a high probability your laptop has more than one core, start by making use of every core – pyCthon Aug 9 '13 at 0:06
• Two ideas; (i) don't run lm(...), use $(X'X)^{-1}X'Y$. (ii) every so often do a write.csv' or a save, and rm() to clear memory, (iii) run the as.character on the whole vector of dates instead of on a single date in each loop iteration.. – Jase Aug 9 '13 at 13:52
• Done for character cast. Thanks! However, the use of vector/matrix multiplication instead of lm() might induce more prior calculations: note that the vector Y might have NA's, while the vector X can not take NA values. That means that prior to tcomputing the OLS betas using the matrix form, we need to match the index of non-NA values of X, with the relevant values of Y, so that they are time-aligned. The match() function that would help achieve that would take time in itself.. – Mayou Aug 9 '13 at 14:43
• For the non-NA matching try data[complete.cases(data),] – Jase Aug 10 '13 at 7:43
• @pyCthon Since all the regressions are independent, it seems that what you suggested might really improve the runtime. However, I am not familiar with parallel computing nor its implementation in R. Although I am an extensive R user, my knowledge isn't at the advanced level. – Mayou Aug 13 '13 at 14:02

This isn't exactly what I would call advanced but running each regression on a separate core in a parallel foreach loop would help

http://cran.r-project.org/web/packages/foreach/foreach.pdf

• Also, in order to enable parallel computing using foreach, you need to call and register the package "doMC". The problem is that the multicore functionality currently only works with operating systems that support the fork system call (which means that Windows isn't supported). As I work on a Windows machine, this option is unfortunately impossible. Note that if the doMC package is not enabled, the foreach loop would perform the operations sequentially rather in in parallel. – Mayou Aug 14 '13 at 16:02
• I found a way around that by using "doParallel". It is a "parallel backend" for the foreach package, applicable in a Windows machine. By the way, I have now successfully reduced the sequential computations to 30min. I am sure that parallelizing would make the computations even faster. But it seems that the entire code has to be modified to accommodate the parallel set-up for the foreach %dopar% loop. – Mayou Aug 14 '13 at 16:58
• @Mariam nice, next steps are to figure out how to profile in R and find out which parts are taking the longest and focus on that – pyCthon Aug 16 '13 at 2:55

It appears that you are re-running the regression with each new data point. Instead, you should use an update/online formula (see an excellent answer by the famous Dr. Huber at stats.se).

You can find an implementation in the R package biglm. If it doesn't have all the features you need (no windowing out of old data) you can at least adapt it and use it to unit test your own work.

• +1: Seems to be an interesting package for Big Data applications with R. – vonjd Aug 11 '13 at 12:20
• The fact that biglm() is appending the datasets at each update really doesn't help.. As I need to delete the previous dataset from memory and use a new chunk of data everytime, the runtime doesn't improve much from the use of plain-vanilla lm() – Mayou Aug 12 '13 at 17:28
• I'm pretty sure biglm is not appending datasets. It sounds like you have a bug. – Brian B Aug 12 '13 at 17:54
• Well, I tried to apply biglm() using the dataset "trees" in R. library(biglm) data(trees) ff<-log(Volume)~log(Girth)+log(Height) chunk1<-trees[1:10,] chunk2<-trees[11:20,] chunk3<-trees[21:31,] a <- biglm(ff,chunk1) b <- update(a,chunk2) c <- biglm(a, chunk2) d <- biglm(a, rbind(chunk1,chunk2)). Technically, if biglm() doesn't append previous dataset, the coefficients from summary(b) and summary(c) should be the same. They are not. Notice that the coefficients from summary(b) and summary(d) are however the same, which proves that biglm() does append the dataset when update() is used.. – Mayou Aug 12 '13 at 19:25
• Could you please copy-paste the code I have just posted in my previous comment, and confirm the results that I have stated? Thanks – Mayou Aug 12 '13 at 19:29

It's really important to vectorize operations as much as possible when working with big data in R when speed is a consideration. The code below is an example of multiple regression performed on a matrix with 1000 rows and 10000 columns with the independent variables of interest in each column. The same 5 covariates are also controlled for in every model. It should take less than 10 seconds to run.

n <- 1000
m <- 10000
k <- 5
S <- matrix(2*rnorm(n*m), n, m) # matrix containing (simulated) independent variables of     interest in each column
y <- rnorm(n)
X0 <- matrix(rnorm(n*k), n, k) # # matrix containing (simulated) covariates
X <- cbind(1, X0)

U1 <- crossprod(X, y)
U2 <- solve(crossprod(X), U1)
ytr <- y - X %*% U2
U3 <- crossprod(X, S)
U4 <- solve(crossprod(X), U3)
Str <- S - X %*% U4
Str2 <- colSums(Str ^ 2)

b <- as.vector(crossprod(ytr, Str) / Str2) # Beta's for each column in S after controlling for covariates
## calculate residual error
sig <- (sum(ytr ^ 2) - b ^ 2 * Str2) / (n - k - 2)
## calculate standard error for beta
err <- sqrt(sig * (1 / Str2))
p <- 2 * pnorm(-abs(b / err))
logp <- -log10(p) # -log10 p-values for each Beta estimate


Convert the problem to a matrix format, and if possible use something like MATLAB because R is significantly slower for matrices - such as the MATLAB's index() function is super fast compared to R's match() function,

Once in matrix format, use diligent use of the expression written by Jase in the comments. There is also the fastmatch package if you want to stick to R, provided the data is sorted a-priori.

Also, a separate thought: Fama-McBeth regressions are usually run over time cross-sectional, than over securities in a time-series.