Q:

Running many multiple regressions at once in R

I want to run time series regressions (Fama-French three factor with my new factor). I have following tables.
Table01

    Date     Port_01  Port_02 --------- Port_18
    01/1965     0.85    0.97               1.86
    02/1965     8.96    7.2                0.98
    03/1965     8.98    7.9                8.86 

Table 02

    Date        Market   SMB    HML     WXO
    01/1965      0.85    0.97    0.86    0.87
    02/1965      8.96    7.2     0.98    0.79
    03/1965      8.98     7.9    8.86    0.86

Table 01 is my Y variable and Table 02 is my independent variables. I want to regress 18 portfolios on factors on Table 02. I want the intercepts to be stored in a vector (N x 1). and residuals to be stored in a matrix (T x N), check the image. I know how to run regression individually. But It there must be a better way to do. Please help me. Thanks in advance. 

A:

In short, yes, there's an easy way to do the calculations for many regressions with the same X-matrix but different $y$'s, especially if you have software that gives you access to the calculations it does.
More specifically, typically regressions are done these days using QR decomposition. Similar manipulations are available for other methods (such as a regression that was done via a Choleski decomposition on an $X^\top X$ matrix, for example, or for other decompositions of the $X$ matrix).
With the least squares normal equations you seek to solve $X^\top X\, \hat\beta=X^\top y$
A $QR$ decomposition would write $X=QR$ where $Q$ is an orthogonal matrix and R is upper triangular.
So this gives $R^\top Q^\top Q R \hat\beta = R^\top Q^\top y$
But $Q^\top Q=I$, so this reduces to solving
$R^\top R \hat\beta = R^\top Q^\top y$
Premultiplying by $(R^\top)^{-1}$, this can be further reduced to solving $R \hat\beta = Q^\top y$
Now $R$ and $Q$ are the same for every regression, it's only $y$ that changes, and because $R$ is upper triangular, the sequence of back-substitutions to solve for it is very fast and simple.
The $QR$ decomposition itself is relatively expensive but you only need to do it once.
Other regression quantities may be similarly obtained. For example, the variance-covariance matrix of parameter estimates is (apart from a scaling factor of $\hat{\sigma}^2$) the same for each regression.
Many packages will happily give you the QR decomposition and some will give it as a by product of a regression. (In R, the qr function performs QR decomposition and lm returns the QR decomposition of the X matrix in the regression in a compressed form -- the Q matrix is not explicitly given, but exists in an implicit form. See the help on both for more information.)
I wouldn't generally suggest explicitly forming $R^{-1}Q$ in general; solving the system once you have $Qy$ is relatively pretty fast, if you do it sensibly (i.e. using a solver that takes advantage of the triangular matrix -- rather than trying to form an inverse of $R$ for example). 

