[Gmm-commits] r224 - in pkg/momentfit: R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 8 22:48:13 CET 2024
Author: chaussep
Date: 2024-02-08 22:48:13 +0100 (Thu, 08 Feb 2024)
New Revision: 224
Modified:
pkg/momentfit/R/gmmfit-methods.R
pkg/momentfit/R/momentModel-methods.R
pkg/momentfit/R/momentModel.R
pkg/momentfit/R/rModel-methods.R
pkg/momentfit/R/sgmmfit-methods.R
pkg/momentfit/R/weak.R
pkg/momentfit/man/momentStrength-methods.Rd
pkg/momentfit/man/weakTest.Rd
pkg/momentfit/vignettes/empir.bib
pkg/momentfit/vignettes/weak.Rmd
pkg/momentfit/vignettes/weak.pdf
Log:
changed a few function to complete the MOP test. Starting the Lewis and Mertens (2022) test
Modified: pkg/momentfit/R/gmmfit-methods.R
===================================================================
--- pkg/momentfit/R/gmmfit-methods.R 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/R/gmmfit-methods.R 2024-02-08 21:48:13 UTC (rev 224)
@@ -182,13 +182,8 @@
lower.tail = FALSE))
df.adj <- attr(v, "type")$df.adj
stest <- specTest(object, df.adj=df.adj)
- vcovType <- switch(object at model@vcov,
- HAC="HAC",
- iid="OLS",
- MDS="HC",
- CL="CL")
strength <- if (testStrength){
- momentStrength(object at model, coef(object), vcovType)
+ momentStrength(object at model, coef(object))
} else { list(strength=NULL, mess=NULL) }
dimnames(coef) <- list(names(par),
c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
Modified: pkg/momentfit/R/momentModel-methods.R
===================================================================
--- pkg/momentfit/R/momentModel-methods.R 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/R/momentModel-methods.R 2024-02-08 21:48:13 UTC (rev 224)
@@ -742,14 +742,15 @@
})
setMethod("momentStrength", signature("linearModel"),
- function(object, theta, vcovType=c("OLS","HC","HAC","CL")){
+ function(object, theta){
spec <- modelDims(object)
+ vcovType <- object at vcov
getF <- function(i)
{
resu <- lm(X[, i] ~ Z - 1)
- v <- switch(vcovType, OLS = vcov(resu),
- HC = vcovHC(resu, "HC1"),
- HAC = vcovHAC(resu),
+ v <- switch(vcovType, iid = vcov(resu),
+ MDS = do.call(vcovHC, c(object at vcovOptions, list(x = resu))),
+ HAC = do.call(vcovHC, c(object at vcovOptions, list(x = resu))),
CL = do.call(vcovCL, c(object at vcovOptions, list(x = resu))))
v <- v[!exoInst, !exoInst]
b <- coef(resu)[!exoInst]
@@ -759,7 +760,6 @@
}
EndoVars <- !(spec$parNames %in% spec$momNames)
exoInst <- spec$momNames %in% spec$parNames
- vcovType <- match.arg(vcovType)
if (all(!EndoVars)) {
fstats <- NULL
mess <- "No endogenous variables: no strength measure"
@@ -1154,8 +1154,12 @@
}
}
}
+ newVcov <- FALSE
if (!is.null(arg[["vcov"]]) && !object at smooth)
+ {
+ newVcov <- arg[["vcov"]] != object at vcov
object at vcov <- arg[["vcov"]]
+ }
if (object at vcov == "HAC" || object at smooth)
{
if (is.null(arg$vcovOptions))
@@ -1170,6 +1174,19 @@
if (object at smooth && !identical(arg$vcovOptions, object at vcovOptions))
chk <- TRUE
object at vcovOptions <- arg$vcovOptions
+ } else {
+ if (is.null(arg$vcovOptions))
+ arg$vcovOptions <- list()
+ if (newVcov)
+ object at vcovOptions <- list()
+ if (length(object at vcovOptions))
+ {
+ tmp <- c(arg$vcovOptions, list(object=object at vcovOptions))
+ arg$vcovOptions <- do.call(update, tmp)
+ }
+ arg$vcovOptions <- .getVcovOptions(object at vcov, NULL, arg$vcovOptions,
+ FALSE)
+ object at vcovOptions <- arg$vcovOptions
}
if (is.null(arg$survOptions))
arg$survOptions <- list()
Modified: pkg/momentfit/R/momentModel.R
===================================================================
--- pkg/momentfit/R/momentModel.R 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/R/momentModel.R 2024-02-08 21:48:13 UTC (rev 224)
@@ -56,6 +56,14 @@
}
if (option$type != "HC0")
stop("Only meatCL with type HC0 is allowed")
+ } else if (type=="MDS") {
+ option <- list(type="HC3")
+ if (length(addO) > 0)
+ {
+ if (!all(names(addO) %in% names(option)))
+ stop(paste("Wrong options for vcov of type", type))
+ option[names(addO)] <- addO
+ }
} else {
option <- list()
}
Modified: pkg/momentfit/R/rModel-methods.R
===================================================================
--- pkg/momentfit/R/rModel-methods.R 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/R/rModel-methods.R 2024-02-08 21:48:13 UTC (rev 224)
@@ -790,7 +790,7 @@
### the model has been modified.
setMethod("momentStrength", "rlinearModel",
- function(object, theta, vcovType = c("OLS", "HC", "HAC")) {
+ function(object, theta) {
fstats <- NULL
mess <- "No strength measure available for restricted models"
list(strength=fstats, mess=mess)
Modified: pkg/momentfit/R/sgmmfit-methods.R
===================================================================
--- pkg/momentfit/R/sgmmfit-methods.R 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/R/sgmmfit-methods.R 2024-02-08 21:48:13 UTC (rev 224)
@@ -199,11 +199,9 @@
names(coef) <- eqnNames
df.adj <- attr(v, "type")$df.adj
stest <- specTest(object, df.adj = df.adj)
- vcovType <- switch(object at model@vcov, HAC = "HAC", iid = "OLS",
- MDS = "HC")
strength <- lapply(1:neqn, function(i) {
if (inherits(object at model, "slinearModel") & testStrength)
- momentStrength(object at model[i], par[[i]], vcovType)
+ momentStrength(object at model[i], par[[i]])
else
list(strength=NULL, mess=NULL)})
wSpec <- object at wObj@wSpec
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/R/weak.R 2024-02-08 21:48:13 UTC (rev 224)
@@ -396,7 +396,7 @@
max(abs(ne))/BM
}
ew2 <- eigen(w$w2)$value
- maxf <- if (type == "LIML") max(ew2)/sum(ew2) else 1-2*min(ew2)/sum(ew2)
+ maxf <- if (type == "LIML") max(ew2)/sum(ew2) else abs(1-2*min(ew2)/sum(ew2))
b <- 1
i <- 1
while (TRUE)
@@ -418,9 +418,9 @@
Be/tau
}
-
MOPtest <- function(object, tau=0.10, size=0.05, print=TRUE,
- estMethod = c("TSLS", "LIML"), simplified = TRUE, ...)
+ estMethod = c("TSLS", "LIML"), simplified = TRUE,
+ digits = max(3L, getOption("digits") - 3L), ...)
{
estMethod <- match.arg(estMethod)
if (!inherits(object, "linearModel"))
@@ -443,18 +443,27 @@
if (!is.null(X1))
{
fitX1 <- lm.fit(X1, Z2)
- Z2 <- fitX1$residuals
+ Z2 <- as.matrix(fitX1$residuals)
X2 <- qr.resid(fitX1$qr, X2)
y <- qr.resid(fitX1$qr, y)
}
Z2 <- qr.Q(qr(Z2))*sqrt(nrow(Z2))
colnames(Z2) <- paste("Z", 1:ncol(Z2), sep="")
+ if ((ncol(Z2)-ncol(X2))<2)
+ if (!simplified)
+ {
+ warning(paste("The generalized test is for models with 2 ",
+ "and more over-identifying restrictions:\n",
+ " simplified is changed to TRUE", sep=""))
+ simplified <- TRUE
+ }
if (simplified)
{
if (estMethod == "LIML")
{
- warning("The simplified test is not defined for LIML")
- return(invisible())
+ warning(paste("The simplified test is not defined for LIML\n",
+ "estMethod changed to TSLS", sep=""))
+ estMethod <- "TSLS"
}
g <- reformulate(colnames(Z2), colnames(X2), FALSE)
h <- reformulate(colnames(Z2), NULL, FALSE)
@@ -475,9 +484,10 @@
vcovOptions = object at vcovOptions)
v <- cbind(y-Z2%*%b1, X2-Z2%*%b2)
omega <- crossprod(v)/nrow(v)
- w <- list(w1 = w[1:ncol(Z2), 1:ncol(Z2)],
- w2 = w[(ncol(Z2)+1):ncol(w), (ncol(Z2)+1):ncol(w)],
- w12 = w[(ncol(Z2)+1):ncol(w), 1:ncol(Z2)],
+ w <- vcov(m, list(b1,b2))
+ w <- list(w1 = w[1:ncol(Z2), 1:ncol(Z2), drop=FALSE],
+ w2 = w[(ncol(Z2)+1):ncol(w), (ncol(Z2)+1):ncol(w), drop=FALSE],
+ w12 = w[1:ncol(Z2), (ncol(Z2)+1):ncol(w), drop=FALSE],
omega = crossprod(v)/nrow(v))
x <- getMOPx(w, tau, estMethod, ...)
}
@@ -485,7 +495,6 @@
se <- sum(ev)
se2 <- sum(ev^2)
me <- max(ev)
-
## Z'Y = Pi x n, so Y'Z'ZY = sum(Pi^2)*n^2
## there for Y'Z'ZY/se/n = sim(Pi^2)*n/se
Feff <- sum(b2^2)/se*spec$n
@@ -499,14 +508,16 @@
return(c(Feff=Feff, Keff=Keff, x=x, critValue=crit, pvalue=pv))
cat("Montiel Olea and Pflueger Test for Weak Instruments\n")
cat("****************************************************\n", sep="")
- cat(ifelse(simplified, "Simplified Test", "Generalized Test\n"))
+ cat(ifelse(simplified, "Simplified Test", "Generalized Test"),
+ " for ", estMethod, "\n", sep="")
cat("Type of LS covariance matrix: ", vcov, "\n", sep="")
cat("Number of included Endogenous: ", ncol(X2), "\n", sep="")
cat("Effective degrees of freedom: ", Keff, "\n", sep="")
cat("x: ", x, "\n", sep="")
cat("Statistics: ", formatC(Feff, ...), "\n", sep="")
- cat(paste("Critical Value (size=",size,"): ", formatC(crit, ...), "\n", sep=""))
- cat(paste("P-Value: ", formatC(pv, ...), "\n\n", sep=""))
+ cat(paste("Critical Value (size=",size,"): ", formatC(crit, digits=digits),
+ "\n", sep=""))
+ cat(paste("P-Value: ", formatC(pv, digits=digits), "\n\n", sep=""))
invisible()
}
@@ -530,7 +541,7 @@
if (!is.null(X1))
{
fitX1 <- lm.fit(X1, Z2)
- Z2 <- fitX1$residuals
+ Z2 <- as.matrix(fitX1$residuals)
X2 <- qr.resid(fitX1$qr, X2)
y <- qr.resid(fitX1$qr, y)
}
@@ -547,8 +558,8 @@
v <- cbind(y-Z2%*%b1, X2-Z2%*%b2)
omega <- crossprod(v)/nrow(v)
w <- vcov(m, list(b1,b2))
- w1 <- w[1:ncol(Z2), 1:ncol(Z2)]
- w2 <- w[(ncol(Z2)+1):ncol(w), (ncol(Z2)+1):ncol(w)]
- w12 <- w[(ncol(Z2)+1):ncol(w), 1:ncol(Z2)]
+ w1 <- w[1:ncol(Z2), 1:ncol(Z2), drop=FALSE]
+ w2 <- w[(ncol(Z2)+1):ncol(w), (ncol(Z2)+1):ncol(w), drop=FALSE]
+ w12 <- w[(ncol(Z2)+1):ncol(w), 1:ncol(Z2), drop=FALSE]
list(w1=w1,w2=w2,w12=w12,omega=omega)
}
Modified: pkg/momentfit/man/momentStrength-methods.Rd
===================================================================
--- pkg/momentfit/man/momentStrength-methods.Rd 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/man/momentStrength-methods.Rd 2024-02-08 21:48:13 UTC (rev 224)
@@ -9,11 +9,10 @@
\alias{momentStrength,nonlinearModel-method}
\title{ ~~ Methods for Function \code{momentStrength} in Package \pkg{momentfit} ~~}
\description{
-It produces measures of the strength of the moment conditons.
+It produces measures of the strength of the moment conditions.
}
\usage{
-\S4method{momentStrength}{linearModel}(object, theta,
-vcovType=c("OLS","HC","HAC","CL"))
+\S4method{momentStrength}{linearModel}(object, theta)
}
\arguments{
\item{object}{An object of class \code{"linearModel"}}
@@ -20,16 +19,6 @@
\item{theta}{Coefficient vector at which the strength must be
measured. It does not impact the measure for objects of class
\code{linearModel}.}
- \item{vcovType}{Type of covariance matrix used to
- compute the F-test of the first-stage regression. For \code{HC},
- the function \code{\link{vcovHC}} is used with "HC1", and
- \code{\link{vcovHAC}} is used with the default setup is "HAC" is
- chosen. In \code{summary} for \code{gmmfit} objects, it is adjusted
- to the type of covariance that is set in the object. For type
- \code{CL}, clustered covariance matrix is computed. The options are
- the one included in the \code{vcovOptions} slot of the object (see
- \code{\link{meatCL}}). The object must have be defined with clusters
- for that to work. See \code{\link{momentModel}}.}
}
\section{Methods}{
\describe{
@@ -57,12 +46,27 @@
Not implemented yet.
}
}}
+
+\details{
+For now, the method only exists for linear models. It returns
+the F-statistics from the first stage regression. The type of covariance
+matrix used to compute the statistics depends on the specification of
+the model. If the argument \code{vcov} of the model is set to
+\code{"iid"}, a non robust estimator is used. If it is set to
+\code{"MDS"}, \code{"HAC"}, or \code{"CL"}, the appropriate robust
+estimator is used. To use a different type, use the method \code{update}
+to change the argument \code{vcov} of the model object. See the
+example below.
+}
+
\examples{
data(simData)
theta <- c(beta0=1,beta1=2)
-model1 <- momentModel(y~x1, ~z1+z2, data=simData)
+model1 <- momentModel(y~x1, ~z1+z2, data=simData, vcov="iid")
momentStrength(model1)
+## changing the type of vcov to get robust tests
+momentStrength(update(model1, vcov="MDS"))
}
\keyword{methods}
Modified: pkg/momentfit/man/weakTest.Rd
===================================================================
--- pkg/momentfit/man/weakTest.Rd 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/man/weakTest.Rd 2024-02-08 21:48:13 UTC (rev 224)
@@ -21,7 +21,8 @@
CDtest(object, print=TRUE, SWcrit=FALSE, ...)
MOPtest(object, tau=0.10, size=0.05, print=TRUE,
- estMethod = c("TSLS", "LIML"), simplified = TRUE, ...)
+ estMethod = c("TSLS", "LIML"), simplified = TRUE,
+ digits = max(3L, getOption("digits") - 3L), ...)
SYTables(object, print=TRUE, SWcrit=FALSE)
@@ -56,10 +57,14 @@
\code{FALSE}
}
+ \item{digits}{The number of digits to print when \code{print} is set
+ to \code{TRUE}.}
+
\item{SWcrit}{If true, the critical values and the test are ajusted to
Sanderson and Windmeijer (2016). See details.}
- \item{...}{Arguments passed to \code{\link{formatC}} to print the statistic.
+ \item{...}{Arguments passed to \code{\link{formatC}} to print the
+ statistic. For \code{MOPtest}, they are passed to \code{momentfit:::getMOPx}.
}
}
Modified: pkg/momentfit/vignettes/empir.bib
===================================================================
--- pkg/momentfit/vignettes/empir.bib 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/vignettes/empir.bib 2024-02-08 21:48:13 UTC (rev 224)
@@ -1,3 +1,11 @@
+ at article{lewis-mertens22,
+ AUTHOR={Lewis, D.J. and Mertens, K.},
+ TITLE={A robust test for weak instruments with Multiplle Endogenous Regressors},
+ JOURNAL={FRB of Dallas Working Paper No. 2208},
+ YEAR={2022},
+ url={https://ssrn.com/abstract=4144103}
+ }
+
@article{montiel-olea-pflueger13,
AUTHOR={Montiel Olea, J. and Pflueger, C.},
TITLE={A robust test for weak instruments},
Modified: pkg/momentfit/vignettes/weak.Rmd
===================================================================
--- pkg/momentfit/vignettes/weak.Rmd 2024-02-07 21:24:21 UTC (rev 223)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-02-08 21:48:13 UTC (rev 224)
@@ -732,6 +732,25 @@
It seems that we are rejecting the weak instrument hypothesis for TSLS
and not for LIML. This is surprising. We'll need to investigate.
+As mentioned by the authors, the efficient F is the robust F when the
+model is just identified. The model `mod4` created above is
+just-identified, but we need to change the argument `vcov` to `"MDS"`:
+
+```{r}
+mod4 <- update(mod4, vcov="MDS")
+MOPtest(mod4, estMethod="TSLS", print=FALSE)["Feff"]
+```
+
+We can compare with the first stage F computed by `momentStrength`. As
+we can see, it is the same as long as we choose the HC0 type.
+
+```{r}
+momentStrength(update(mod4, vcovOptions=list(type="HC0")))$strength
+```
+
+## Testing for weak instruments: @lewis-mertens22
+
+
## Data Generating Process (for later use)
The following function is used to generate dataset with $k$
Modified: pkg/momentfit/vignettes/weak.pdf
===================================================================
(Binary files differ)
More information about the Gmm-commits
mailing list