Theil’s Blus Residuals and R Tools for Testing and Removing Autocorrelation and Heteroscedasticity
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=2412740It 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))