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.
 A: 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.
A: 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://www.jstor.org/pss/2627674 
http://en.wikipedia.org/wiki/Exponential_smoothing
A: 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}$.
A: 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}
A: 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.
A: 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.
A: 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..
