# Exponentially weighted moving linear regression

I have a problem where I need to calculate linear regression as samples come in. Is there a formula that I can use to get the exponentially weighted moving linear regression? Not sure if that's what you would call it though.

• I say exponentially because I want to control the weight of old samples just like using weighted moving average. But I want the old samples to be less weighted, exponentially less, than newer samples Apr 24, 2011 at 18:57

Sounds like what you want to do is a two-stage model. First transform your data into exponentially smoothed form using a specified smoothing factor, and then input the transformed data into your linear regression formula.

http://en.wikipedia.org/wiki/Exponential_smoothing

• FYI, The Crane and Crotty reference (jstor.org/pss/2627674) is a valid URL, but when trying to get the PDF I get a message that "the file is damaged and could not be repaired". Apr 25, 2011 at 17:34

Sure, just add a weights= argument to lm() (in case of R):

R> x <- 1:10    ## mean of this is 5.5
R> lm(x ~ 1)    ## regression on constant computes mean

Call:
lm(formula = x ~ 1)

Coefficients:
(Intercept)
5.5

R> lm(x ~ 1, weights=0.9^(seq(10,1,by=-1)))

Call:
lm(formula = x ~ 1, weights = 0.9^(seq(10, 1, by = -1)))

Coefficients:
(Intercept)
6.35

R>


Here is give 'more recent' (i.e., higher) values more weight and the mean shifts from 5.5 to 6.35. The key, if any, is the $\lambda ^ \tau$ exponential weight I compute on the fly; you can change the weight factor to any value you choose and depending on how you order your data you can also have the exponent run the other way.

You can do the same with regression models involving whichever regressors you have.

• Never heard of R. Seems to be what I want, but I don't understand all of the syntax. Can you explain how that would work with this linear regression formula? (NΣXY - (ΣX)(ΣY)) / (NΣX^2 - (ΣX)^2). The formula works perfectly for my problem, but I'm writing a research paper and I need a more acceptable exponential moving form of this to work as samples come in. Apr 24, 2011 at 21:47
• It seems like the seq command in that syntax means it keeps track of the last 10 samples. I need something like exponentially weighted moving average that decreases all past weights approaching but never reaching 0 Apr 24, 2011 at 21:49
• I used 10 for the example; replace that with N where you determine N as the length of your dataset. Also, this is simply an application of weighted least squares which any decent regression textbook will cover in more detail. Apr 25, 2011 at 1:19

If you are looking for an equation of the form

$$y=\alpha_n + \beta_n x$$

after $n$ pieces of data have come in, and you are using an exponential factor $k \ge 1$ then you could use

$$\beta_n = \frac{\left(\sum_{i=1}^n k^i\right) \left(\sum_{i=1}^n k^i X_i Y_i\right) - \left(\sum_{i=1}^n k^i X_i\right) \left(\sum_{i=1}^n k^i Y_i\right) }{ \left(\sum_{i=1}^n k^i\right) \left(\sum_{i=1}^n k^i X_i^2\right) - \left(\sum_{i=1}^n k^i X_i \right)^2}$$

and

$$\alpha_n = \frac{\left(\sum_{i=1}^n k^i Y_i\right) - \beta_n \left(\sum_{i=1}^n k^i X_i\right)}{\sum_{i=1}^n k^i} .$$

If rounding or speed become issues, this can be recast in other forms. It may also be worth knowing that for $k>1$ you have $\sum_{i=1}^n k^i = \frac{k(k^n - 1)}{k-1}$.

• This is nice, but is there an online, recursive formulation? i.e. Can you write $a_n, b_n$ in terms of $a_{n-1}, b_{n-1}, x_n, y_n, k$? May 10, 2018 at 9:15
• @Peter: probably not easily, but instead of keeping all the various $X_i$ and $Y_i$ you can just store four or five running sums May 10, 2018 at 9:40

Yes you can. The method you are looking for is called exponentially weighted least squares method. It is a variation on the recursive least squares method: \begin{align} Θ ̂(k+1)&=Θ ̂(k)+K[z(k+1)-x^T (k+1) Θ ̂(k)] \\ K(k+1) &= D(k) x(k+1) [λ+x^T (k+1)D(k)x(k+1)]^(-1) \\ D(k+1) &=\frac 1 λ \bigg(D(k)-D(k)x(k+1)\bigg[λ+x^T (k+1)D(k)x(k+1)\bigg]^{-1} x^T (k+1)D(k)\bigg) \end{align} $0.9<λ<1$ typically.

