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

Question

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,

Answer 1

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

Answer 2

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.

Answer 3

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

Answer 4

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.