Beta's memory length in weighted least squares regression

You can use least squares regression to estimate the trend of a signal. If your inputs are time indices and your target is the signal itself, then the regression coefficient \(\beta\) is the gradient of the best fit line.

Unlike a simple rise-over-run (the value now minus the value from some chosen time ago, divided by the gap¹), weighted regression uses all the data, so it is far less thrown by noise.² The choice of weights determines the effective timescale over which the trend is calculated.

If I do a vanilla rise-over-run calculation I know exactly what my window size/lookback period/memory length is: if I look at the signal now and compare it to the value of the signal 1 minute ago, I know I am looking at the signal over a 1 minute history.

Ok, but if I am trying to estimate trends using this regression approach, what is my memory then?

What I want to tell you about is that the standard memory lengths for weighted linear regression are for the estimate of the mean, the smoothed value of the signal itself, and that if we are calculating a trend we need to do a slightly different calculation. The punchline is:

For exponential weighting the trend memory length is about 1.8 times longer than the mean’s memory length. (Exactly, it is \(6^{1/3} \approx 1.817\) times longer.)

For example, if you get a new datapoint every hour, \(\lambda = 0.9\) gives a memory of 19 hours for the mean — but a trend memory of about 35 hours.

Some Notation

Let’s add a bit of notation here. We are given \(n\) pairs of values \(\{(t_1, y_1), (t_2, y_2), \ldots, (t_n, y_n)\}\) and we model the relationship in our data as

where \(\epsilon_i\) is independent and identically distributed noise with variance \(\sigma^2\). In weighted linear regression we weight each datapoint by a factor \(0 < w_i \leq 1\), and we aim to find the \(\alpha\) and \(\beta\) that minimise the weighted sum of squares,

One reason for weighting is that some datapoints are viewed as being more important than others to regress accurately at — these are given higher weights.

Applying this to a time series, \(t_i\) will be the times at which we gain observations and \(y_i\) the signal. We want local, denoised estimates of the current gradient — we are assuming that the linear model holds only locally — and it is natural to enforce this locality by giving to recent observations high weight and giving to older observations a low weight.

The solution \(\beta^*\) is likely familiar, the weighted covariance of the signal \(y_i\) with the times of the measurements \(t_i\), divided by the variance of the time measurements:

where \(\bar t\) and \(\bar y\) are the weighted means of the times and the signal.

Things are made somewhat simpler if all our observations come at regular intervals. We can then directly identify \(t_i = i\), so \(i\) is simply the time index of the observation, and the most recent observation is \(y_n\). We restrict ourselves to this case.

Exponential weights

Exponential weighting is the obvious go-to: we give the most recent observation weight \(1\) and each observation further back in time gains a factor of \(\lambda\),

so the further back in time you are the greater the power \(\lambda\) is raised to. Under this rule, we pay more heed to recently gathered data. (Note that when a new observation comes in it’s not just that we add it to the analysis and do an update; we also retrospectively down-weight all previous data by a factor of \(\lambda\).)

The exponentially weighted moving average (EWMA) is \(\bar y_n = \bar y_{n-1} + (1-\lambda)\left(y_n - \bar y_{n-1}\right)\).³

Two helpful initial results: Using \(j = n - I\), the age of an observation, we can write the EWMA without recurrence as

and the slope \(\hat\beta_\lambda\), for long sequences, is

where \(\bar j = \lambda/(1-\lambda)\) is the average age of the data.

Our goal is to get the memory length of \(\hat\beta_\lambda\).

It’s a standard result that the memory length of the EWMA itself is

Before we handle the slope, we should remind ourselves where the simple EWMA’s memory length comes from. If you have a simple sliding window the memory is just the number of points in that sliding window. But if you have some procedure like an exponential weighting that touches every past point there is no clean in-or-out count of points.

But what you can do is analyse the noise you have in the data from a sliding window, calculate the noise from your estimator, and find the equivalent vanilla sliding window that gives you the noise you have in your estimator. That equivalent window size is your effective memory length. Crudely, an average of \(N\) points has variance \(\sigma^2/N\); and the variance of the weighted estimator (a function of \(\lambda\)) will match the variance of that sliding window estimator for some value of \(N\).

