[R] Newey West and Singular Matrix

[ivo welch]

dear R experts:  I am writing my own little newey-west standard error
function, with heteroskedasticity and arbitrary x period
autocorrelation corrections.  including my function in this post here
may help others searching for something similar.  it is working quite
well, except on occasion, it complains that

  Error in solve.default(crossprod(x.na.omitted, x.na.omitted)) :
    system is computationally singular: reciprocal condition number =
3.63797e-23

I know that lm can do the inversion, so I presume that there is a more
stable way than qr.solve .  I looked into lm, then into lm.fit, and it
seems to invoke dqrls .  is this the recommended way, or is there a
higher-level more stable matrix inversion routine that I could use?

help is, as always, appreciated.  (also, if you see something else
silly in my code, let me know, please.)

regards,

/iaw


se.neweywest <- function( lmobject.withxtrue, ar.terms =0 ) {
  assert( (class(lmobject.withxtrue)=="lm"),
         "[se.white] works only on 'lm' objects, not on ",
class(lmobject.withxtrue), "objects \n" )
  x.na.omitted <- lmobject.withxtrue$x
  assert( class(x.na.omitted)=="matrix", "[se.white] internal
error---have no X matrix.  did you invoke with 'x=TRUE'?\n")
  r.na.omitted <- residuals( lmobject.withxtrue )

  diagband.matrix <- function( m, ar.terms ) {
    nomit <- m-ar.terms-1
    mm <- matrix( TRUE, nrow=m, ncol=m )
    mm[1:nomit,(ncol(mm)-nomit+1):ncol(mm)] <- (lower.tri(
matrix(TRUE, nrow=nomit, ncol=nomit) ))
    mm[(ncol(mm)-nomit+1):ncol(mm),1:nomit] <- (upper.tri(
matrix(TRUE, nrow=nomit, ncol=nomit) ))
    mm
  }

  ##    V(b) = inv(X'X) X' diag(e^2) X inv(X'X)
  invx <- qr.solve( crossprod( x.na.omitted, x.na.omitted ) )
  if (!ar.terms) resid.matrix <- diag( r.na.omitted^2 ) else {
    full <- r.na.omitted %*% t(r.na.omitted)

    ## the following is not particularly good.  see, we could zero out also
    ## items which are multiplications with a missing residual for example,
    ## if we do an ar1 correction, and obs 5 is missing, then the AR term on
    ## 4 and 6 could be set to 0.  right now, we just adjust for an add'l
    ## term.

    maskmatrix <- diagband.matrix( length(r.na.omitted), ar.terms )
    resid.matrix <- full * maskmatrix
  }

  invx.x <- invx %*% t(x.na.omitted)
  vmat <-  invx.x %*% resid.matrix %*% t(invx.x)
  sqrt(diag(vmat))
}