Solved – In weighted least squares, how to weight the residuals to get an accurate “z score”

regression

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.
enter image description here

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.

Best Answer

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.

Related Question