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@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.
______________________________________________
R-help@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.
