6
$\begingroup$

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.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ 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 $\endgroup$
    – brandon
    Apr 24, 2011 at 18:57

7 Answers 7

Reset to default
6
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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". $\endgroup$
    – zbicyclist
    Apr 25, 2011 at 17:34
7
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – brandon
    Apr 24, 2011 at 21:47
  • $\begingroup$ 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 $\endgroup$
    – brandon
    Apr 24, 2011 at 21:49
  • 2
    $\begingroup$ 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. $\endgroup$
    – Dirk Eddelbuettel
    Apr 25, 2011 at 1:19
4
$\begingroup$

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}$.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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$? $\endgroup$
    – Peter
    May 10, 2018 at 9:15
  • 1
    $\begingroup$ @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 $\endgroup$
    – Henry
    May 10, 2018 at 9:40
4
$\begingroup$

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}

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ Welcome to the site, @MohSahx! It would be clearer if you can edit your formulae in latex, especially to revise the symbols like ▒. $\endgroup$
    – Randel
    Jan 3, 2014 at 22:06
  • $\begingroup$ @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. $\endgroup$
    – NBF
    Mar 18, 2021 at 6:10
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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..

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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