R help - Sep 2010 - Newey West and Singular Matrix + library(sandwich)

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

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