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