Theil’s Blus Residuals and R Tools for Testing and Removing Autocorrelation and Heteroscedasticity

March 28, 2014
By

Guest post by Hrishikesh (Rick) D. Vinod, Professor of Economics, Fordham University.

Theil (1968) proposed a transformation of regression residuals so that they are best, unbiased, linear, scalar (BLUS). No R code is available to implement them. I am providing the detailed description of the properties of BLUS residuals to the uninitiated and code. The matrix algebra itself is quite interesting and the power of R in deriving complicated matrix algebra expression involving eigenvalues and eigenvectors is well illustrated in the code. Moreover, the following paper available free at the link

Vinod, Hrishikesh D., Theil’s Blus Residuals and R Tools for Testing and Removing Autocorrelation and Heteroscedasticity (March 21, 2014). Available at SSRN: http://ssrn.com/abstract=2412740

It goes beyond testing for autocorrelation and/or heteroscedasticity with BLUS residuals or Durbin-Watson test or Breusch-Pagan test. It provides links to flexible R code for doing generalized least squares (GLS) to correct for any autocorrelation and heteroscedasticity at the same time. It is function called autohetero
Any improvements to the code to remove loop etc will be appreciated and acknowledged. The BLUS function when regressing y on a matrix of regressors called x (e.g. reg=lm(y~x)) is

# blus residuals
blus = function(y, x) {
T = length(y)
u = resid(lm(y ~ x))
ane = rep(1, T)
bigx = cbind(ane, x)
p = NCOL(bigx)
u0 = u[1:p]
# print(u0)
u1 = u[(p + 1):T]
X0 = bigx[1:p, 1:p]
X0inv = solve(X0)
# print(X0inv)
X1 = bigx[(p + 1):T, ]
xtxinv = solve(t(bigx) %*% bigx)
mtx = X0 %*% xtxinv %*% t(X0)
ei = eigen(mtx, symmetric = TRUE)
disq = ei$values
# print(c("disq=", disq),quote=FALSE)
di = sqrt(disq)
c1 = di/(1 + di) # p constants
q = ei$vectors # matrix of eigenvectors
sum1 = matrix(0, p, p) #initialize to 0 before summing
for (i in 1:p) {
qi = q[, i] # ith column is the eigenvector
mtx2 = qi %*% t(qi) #outer product
# print(mtx2)
sum1 = sum1 + (c1[i] * mtx2) #p by p matrix
# print(sum1)
} #end of the for loop
m1 = X0inv %*% sum1 # a p by p matrix
v1 = m1 %*% u0 # p by 1 vector
u1blus = u1 - X1 %*% v1 #(T-p) by 1 vector
return(u1blus)
}
# Examples blus(y,x) or b1=blus(y,cbind(x1,x2))



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

5 Responses to Theil’s Blus Residuals and R Tools for Testing and Removing Autocorrelation and Heteroscedasticity

  1. Jeffrey Evans on March 28, 2014 at 5:36 pm

    The resulting vector is not the same length as the data specified in the model (data n=155, blus n=153).
    Coincidentally, the vector is lacking 3 observations which is the same as the diagonal of the eigenvector matrix.

    I wanted to look at this in relation to the Lagrange spatial dependence test so I tested this code on the meuse dataset in the sp library.

    • Rick Vinod on March 29, 2014 at 9:32 am

      If you read the theory of BLUS residuals, you will
      see that the length of the blus residuals vector
      is T-p where p are the number of estimated coefficients
      including the intercept. My function is behaving
      the way it is supposed to!

  2. Luke Hartigan on March 28, 2014 at 9:17 pm

    Hi,

    Great post. I have a couple of suggestions that might be worth considering:

    Instead of: u = resid(lm(y ~ x))

    try: u <- qr.resid(qr(x), y))

    Instead of:
    ane = rep(1, T)
    bigx = cbind(ane, x)

    Try:
    bigx <- cbind(1, x) # use R's recycling rules

    Instead of: X0inv = solve(X0)

    Try: X0inv <- qr.solve(X0)

    Instead of: xtxinv = solve(t(bigx) %*% bigx)

    Try: xtxinv <- qr.solve(crossproduct(bigx))

    Instead of: mtx2 = qi %*% t(qi) #outer product

    Try: mtx2 <- tcrossprod(qi)

    Best regards,
    Luke

  3. Rick Vinod on March 29, 2014 at 9:36 am

    Dear Luke
    These are great suggestions. Many thanks.
    These should be useful to many R programmers in other contexts.

    R is great. One learns something new all the time.

  4. Rick Vinod on March 29, 2014 at 10:21 am

    #blus residuals revised
    blus=function(y,x){
    #author H. D. Vinod
    #suggestions by Luke Hartigan March 28,2014
    n=length(y)
    bigx=cbind(1, x)
    #u=resid(lm(y~x))
    u=qr.resid(qr(bigx), y)
    p=NCOL(bigx)
    u0=u[1:p]
    u1=u[(p+1):n]
    X0=bigx[1:p,1:p]
    X0inv=qr.solve(X0)
    X1=bigx[(p+1):n,]
    xtxinv = qr.solve(crossprod(bigx))
    mtx=X0 %*% xtxinv %*% t(X0)
    ei=eigen(mtx,symmetric=TRUE)
    disq=ei$values
    #print(c(“disq=”, disq),quote=FALSE)
    di=sqrt(disq)
    c1=di/ (1+di)# p constants
    q= ei$vectors# matrix of eigenvectors
    sum1=matrix(0,p,p)#initialize to 0 before summing
    for (i in 1:p){
    qi=q[,i]# ith column is the eigenvector
    #mtx2=qi %*% t(qi)#outer product
    mtx2=tcrossprod(qi)#outer product
    sum1=sum1+(c1[i]*mtx2)#p by p matrix
    }#end of the for loop

    m1=X0inv %*% sum1# a p by p matrix
    v1=m1 %*% u0# p by 1 vector
    u1blus=u1-X1 %*% v1 #(n-p) by 1 vector
    return(u1blus)}
    #Examples set.seed(2);y=runif(12);x=runif(12)
    #blus(y,x)
    #b1=blus(y,cbind(x1,x2))