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

[Achim Zeileis]

On Thu, 23 Sep 2010, ivo welch wrote:

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

In general it means, that you do not look at a (typically multivariate) 
series Y, but at the residuals of a VAR(p) model for Y. In case of 
HAC covariances, people do not use the empirical estimating functions 
directly, but the residuals of a VAR(1) model.

> and the vignette did not explain it.  "?coeftest" did
> not work after I loaded the library.

Presumably because you did load "sandwich" but did not load the "lmtest" 
package in which coeftest() is provided.

> automatic bandwidth selection can be a good thing, but is not always.

As my short command in the previous mail showed, it can be conveniently 
set to any other value, just by specifying the "lag" argument.

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

But the code would be even shorter without it ;-) Just kidding...

hth,
Z

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