Its a method developed to account for time varying parameters but are still in a linear format. which comes from the cost function: $$J(Θ)=1/2 ∑_(i=k-m)^k▒〖λ^(k-i) [z(i)-x^T (i)Θ]〗^2$$

Ordinary Least squares is calculated from the following for comparison:

the cost function being: $$J(Θ)=1/2 ∑_(i=i)^k▒[z(i)-x^T (i)Θ]^2$$ with \begin{align} Θ(k) &= D(k) X_k^T Z_k \\ Cov[Θ ̂(k)] &= σ^2 D(k) \\ D(k) &= [X_k^T X_k ]^{-1} \end{align}

• Welcome to the site, @MohSahx! It would be clearer if you can edit your formulae in latex, especially to revise the symbols like ▒. Jan 3, 2014 at 22:06
• @MohSahx it's clear that in your latex, you forgot to use {}. So $$x^{(t-k)}$$ will put the entire (t-k) parentheses in superscript but $x^(t-k)$ only puts the left parenthesis ( in the superscript.
– NBF
Mar 18, 2021 at 6:10

I'm not sure of the actual relationship of this to exponentially weighted moving linear regression, but a simple online formula for estimating an exponentially-weighted slope and offset is called Holt-Winters double exponential smoothing. From the Wikipedia page:

Given a time series $x_0 ... x_t$, and smoothing parameters $\alpha \in (0,1], \beta \in (0, 1]$, Initialize with:

\begin{align} s_1 &= x_1 \\ b_1 &= x_1 - x_0 \end{align}

And then for $t>1$: \begin{align} s_t &= (1-\alpha) (s_{t-1}+b_{t-1}) + \alpha x_t \\ b_t &= (1-\beta) b_{t-1} + \beta (s_t - s_{t-1}) \end{align}

Where $b_t$ is an estimated slope and $s_t$ is an estimated y-intercept at time t.

Maybe a statistically-inclined person can comment on how close this is to the solution of exponentially weighted moving linear regression.

If you form the Transfer Function Model y(t)=W(B)*X(t)+[THETA(B)/PHI(B)]*a(t) the operator [THETA(B)/PHI(B)] is the "smoothing component". For examnple if PHI(B)=1.0 and THETA(B)=1-.5B this would imply a set of weights of .5,.25,.125,... . in this way you could provide the answer to optimizing the "weighted moving linear regression" rather than assuming it's form.

First time here, first time posting, probably incorrect, but bare with me. So the classical linear regression calculation is as follows: \begin{align} y=\alpha + \beta \cdot x \\\\ \alpha=\frac{\left(\sum_{i}^N Y_i\right)\left(\sum_{i}^N X_i^2\right) - \left(\sum_{i}^N X_i\right)\left(\sum_{i}^N X_i\cdot Y_i\right)}{N\left(\sum_{i}^N X_i^2\right)-\left(\sum_{i}^N X_i\right)^2} \\ \beta =\frac{N\left(\sum_{i}^N X_i\cdot Y_i\right) - \left(\sum_{i}^N X_i\right)\left(\sum_{i}^N Y_i\right)}{N\left(\sum_{i}^N X_i^2\right)-\left(\sum_{i}^N X_i\right)^2} \end{align}

I used the assumption the sums in the aforementioned equations could be estimated with there individual EMA's.

