[R] Newey West and Singular Matrix

[Achim Zeileis]

On Wed, 22 Sep 2010, ivo welch wrote:

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

I typically leverage chol2inv(). In "strucchange" I've written a small 
convenience function solveCrossprod() that provides two different 
approaches (plus the naive method). The man page has a few comments. The 
function defintions are all one-liners, though.

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

1. There is no assert() function in base R.
2. The error message of se.neweywest() refers to se.white.
3. A much more flexible and powerful solution of this is included
    in package "sandwich", see the vignette for details. The code
      sqrt(diag(NeweyWest(lm_object, lag = 0, prewhite = 0)))
    replicates se.neweywest(lm_object) but has the following advantages:
    it also does automatic bandwidth selection, it does not require setting
    "x = TRUE", it incorporates other kernel weighting functions, supports
    prewhitening etc.

Best,
Z

> 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))
> }
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>