How can I efficiently calculate an OLS fit for N multiple variables for a rolling window?

I've worked out how to do it for 1 and 2 variable linear fits, I'd like to extend to the general case of N variables if possible (or at least to 3). I don't want to calculate the whole regression for each individual window, I want to update my calculation with each new data point.

For example, for 1 variable fit of the form: $y_i = \alpha + \beta x_i + \epsilon_i$ (where $\epsilon$ is the error term)

$$\beta = (n\sum_{i=1}^{n} x_i y_i - \sum_{i=1}^{n} x_i \sum_{i=1}^{n} y_i) / (n \sum_{i=1}^{n} x_i^2 - (\sum_{i=1}^{n} x_i)^2)$$

$$\alpha = (\sum_{i=1}^{n} y_i - \beta \sum_{i=1}^{n} x_i)/n$$

In this case I only need to maintain rolling values of $\sum x$, $\sum y$, $\sum xy$, and $\sum x^2$ and I can easily calculate the fit parameters quickly at each window step without redoing the whole thing.

I achieved the same for a 2 variable multiple fit of the form: $y_i = \alpha + \beta_{1} x_{i1} + \beta_{2} x_{i2} + \epsilon_i$

It's pretty messy: $$\sum X_1Y = \sum_{i=1}^{n} x_{i1}y_i - \frac{1}{n}(\sum_{i=1}^{n}x_{i1} \sum_{i=1}^{n} y_i)$$ $$\sum X_2Y = \sum_{i=1}^{n} x_{i2}y_i - \frac{1}{n}(\sum_{i=1}^{n}x_{i2} \sum_{i=1}^{n} y_i)$$ $$\sum X_1X_2 = \sum_{i=1}^{n} x_{i1}x_{i2} - \frac{1}{n}(\sum_{i=1}^{n}x_{i1} \sum_{i=1}^{n} x_{i2})$$ $$\sum X_1X_1 = \sum_{i=1}^{n} x_{i1}^2 - \frac{1}{n}(\sum_{i=1}^{n}x_{i1})^2$$ $$\sum X_2X_2 = \sum_{i=1}^{n} x_{i2}^2 - \frac{1}{n}(\sum_{i=1}^{n}x_{i2})^2$$

we then get $$d = (\sum X_1X_1 \sum X_2X_2) - (\sum X_1X_2)^2$$ $$\beta_1 = \frac{1}{d}(\sum X_2X_2 \sum X_1Y - \sum X_1X_2 \sum X_2Y)$$ $$\beta_2 = \frac{1}{d}(\sum X_1X_1 \sum X_2Y - \sum X_1X_2 \sum X_1Y)$$

Again, here I only need to maintain rolling sums of $\sum x_1$, $\sum x_2$, $\sum y$, $\sum x_1y$, $\sum x_2y$, $\sum x_1x_2$, $\sum x_1^2$, and $\sum x_2^2$ to calculate the regression, I don't need to do the fit individually for each window, and so it's really quick. (Note I only care about the $\beta$ values).

I now want to extend this to 3 variables (and possibly more). I've started working through the algebra and it's hairy. I feel sure I must be missing a trick or a pattern here that makes the extension to 3+ variables relatively easy.

Can anyone suggest how I can extend my technique above to 3 variables (or possibly more) or otherwise to efficiently calculate the $\beta$ for a rolling window? Thanks!

  • $\begingroup$ Have you heard of the Sherman-Morrison formula? de.wikipedia.org/wiki/Sherman-Morrison-Woodbury-Formel $\endgroup$ – Christoph Hanck Nov 25 '19 at 13:42
  • $\begingroup$ Also, such packages may be useful: cran.r-project.org/web/packages/roll/roll.pdf $\endgroup$ – Christoph Hanck Nov 25 '19 at 13:42
  • $\begingroup$ Thanks @ChristophHanck, I'd not heard of the Sherman-Morrison-Woodbury formula. I'll have to think through how I might implement that, though the mention of numerical stability issues is a bit of a concern! $\endgroup$ – Justin Nov 26 '19 at 10:04
  • $\begingroup$ Regarding existing packages, unfortunatley I'm not using R. I did look at some R source code (in another package) and now I've looked through the source code for roll but I can't work out what it does, roll_lm_z looks key, but I can't decipher how it works. Nonetheless I would like to code it myself. I guess the question boils down to how have other people done this?! $\endgroup$ – Justin Nov 26 '19 at 10:10

The formulas may be completely correct but not so easy to extend. It'd be better to consider appropriate matrix algebra. This paper (especially the first 7pages) has all the formulas necessary for an increasing-window regression via fast updates (Hostetter, Recursive Estimation). If you download the reference you'll see it's basically an application of this 'Sherman-Morrison-Woodbury' formula for updating inverse matrices.

The regression they consider is $$z = Hx + v$$ (funny notation but it's a Signal Processing publication). The updating equations are Kalman-filter-like but simpler:

The variables are $x$ is the regression coefficients and $\hat{x}$ our estimate, $h$ is the observed independent data. $z$ is the observed (scalar) dependent data and $P = (H'H)^{-1}$ is part of the covariance of the estimator $\hat{x}$.

Measurement model (for scalar measurements)

$$z_{k+1}= h^T (k + 1 ) x ( k + 1 ) + v_{k+1} $$

Predictor-corrector $$\hat{x}(k + 1) = \hat{x}(k) + K(k + 1)[z_{k+1} - h^T(k + 1)\hat{x}(k)]$$

Corrector gain $$K(k + 1) = P(k)h(k + 1)c(k + 1)$$

Gain quantities $$c(k + 1) = [h^T(k + 1)P(k)h(k + 1) + 1 ] ^{-1}$$

$$P(k + 1) = P(k) - P(k)h(k + 1)c(k + 1)h^T(k + 1)P(k)$$ or $$P(k+1) = [ I - K (k+ 1)h^T(k+ 1)]P(k)$$

Basically, it is more straightforward to deal with your updating equations when using matrix algebra

  • $\begingroup$ Thanks @nbf, I'll have to work through this. However you mention a paper but there's no link or reference? $\endgroup$ – Justin Sep 22 '20 at 15:53
  • $\begingroup$ The link to the Hostetter chapter is in my posting above. $\endgroup$ – NBF Sep 22 '20 at 15:55
  • $\begingroup$ Gene H Hostetter Recursive Estimation in Handbook of Digital Signal Processing (by Douglas Elliott) 1987 $\endgroup$ – NBF Sep 22 '20 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.