In weighted least squares, how do I weight the residuals to get an accurate "z score"

I am regressing spreads in yield curves in certain countries, in the chart below, the Spanish 2-5-10 spread, against the Italian 2-5-10 spread.

I want recent data to count more, so I weight the inputs using a decay weighting scheme with a 1 year halflife.

A "simple" regression line and the weighted regression line are shown.

I want to calculate the perpendicular distance of the current point (green) from the regression line, in number of standard errors. In the unweighted regression, in R, I will simply say:

> l <- lm(SP ~ IT, data = ss)
> last(l$residuals) / sd(l$residuals)
2013-10-28
-0.1817122


Which gives -.18 standard errors away from the regression line.

How do I do this same analysis for the weighted regression though? I am sure the following is incorrect:

> decay
function(len, halflife, sumone = TRUE) {
#function generates an exponentially decaying series
t <- len:1 # generate a series of numbers reverse order so biggest weights last
lambda <- log(2) / halflife #figure out the lambda for the halflife
w <- exp(-lambda * t) #create the weights series
if(sumone) w <- w / sum(w) #normalise sum to 1 if necessary
return(w)
}
> d <- decay(nrow(ss), 260)
> ld <- lm(SP ~ IT, data = ss, weights = d)
> last(ld$residuals) / sd(ld$residuals)
2013-10-28
-0.3667876


I should surely weight the residuals somehow, before doing the above, is that correct? Could I for example take the weighted standard deviation of the residuals that is:

> last(ld$residuals) / wt.sd(ld$residuals, d)
2013-10-28
-0.39717


where my wt.sd function looks like this:

> wt.sd
function (x, wt) {
return(sqrt(wt.var(x, wt)))
}

> wt.var
function (x, wt) {
s = which(is.finite(x + wt))
wt = wt[s]
x = x[s]
xbar = wt.mean(x, wt)
return(sum(wt * (x - xbar)^2) * (sum(wt)/(sum(wt)^2 - sum(wt^2))))
}


Basically, I want to know how to find the distance from the weighted regression line, in standard errors, accounting for the weights.

You can do this with the lm and associated functions, but you need to be a little careful about how you construct your weights.

Here's an example / walkthrough. Note that the weights are normalized so that the average weight = 1. I'll follow with what happens if they aren't normalized. I've deleted a lot of the less relevant printout associated with various functions.

x <- rnorm(1000)
y <- x + rnorm(1000)
wts <- rev(0.998^(0:999)) # Weights go from 0.135 to 1
wts <- wts / mean(wts)    # Now we normalize to mean 1
> summary(unwtd_lm <- lm(y~x))

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.04238    0.031    ---
x            1.03071    0.03268  31.539   <2e-16 ***
Residual standard error: 1.01 on 998 degrees of freedom

> summary(wtd_lm <- lm(y~x, weights=wts))

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.03436    0.03227   1.065    0.287
x            1.03869    0.03295  31.524   <2e-16 ***
Residual standard error: 1.02 on 998 degrees of freedom


You can see that with this much data we don't have much difference between the two estimates, but there is some.

Now for your question. It's not clear whether you want the distance in standard errors where the standard errors are for fitted values or for prediction, so I'll show both. Let us say we are doing this for the value $x = 1$ and the target value (green dot) $y = 1.1$):

> y_eval <- 1.10
> wtd_pred <- predict(wtd_lm, newdata=data.frame(x=1), se.fit=TRUE)
> # Distance relative to predictive std. error
> (y_eval-wtd_pred$fit[1]) / sqrt(wtd_pred$se.fit^2 + wtd_pred$residual.scale^2) [1] 0.02639818 > > # Distance relative to fitted std. error > (y_eval-wtd_pred$fit[1]) / wtd_pred$se.fit [1] 0.5945089  where I've deleted the warning message associated with predictive confidence intervals and weighted model fits. Now I'll show you how to do the residual variance calculation. First, if your weights aren't normalized, you will have problems: > wts <- rev(0.998^(0:999)) > summary(wtd_lm <- lm(y~x, weights=wts)) Estimate Std. Error t value Pr(>|t|) (Intercept) 0.03436 0.03227 1.065 0.287 x 1.03869 0.03295 31.524 <2e-16 *** Residual standard error: 0.6707 on 998 degrees of freedom > predict(wtd_lm, newdata=data.frame(x=1), interval="prediction") fit lwr upr 1 1.073049 -0.2461643 2.392262  Note how that residual standard error has gone way down and the prediction confidence interval has really changed, but the coefficient estimates themselves have not. This is because the calculation for the residual s.e. divides by the residual degrees of freedom (998 in this case) without regard for the scale of the weights. Here's the calculation, mostly lifted from the interior of summary.lm: w <- wtd_lm$weights
r <- wtd_lm$residuals rss <- sum(w * r^2) sqrt(rss / wtd_lm$df)
[1] 0.6707338


which you can see matches the residual s.e. in the previous printout.

Here's how you ought to do this calculation if you find yourself in a position where you need to do it by hand, so to speak:

> rss_w <- sum(w*r^2)/mean(w)
> sqrt(rss_w / wtd_lm$df) [1] 1.019937  However, normalizing the weights up front takes care of the need to divide by mean(w) and the various lm-related calculations come out correctly without any further manual intervention. • I am looking for the distance from the regression line, in number of standard errors of the residuals. Usually in a standard lm I would say lm <- (x ~ y); res <- lm$residuals; zscore <- last(res)/se(res). I think you are showing me how to do this in the weighted case above but would you mind clarifying exactly which bit of your answer does that? I think I mean standard errors for fitted values, as per your uncertainty above. Many thanks for your very complete answer. Commented Nov 7, 2013 at 10:20
• The second block of code contains what you want. wtd_pred$fit gives the point on the regression line for the$x\$ specified in the predict statement, the y_eval is the value for which you are calculating the distance, and the denominator is the standard error. The line under the Distance relative to fitted std error comment contains the actual calculation. Commented Nov 7, 2013 at 14:00

First, just on a side note - your yields (x and y) seems to be non-stationary - you see these 2 clouds on your scatter - they probably happened to be at that particular slope to each other just by chance because yields "randomly walked" in those particular directions. If they would randomly walked in other directions, you could have completely different slope (and you probably will have in future - something like 3rd cloud somewhere aside). So when you draw a regression line it may happen to be spurious.

Next, one of interpretations of weighted least squares regression is that residuals in your underlying model have different variances in different regions (heteroscedasticity) and you put more weight on where they have smaller variance to obtain consistent estimates of regression coeffs. But then computing residuals' stdev will not be interpretable as an estimate of their true variance - because they come from different distributions with different variances :)

But if you still want some measure of stdev taking into account relative weights or importance of points - then Mean Squared Weighted Deviation may be a way to go.

• Thanks Kochede - actually, while I hear you on heteroskedasticity, the issue is less about stationarity than about more recent data being more important. Relationships between spreads evolve, and I want recent relationships to count more than old ones. This may lead to heteroskedasticity, or it may not: we could have a situation where the slope of the regression line simply changes (though I doubt it will reverse). In both cases of course, the residuals will be non-normally distributed: agreed. But I want to transform them so that they are. Will look at your suggestion before accepting. Thanks Commented Nov 4, 2013 at 18:43