\begin{align} \tilde x_n= \lambda \cdot x_n + (1-\lambda)\cdot \tilde x_{n-1} \rightarrow N_{ema}\cdot\tilde x_n\approx \sum_i^N x_i\\ \tilde y_n= \lambda \cdot y_n + (1-\lambda)\cdot \tilde y_{n-1} \rightarrow N_{ema}\cdot\tilde y_n \approx \sum_i^N y_i \\ \tilde {\psi}_n= \lambda \cdot x_n^2 + (1-\lambda)\cdot \tilde \psi_{n-1} \rightarrow N_{ema}\cdot\tilde x_n\approx \sum_i^N x_i^2\\ \tilde {\gamma}_n= \lambda \cdot y_n^2 + (1-\lambda)\cdot \tilde \gamma_{n-1} \rightarrow N_{ema}\cdot\tilde y_n\approx \sum_i^N y_i^2\\ \tilde {p}_n= \lambda \cdot (x_n\cdot y_n) + (1-\lambda)\cdot \tilde p_{n-1} \rightarrow N_{ema}\cdot\tilde p_n\approx \sum_i^N x_i\cdot y_i\\ \end{align} Where $$N_{ema}$$ is the approximated windowed sample size with following estimate: $$N\approx N_{ema}=2/\lambda - 1$$

Substitution into the classical equation, the $$N_{ema}$$ cancel and out and give:

\begin{align} \alpha=\frac{\tilde y_n \cdot \tilde \psi_n - \tilde x_n\tilde p_n}{\tilde \psi_n-\tilde x_n^2} \\ \beta =\frac{\tilde p_n - \tilde x_n \cdot \tilde y_n}{\tilde \psi_n-\tilde x_n^2} \end{align}

Python code comparing the two online methods is below:

import matplotlib.pyplot as plt
import numpy as np

def np_arange_fix(a, b, step):
b += (lambda x: step*max(0.1, x) if x < 0.5 else 0)((lambda n: n-int(n))((b - a)/step+1))
return np.arange(a, b, step)

if __name__ == "__main__":
fig1 = plt.figure()
m, bx, by = 1.0, 1.0, 2.0
mu, sigma = 0, 0.1

t = np_arange_fix(0, 3, 0.01)
x = bx + m * t
y = by + m * t

xr = x + np.random.normal(mu, sigma, len(t))
yr = y + np.random.normal(mu, sigma, len(t))

plt.plot(t, x)
plt.plot(t, xr)
plt.plot(t, y)
plt.plot(t, yr)

Nsma = 500
alpha = 2/(Nsma + 1)
print('Alpha= {:.5f}'.format(alpha))

x_sma, x_ema = np.zeros(len(t)), np.zeros(len(t))
y_sma, y_ema = np.zeros(len(t)), np.zeros(len(t))
xy_sma, xy_ema = np.zeros(len(t)), np.zeros(len(t))
xsq_sma, xsq_ema = np.zeros(len(t)), np.zeros(len(t))
ysq_sma, ysq_ema = np.zeros(len(t)), np.zeros(len(t))
r_sma, r_ema = np.zeros(len(t)), np.zeros(len(t))

x_sma_buff = np.zeros(Nsma)
y_sma_buff = np.zeros(Nsma)
xy_sma_buff = np.zeros(Nsma)
xsq_sma_buff = np.zeros(Nsma)
ysq_sma_buff = np.zeros(Nsma)

for i in range(len(t)):
xn = xr[i]
yn = yr[i]
xyn = xn * yn
x2n = xn * xn
y2n = yn * yn

# Buffering for sma
x_sma_buff = np.roll(x_sma_buff, -1, 0)
y_sma_buff = np.roll(y_sma_buff, -1, 0)
xy_sma_buff = np.roll(xy_sma_buff, -1, 0)
xsq_sma_buff = np.roll(xsq_sma_buff, -1, 0)
ysq_sma_buff = np.roll(ysq_sma_buff, -1, 0)

if i == 0:
x_sma_buff = np.ones(Nsma) * xn
y_sma_buff = np.ones(Nsma) * yn
xy_sma_buff = np.ones(Nsma) * xyn
xsq_sma_buff = np.ones(Nsma) * x2n
ysq_sma_buff = np.ones(Nsma) * y2n

x_sma_buff[-1] = xn
y_sma_buff[-1] = yn
xy_sma_buff[-1] = xyn
xsq_sma_buff[-1] = x2n
ysq_sma_buff[-1] = y2n

x_sma[i] = sum(x_sma_buff)/Nsma
y_sma[i] = sum(y_sma_buff)/Nsma
xy_sma[i] = sum(xy_sma_buff)/Nsma
xsq_sma[i] = sum(xsq_sma_buff)/Nsma
ysq_sma[i] = sum(ysq_sma_buff)/Nsma

