[R] Newey West and Singular Matrix + library(sandwich)

[ivo welch]

thank you, achim.  I will try chol2inv.

sandwich is a very nice package, but let me make some short
suggestions.  I am not a good econometrician, so I do not know what
prewhitening is, and the vignette did not explain it.  "?coeftest" did
not work after I loaded the library.  automatic bandwidth selection
can be a good thing, but is not always.

as to my own little function, I like the idea of specifying my choice
of overlapping periods.  for me, the need often arises in a case of
regressions in which the variables are measured over overlapping
intervals, so I know the number of periods to worry about.  besides,
my function has an attractive simplicity to it.  it is short.

assert does indeed not exist in R, but it should be.  YMMV.

    assert <- function( cond, ... ) { if (!cond) cat(...,
file=stderr()); stopifnot(cond) }

thanks again.

/iaw

On Wed, Sep 22, 2010 at 7:41 PM, Achim Zeileis <Achim.Zeileis at uibk.ac.at> wrote:
> 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.
>