Doing this (the variance results are derived in the appendix; both are simple geometric sums): the exponentially weighted mean has⁴

while a sliding-window mean of \(N\) points has variance \(\sigma^2/N\). Matching the two,

and so we get that the matching \(N\) value is \(t_{\mathrm{mem}}^{\mathrm{mean}} = (1+\lambda)/(1-\lambda)\), as expected. This approach comes from Brown (1963), Chapter 8 (pp. 106–111).

You can also get this value without computing any variances at all, by matching the average age of the data under the two weighting schemes (that derivation is in the appendix).

Memory length of the trend

What is the memory length of this gradient estimator? More technically: what look-back window for an unweighted estimator has the same noise properties as the decaying estimator? We match noise: for a given \(\lambda\), at what window length \(N\) do the variances of the two \(\beta\) estimates match?⁵

Since \(\hat\beta_\lambda = \sum_j h_j y_{n-j}\), its variance is \(\sigma^2 \sum_j h_j^2\), which works out to (Appendix C derives it the long way):

while for the sliding window the textbook OLS result is \(\mathrm{Var}(\hat\beta_{\mathrm{OLS}}) = \sigma^2/\sum_i(t_i-\bar t)^2\), and for \(N\) equally spaced points \(\sum_i (t_i-\bar t)^2 = N(N^2-1)/12\) (also in the appendix), so

So working it through, by matching the variance between a sliding window and the weighted approach, we get

For large values of \(t_{\mathrm{mem}}^{\mathrm{trend}}\) we can drop the \(-t_{\mathrm{mem}}^{\mathrm{trend}}\) term and rearrange to get the memory length:

Here is the variance-matching picture, with simulations confirming the algebra:

The sliding-window slope variance (blue) falls as \(\approx 12/N^3\); the exponentially weighted slope variance (red) is a horizontal line for fixed \(\lambda\); the window length where they cross is the trend memory — \(1.817\times\) the mean’s memory (green).

Conclusion

When I first dug out these results I didn’t have a clear grasp of why the variance of the gradient is different from the variance of the value. I think there are two helpful views, both of which come from the expression for the slope as a covariance.

First, think about what goes into a vanilla rise-over-run calculation; it’s the value now minus the value from some time ago. Under the hood something like that has to be happening in the regression, we are adding some values and negating some others. And if we look at the numerator above, \(\mathrm{Cov}_w(t, y)\), we can find it.

This crosses zero at the average age of the data (the centre of the equivalent vanilla window). The negative ‘lobe’ is deepest at the mean memory length itself with weights going well beyond that. This naturally increases the memory length over the mean: we care (negatively) about these longer-ago data points.

That these weights sum to zero is also why the estimator ignores any constant offset in the signal, as a gradient should.

Secondly, \(\Beta\) is a covariance by construction. The slope estimator is the sample covariance of the signal with time, divided by the variance of time: \(\hat\beta_\lambda = \frac{\mathrm{Cov}_w(t, y)}{\mathrm{Var}_w(t)}\), while the mean is just a weighted average. That gives the two memory calculations different denominators. The variance of the mean is \(\frac{\sigma^2}{N}\), but the variance of the slope is \(\frac{\sigma^2}{N \, \mathrm{Var}(t)}\). The variance of the (uniform) ticks matters, and a window’s tick-variance grows like \(N^2\). With \(N\) points we get another factor of \(N\), so \(N^3\) overall. The cube root in the memory length comes simply from this.

The cube root comes from the variance of our timesteps, and results in a simple scaling of the result for the mean: a mean has variance \(\propto 1/N\) but a slope has variance \(\propto 1/N^3\), because \(\sum(t_i-\bar t)^2\) grows like \(N^3\). The \(6^{1/3}\) is just \((12/2)^{1/3}\) — the 12 from the OLS slope variance, the 2 from the weighted one.

References

• Brown, R. G. (1963). Smoothing, Forecasting and Prediction of Discrete Time Series. Prentice-Hall. Available at HathiTrust. The mean’s memory length and the \(2/(N+1)\) rule of thumb are Chapter 8 (pp. 106–111); the weighted slope variance \(2\alpha^3/(1+\beta)^3\,\sigma^2\) (in his notation \(\alpha = 1-\lambda\), \(\beta=\lambda\)) appears in Chapter 10 (p. 156, Table 10.5) and again via the matrix method in Chapter 16 (p. 229). Brown has all the ingredients needed but for whatever reason did not turn the slope variance into an equivalent window.
• Gardner, E. S. (1985). Exponential smoothing: the state of the art. Journal of Forecasting, 4(1), 1–28.