x_ema[i] = alpha * xn + (1-alpha) * x_ema[i-1] if i >= 1 else xn
y_ema[i] = alpha * yn + (1-alpha) * y_ema[i-1] if i >= 1 else yn
xy_ema[i] = alpha * xyn + (1-alpha) * xy_ema[i-1] if i >= 1 else xyn
xsq_ema[i] = alpha * x2n + (1-alpha) * xsq_ema[i-1] if i >= 1 else x2n
ysq_ema[i] = alpha * y2n + (1-alpha) * ysq_ema[i-1] if i >= 1 else y2n

num = xy_sma[i] - x_sma[i] * y_sma[i]
den = np.sqrt((xsq_sma[i] - x_sma[i]*x_sma[i]) * (ysq_sma[i] - y_sma[i]*y_sma[i]))
r_sma[i] = num / den if abs(den) > 1e-4 else 0

num = xy_ema[i] - x_ema[i] * y_ema[i]
den = np.sqrt((xsq_ema[i] - x_ema[i]**2) * (ysq_ema[i] - y_ema[i]**2))
r_ema[i] = num / den if abs(den) > 1e-4 else 0

fig2 = plt.figure()
plt.subplot(4, 1, 1)
plt.plot(t, x_sma, label='x_sma')
plt.plot(t, x_ema, label='x_ema')
plt.legend()

plt.subplot(4, 1, 2)
plt.plot(t, y_sma, label='y_sma')
plt.plot(t, y_ema, label='y_ema')
plt.legend()

plt.subplot(4, 1, 3)
plt.plot(t, xy_sma, label='xy_sma')
plt.plot(t, xy_ema, label='xy_ema')
plt.legend()

plt.subplot(4, 1, 4)
plt.plot(t, xsq_sma, label='xsq_sma')
plt.plot(t, xsq_ema, label='xsq_ema')
plt.legend()

m_sma, b_sma = [], []
m_ema, b_ema = [], []

for i in range(len(t)):
sum_x = x_sma[i] * Nsma
sum_y = y_sma[i] * Nsma
sum_xy = xy_sma[i] * Nsma
sum_x2 = xsq_sma[i] * Nsma

den = Nsma * sum_x2 - sum_x * sum_x
m_num = Nsma * sum_xy - sum_x * sum_y
b_num = sum_y * sum_x2 - sum_x * sum_xy

if abs(den) > 1e-4:
m_sma.append(m_num/den)
b_sma.append(b_num/den)
else:
m_sma.append(0)
b_sma.append(0)

xn = x_ema[i]
yn = y_ema[i]
xyn = xy_ema[i]
x2n = xsq_ema[i]

den = x2n - xn*xn
m_num = xyn - xn * yn
b_num = yn * x2n - xn * xyn

if abs(den) > 1e-4:
m_ema.append(m_num/den)
b_ema.append(b_num/den)
else:
m_ema.append(0)
b_ema.append(0)

fig3 = plt.figure()
ax1 = plt.subplot(2, 2, 1)
plt.plot(t, b_sma, label='b_sma')
plt.plot(t, b_ema, label='b_ema')
plt.legend()

ax2 = plt.subplot(2, 2, 3, sharex=ax1)
plt.plot(t, m_sma, label='m_sma')
plt.plot(t, m_ema, label='m_ema')
plt.legend()

ax3 = plt.subplot(2, 2, 2)
plt.plot(t,  x, label='x={:.1f}*t+{:.1f}'.format(m, bx))
plt.plot(t, xr, label='x={:.1f}*t+{:.1f}+e'.format(m, bx))
plt.plot(t,  y, label='y={:.1f}*t+{:.1f}'.format(m, by))
plt.plot(t, yr, label='y={:.1f}*t+{:.1f}+e'.format(m, by))
plt.legend()

ax4 = plt.subplot(2, 2, 4, sharex=ax2)
plt.plot(t, r_sma, label='r_sma')
plt.plot(t, r_ema, label='r_ema')
plt.legend()
# fig4 = plt.figure()
# plt.plot(xr, yr)

plt.show()


Would be interesting to compare @Peter answer using double exponential smoothing with this as it avoids more calculations and less variables..