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