Appendix

A. The variance of the time indices

For \(N\) equally spaced points \(t_i = 1, \ldots, N\), using standard results for summing natural numbers and the squares of natural numbers (\(\sum i = N(N+1)/2\) and \(\sum i^2 = N(N+1)(2N+1)/6\)):

This is the \(\propto N^3\) growth that everything downstream depends on.

B. Variance of beta, unweighted sliding window

OLS over the most recent \(N\) points: \(\hat\beta = \sum_i (t_i-\bar t)(y_i - \bar y) / \sum_i (t_i-\bar t)^2\), a linear combination of the \(y_i\), so with i.i.d. noise

C. Variance of beta, exponential weights

The clean way to take the long-history limit is to re-index our previous observations when we get a new one: all previous time series points are relabelled as having come from one increment back in the past and the new observation is put at \(t=0\). (This only works because our time series observations come at fixed intervals.) Concretely, rather than having it be that \(t_i = i\), instead we have \(t_i = i - n\) — the regressor for the observation \(j\) steps ago is \(-j\), with weight \(\lambda^j\), \(j = 0, 1, 2, \ldots\)

Intuitively, by both having longer-ago datapoints downweighted and keeping recent, high-weight datapoints close to \(t=0\), all the sums that make up the estimator saturate at certain values as \(n \to \infty\), rather than growing without limit. We work in this limiting case throughout: the steady state the estimator settles into once it has seen plenty of data. (This moving-origin construction is Brown’s “discounted least squares” — and the saturation is also exactly what you want if you ever have to run this estimator online, keeping a running estimate that you update as each observation comes in.)

With design vector \(f(j) = (1, -j)^\top\), the weighted normal equations are \(F\,(\hat \alpha, \hat\beta)^\top = g\) with

The estimator is linear in the data, \((\hat \alpha, \hat\beta)^\top = F^{-1} g\), so its covariance is \(\sigma^2\, F^{-1} K F^{-1}\) where \(K = \sum_j \lambda^{2j} f(j) f(j)^\top\) (the same matrix with \(\lambda \to \lambda^2\), because the weights get squared inside the variance). Everything is built from three geometric sums and their \(\lambda^2\) versions — these are the limiting values as \(n\to\infty\); at finite \(n\) each carries correction terms of order \(\lambda^n\), which die off quickly:

Reading off the elements of \(\sigma^2 F^{-1} K F^{-1}\) (a few lines of algebra, or one line of sympy) gives

and as a sanity check, the same machinery with the constant-only model \(f(j) = 1\) gives \(\mathrm{Var}(\bar y_\lambda) = \frac{1-\lambda}{1+\lambda}\,\sigma^2\), the result for the mean from the main text. Both formulas match Monte Carlo simulation to within sampling error (the dots and dotted lines in the crossing figure above).

D. The mean’s memory length without variances: average age

Brown’s original argument (1963, p. 107) doesn’t mention variance at all. The average age of the data in an exponentially weighted mean is

while a sliding window of \(N\) points has ages \(0, 1, \ldots, N-1\), so average age \((N-1)/2\). Equating:

Same answer as the variance match — for the mean, the two notions of “equivalent window” coincide. (For the slope they would not, which is part of why the trend result is easy to miss.)

E. What if you don’t drop the \(-t_{\mathrm{mem}}^{\mathrm{trend}}\)?

The exact equivalent window solves the depressed cubic

whose real root has a closed form, \(N = \tfrac{2}{\sqrt 3}\cosh\!\left(\tfrac13 \cosh^{-1}\!\big(\tfrac{3\sqrt3}{2}c\big)\right)\), though a perturbative correction to the cube root is helpful: writing \(N = c^{1/3}(1+\delta)\) and keeping first order gives

Numbers:

So the cube-root form is already within ~1% at \(\lambda = 0.5\) (a five-step memory!) and the one-term correction makes it exact for all practical purposes: