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

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.

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!

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

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.

#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))