[Gmm-commits] r221 - in pkg/momentfit: R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 6 16:02:14 CET 2024
Author: chaussep
Date: 2024-02-06 16:02:14 +0100 (Tue, 06 Feb 2024)
New Revision: 221
Modified:
pkg/momentfit/R/momentModel-methods.R
pkg/momentfit/R/weak.R
pkg/momentfit/man/model.matrix-methods.Rd
pkg/momentfit/man/weakTest.Rd
pkg/momentfit/vignettes/weak.Rmd
pkg/momentfit/vignettes/weak.pdf
Log:
working on MOP test for weak instruments
Modified: pkg/momentfit/R/momentModel-methods.R
===================================================================
--- pkg/momentfit/R/momentModel-methods.R 2024-02-02 22:39:34 UTC (rev 220)
+++ pkg/momentfit/R/momentModel-methods.R 2024-02-06 15:02:14 UTC (rev 221)
@@ -122,7 +122,8 @@
setGeneric("model.matrix")
setMethod("model.matrix", signature("linearModel"),
- function(object, type=c("regressors","instruments"))
+ function(object, type=c("regressors","instruments","excludedExo",
+ "includedExo", "includedEndo"))
{
type <- match.arg(type)
if (type == "regressors")
@@ -129,9 +130,37 @@
{
ti <- attr(object at modelF, "terms")
mat <- as.matrix(model.matrix(ti, object at modelF)[,])
- } else {
+ } else if (type == "instruments"){
ti <- attr(object at instF, "terms")
mat <- as.matrix(model.matrix(ti, object at instF)[,])
+ } else {
+ X <- model.matrix(object)
+ Z <- model.matrix(object, "instruments")
+ z1 <- colnames(Z) %in% colnames(X)
+ endo <- modelDims(object)$isEndo
+ if (type == "excludedExo")
+ {
+ if (all(!z1))
+ {
+ warning("No excluded Exogenous Variables")
+ return(NULL)
+ }
+ mat <- Z[,!z1,drop=FALSE]
+ } else if (type == "includedExo") {
+ if (all(endo))
+ {
+ warning("No included Exogenous Variables")
+ return(NULL)
+ }
+ mat <- X[,!endo,drop=FALSE]
+ } else {
+ if (all(!endo))
+ {
+ warning("No included endogenous variables")
+ return(NULL)
+ }
+ mat <- X[,endo,drop=FALSE]
+ }
}
mat
})
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-02-02 22:39:34 UTC (rev 220)
+++ pkg/momentfit/R/weak.R 2024-02-06 15:02:14 UTC (rev 221)
@@ -363,7 +363,7 @@
## Montiel Olea and Pflueger (2013)
-MOPtest <- function(object, tau=0.10, print=TRUE, ...)
+MOPtest <- function(object, tau=0.10, size=0.05, print=TRUE, ...)
{
if (!inherits(object, "linearModel"))
stop("object must be of class linearModel")
@@ -396,19 +396,22 @@
dat <- data.frame(cbind(X2, Z))
m <- momentModel(g, h, data=dat, vcov=object at vcov,
vcovOptions=object at vcovOptions)
- v <- vcov(gmmFit(m), !(vcov=object at vcov=="iid"))
+ Pi <- lm.fit(Z, X2)$coefficients
+ v <- vcov(m, Pi)
ev <- eigen(v)$val
se <- sum(ev)
se2 <- sum(ev^2)
me <- max(ev)
- Feff <- sum(crossprod(X2, Z)^2)/se
+ Feff <- sum(crossprod(X2, Z)^2)/se/spec$n
x <- 1/tau
Keff <- se^2*(1+2*x)/(se2+2*x*se*me)
+ crit <- qchisq(1-size, Keff, Keff*x)/Keff
+ pv <- 1-pchisq(Feff*Keff, Keff, Keff*x)
vcov <- object at vcov
if (vcov=="MDS")
vcov <- "HCCM"
if (!print)
- return(c(Feff=Feff, Keff=Keff, x=x))
+ 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("Type of LS covariance matrix: ", vcov, "\n", sep="")
@@ -415,6 +418,47 @@
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\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=""))
invisible()
}
+
+getMOPW <- function(object)
+{
+ spec <- modelDims(object)
+ if (sum(spec$isEndo)<1)
+ {
+ warning("The model does not contain endogenous variables")
+ return(NA)
+ }
+ if (sum(spec$isEndo)>1)
+ {
+ warning("The MOP test is defined for models with only one endogenous variable")
+ return(NA)
+ }
+ spec <- modelDims(object)
+ Z2 <- model.matrix(object, "excludedExo")
+ X1 <- model.matrix(object, "includedExo")
+ X2 <- model.matrix(object, "includedEndo")
+ y <- modelResponse(object)
+ fitX1 <- lm.fit(X1, Z2)
+ Z2 <- fitX1$residuals
+ X2 <- qr.resid(fitX1$qr, X2)
+ y <- qr.resid(fitX1$qr, y)
+ Z2 <- qr.Q(qr(Z2))
+ Z <- rbind(cbind(Z2, matrix(0, nrow(Z2), ncol(Z2))),
+ cbind(matrix(0, nrow(Z2), ncol(Z2)), Z2))
+ colnames(Z) = paste("Z", 1:ncol(Z), sep="")
+ dat <- as.data.frame(cbind(y=c(y,X2),Z))
+ g <- reformulate(colnames(Z), "y", FALSE)
+ h <- reformulate(colnames(Z), NULL, FALSE)
+ m <- momentModel(g, h, data = dat, vcov=object at vcov,
+ vcovOptions = object at vcovOptions)
+ b <- lm.fit(Z,dat$y)$coefficients
+ w <- vcov(m, b)*2
+ w11 <- w[1:ncol(Z2), 1:ncol(Z2)]
+ w22 <- w[(ncol(Z2)+1):ncol(w), (ncol(Z2)+1):ncol(w)]
+ w21 <- w[(ncol(Z2)+1):ncol(w), 1:ncol(Z2)]
+ list(w11=w11,w22=w22,w21=w21)
+}
Modified: pkg/momentfit/man/model.matrix-methods.Rd
===================================================================
--- pkg/momentfit/man/model.matrix-methods.Rd 2024-02-02 22:39:34 UTC (rev 220)
+++ pkg/momentfit/man/model.matrix-methods.Rd 2024-02-06 15:02:14 UTC (rev 221)
@@ -17,7 +17,7 @@
}
\usage{
\S4method{model.matrix}{linearModel}(object,
-type=c("regressors","instruments"))
+type=c("regressors","instruments","excludedExo", "includedExo", "includedEndo"))
\S4method{model.matrix}{rlinearModel}(object,
type=c("regressors","instruments"))
\S4method{model.matrix}{nonlinearModel}(object,
Modified: pkg/momentfit/man/weakTest.Rd
===================================================================
--- pkg/momentfit/man/weakTest.Rd 2024-02-02 22:39:34 UTC (rev 220)
+++ pkg/momentfit/man/weakTest.Rd 2024-02-06 15:02:14 UTC (rev 221)
@@ -20,7 +20,7 @@
CDtest(object, print=TRUE, SWcrit=FALSE, ...)
-MOPtest(object, tau=0.10, print=TRUE, ...)
+MOPtest(object, tau=0.10, size = 0.05, print=TRUE, ...)
SYTables(object, print=TRUE, SWcrit=FALSE)
@@ -38,6 +38,8 @@
}
\item{tau}{The desired bias.}
+
+ \item{size}{The size of the test for the critical value.}
\item{print}{
If code{TRUE}, the result is printed and nothing is returned. See
Modified: pkg/momentfit/vignettes/weak.Rmd
===================================================================
--- pkg/momentfit/vignettes/weak.Rmd 2024-02-02 22:39:34 UTC (rev 220)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-02-06 15:02:14 UTC (rev 221)
@@ -546,9 +546,9 @@
only option for this test, the procedure is:
- Replace y by $M_1y$, $X_2$ by $M_1X_2$ and $Z_2$ by $M_1Z_2$ and
- normalize the matrix $Z_2$ so that all vectors are
- orthonormal. Then, the model become $$y=X_2\beta_2 + u$$ and $$X_2 =
- Z_2\Pi_2+e\,.$$
+ normalize the matrix $Z_2$ so that all vectors are orthonormal (we
+ replace $Z$ by $Q$ from its QR decomposition). Then, the model becomes
+ $$y=X_2\beta_2 + u$$ and $$X_2 = Z_2\Pi_2+e\,.$$
- Obtain the robust covariance matrix estimate of $\hat\Pi_2$, $\hat{W}_2$.
@@ -580,8 +580,129 @@
MOPtest(mod2)
```
+The above procedure is the simplified version of the test. We start
+exploring the generalized test. First, we need an estimate of the
+matrix $W$. Given the structure of $Z$, the robust covariance matrix
+of the OLS estimators is the covariance matrix of the moment
+conditions, because when $Z'Z=I$, the OLS estimator of $y=Z\beta+u$ is
+$\hat\beta = Z'y=\beta+Z'u$. Therefore, the variance of $\hat\beta$ is
+the variance of the moment condition function $Z'u$. The reduced form
+for our model is:
+\begin{eqnarray*}
+y &=& X_2\beta_2 + u = Z_2[\Pi_2\beta_2] + v\\
+X_2 &=& Z_2 \Pi_2 + e
+\end{eqnarray*}
+
+For example, we can compute $W_2$ above
+as follows for `mod2`.
+
+- We extract $Z_2$, $X_1$, $X_2$ and $y$,
+
+```{r}
+## get Z
+spec <- modelDims(mod2)
+Z2 <- model.matrix(mod2, "excludedExo")
+X1 <- model.matrix(mod2, "includedExo")
+X2 <- model.matrix(mod2, "includedEndo")
+y <- modelResponse(mod2)
+```
+
+- We project $X_1$ off $Z_2$, $X_2$ and $y$ and normalize $Z$ using
+ its QR decomposition:
+
+```{r}
+fitX1 <- lm.fit(X1, Z2)
+Z2 <- fitX1$residuals
+X2 <- qr.resid(fitX1$qr, X2)
+y <- qr.resid(fitX1$qr, y)
+Z2 <- qr.Q(qr(Z2))
+```
+
+To compute $\hat{W}_2$, we can use the tools already included in the
+package. We just need to create a linear model. We set the argument
+`vcov` to `"MDS"` to obtain a robust to heteroskedasticity
+$\hat{W}_2$:
+
+```{r}
+colnames(Z2) = paste("Z", 1:ncol(Z2), sep="")
+dat <- as.data.frame(cbind(X2,Z2))
+g <- reformulate(colnames(Z2), colnames(X2), FALSE)
+h <- reformulate(colnames(Z2), NULL, FALSE)
+m2 <- momentModel(g,h,data=dat,vcov="MDS")
+```
+
+We need to compute the OLS estimator and then use the `vcov` method
+for `linearModel` objects to estimate the asymptotic variance of
+$Z_2'e/\sqrt{n}$:
+
+```{r}
+b <- lm.fit(Z2,X2)$coef
+w2 <- vcov(m2, b)
+```
+
+This is the only part of $W$ we need to estimate for the simplified
+version of the test. For the general test, we also need to estimate
+$W_{1}$, which is the asymptotic variance of $Z'v/\sqrt{n}$, and
+$W_{12}$ which is the asymptotic covariance between $Z_2'e/\sqrt{n}$
+and $Z_2'n/\sqrt{n}$. This can be done by writing the model as:
+
+\[
+\begin{pmatrix}
+y\\X_2
+\end{pmatrix} =
+\begin{pmatrix}
+Z_2 & 0\\
+0 & Z_2
+\end{pmatrix}\begin{pmatrix}
+\Pi_2\beta_2\\\Pi_2
+\end{pmatrix}+\begin{pmatrix}v\\ e
+\end{pmatrix}
+\]
+
+We can obtain $\hat{W}$ in one line:
+
+```{r}
+Z <- rbind(cbind(Z2, matrix(0, nrow(Z2), ncol(Z2))),
+ cbind(matrix(0, nrow(Z2), ncol(Z2)), Z2))
+colnames(Z) = paste("Z", 1:ncol(Z), sep="")
+dat <- as.data.frame(cbind(y=c(y,X2),Z))
+g <- reformulate(colnames(Z), "y", FALSE)
+h <- reformulate(colnames(Z), NULL, FALSE)
+m <- momentModel(g,h,data=dat,vcov="MDS")
+```
+
+We can then compute $\hat{W}$ the same way we computed $W_2$:
+
+
+```{r}
+b <- lm.fit(Z,dat$y)$coef
+w <- vcov(m, b)*2
+w
+```
+
+Note that we need to adjust the sample size. The way `m` is defined,
+the sample is multiplied by 2. Since we divide by twice the sample
+size to compute the estimator, we need to multiply the estimated W by
+2. We can see that $\hat{W}_2$ is the $2\times 2$ bottom left block of
+ $\hat W$:
+
+```{r}
+w2
+```
+
+This is what te function `getW` computes. For now, it is not exported,
+ so we need to run it using `momentfit:::`. The function returns
+ the different $\hat{W}$'s separately for convenience.
+
+
+```{r}
+res <- momentfit:::getMOPW(mod2)
+res
+```
+
+
## 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