[Gmm-commits] r223 - in pkg/momentfit: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 7 22:24:21 CET 2024
Author: chaussep
Date: 2024-02-07 22:24:21 +0100 (Wed, 07 Feb 2024)
New Revision: 223
Modified:
pkg/momentfit/NAMESPACE
pkg/momentfit/R/weak.R
pkg/momentfit/man/weakTest.Rd
pkg/momentfit/vignettes/weak.Rmd
pkg/momentfit/vignettes/weak.pdf
Log:
working on MOP test... added the generalized and simplified tests
Modified: pkg/momentfit/NAMESPACE
===================================================================
--- pkg/momentfit/NAMESPACE 2024-02-07 17:01:28 UTC (rev 222)
+++ pkg/momentfit/NAMESPACE 2024-02-07 21:24:21 UTC (rev 223)
@@ -16,7 +16,7 @@
"D", "numericDeriv", "sd", "optim", "lm", "pf", "coef", "update",
"fitted", "lm.fit", "pchisq", "pnorm", "printCoefmat", "anova",
"model.frame", "reformulate", "formula", "nlminb", "kernapply",
- "constrOptim", "kernel", "confint", "qnorm", "uniroot", "getCall", "qchisq")
+ "constrOptim", "kernel", "confint", "qnorm", "uniroot", "getCall", "qchisq", "optimize")
importFrom("sandwich", "vcovHAC", "estfun","kernHAC","vcovCL", "meatCL",
"bread","bwAndrews","bwNeweyWest","weightsAndrews",
"weightsLumley", "vcovHC")
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-02-07 17:01:28 UTC (rev 222)
+++ pkg/momentfit/R/weak.R 2024-02-07 21:24:21 UTC (rev 223)
@@ -363,8 +363,66 @@
## Montiel Olea and Pflueger (2013)
-MOPtest <- function(object, tau=0.10, size=0.05, print=TRUE, ...)
+# computing x for the generalized test
+
+getMOPx <- function(w, tau, type = c("TSLS", "LIML"), e=0.0001, nP = 10000,
+ maxi = 1000)
{
+ type <- match.arg(type)
+ fb <- function(beta, w, type)
+ {
+ s1 <- w$w1-2*beta*w$w12+beta^2*w$w2
+ sig1 <- w$omega[1,1]-2*beta*w$omega[1,2]+beta^2*w$omega[2,2]
+ s2 <- w$w2
+ sig2 <- w$omega[2,2]
+ s12 <- w$w12-beta*w$w2
+ sig12 <- w$omega[1,2]-beta*w$omega[2,2]
+ ts1 <- sum(diag(s1))
+ ts2 <- sum(diag(s2))
+ ts12 <- sum(diag(s12))
+ if (type == "LIML")
+ {
+ tmp <- 2*s12-sig12*s1/sig1
+ tmp <- 0.5*tmp + 0.5*t(tmp)
+ lam <- eigen(tmp)$value[c(1,nrow(tmp))]
+ num <- ts12 - ts1*sig12/sig1 - lam
+ ne <- num/ts2
+ } else {
+ lam <- eigen(0.5*s12+0.5*t(s12))$value[c(1,nrow(s12))]
+ num <- 1-2*lam/ts12
+ ne <- num*ts12/ts2
+ }
+ BM <- sqrt(ts1/ts2)
+ max(abs(ne))/BM
+ }
+ ew2 <- eigen(w$w2)$value
+ maxf <- if (type == "LIML") max(ew2)/sum(ew2) else 1-2*min(ew2)/sum(ew2)
+ b <- 1
+ i <- 1
+ while (TRUE)
+ {
+ crit <- min(abs(fb(b, w, type)/maxf - 1),
+ abs(fb(-b, w, type)/maxf - 1))
+ if (crit <= e)
+ break
+ i <- i+1
+ b <- b*2
+ if (i>maxi)
+ {
+ warning("max iteration to find Brange reached")
+ break
+ }
+ }
+ res <- optimize(fb, c(-b,b), w=w, type=type, maximum=TRUE)
+ Be <- res$objective
+ Be/tau
+}
+
+
+MOPtest <- function(object, tau=0.10, size=0.05, print=TRUE,
+ estMethod = c("TSLS", "LIML"), simplified = TRUE, ...)
+{
+ estMethod <- match.arg(estMethod)
if (!inherits(object, "linearModel"))
stop("object must be of class linearModel")
spec <- modelDims(object)
@@ -391,21 +449,46 @@
}
Z2 <- qr.Q(qr(Z2))*sqrt(nrow(Z2))
colnames(Z2) <- paste("Z", 1:ncol(Z2), sep="")
- g <- reformulate(colnames(Z2), colnames(X2), FALSE)
- h <- reformulate(colnames(Z2), NULL, FALSE)
- dat <- data.frame(cbind(X2, Z2))
- m <- momentModel(g, h, data=dat, vcov=object at vcov,
- vcovOptions=object at vcovOptions)
- Pi <- c(crossprod(Z2,X2))/nrow(Z2)
- v <- vcov(m, Pi)
- ev <- eigen(v)$val
+ if (simplified)
+ {
+ if (estMethod == "LIML")
+ {
+ warning("The simplified test is not defined for LIML")
+ return(invisible())
+ }
+ g <- reformulate(colnames(Z2), colnames(X2), FALSE)
+ h <- reformulate(colnames(Z2), NULL, FALSE)
+ dat <- data.frame(cbind(X2, Z2))
+ m <- momentModel(g, h, data=dat, vcov=object at vcov,
+ vcovOptions=object at vcovOptions)
+ b2 <- c(crossprod(Z2,X2))/nrow(Z2)
+ w <- list(w2=vcov(m, b2))
+ x <- 1/tau
+ } else {
+ b1 <- crossprod(Z2,y)/spec$n
+ b2 <- crossprod(Z2,X2)/spec$n
+ g <- list(reformulate(colnames(Z2), "y", FALSE),
+ reformulate(colnames(Z2), colnames(X2), FALSE))
+ h <- reformulate(colnames(Z2), NULL, FALSE)
+ dat <- as.data.frame(cbind(y=y,X2,Z2))
+ m <- sysMomentModel(g=g, list(h), data = dat, vcov=object at vcov,
+ 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)],
+ omega = crossprod(v)/nrow(v))
+ x <- getMOPx(w, tau, estMethod, ...)
+ }
+ ev <- eigen(w$w2)$val
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(Pi^2)/se*spec$n
- x <- 1/tau
+ Feff <- sum(b2^2)/se*spec$n
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)
@@ -416,6 +499,7 @@
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("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="")
@@ -426,7 +510,7 @@
invisible()
}
-getMOPW <- function(object)
+getMOPw <- function(object)
{
spec <- modelDims(object)
if (sum(spec$isEndo)<1)
Modified: pkg/momentfit/man/weakTest.Rd
===================================================================
--- pkg/momentfit/man/weakTest.Rd 2024-02-07 17:01:28 UTC (rev 222)
+++ pkg/momentfit/man/weakTest.Rd 2024-02-07 21:24:21 UTC (rev 223)
@@ -20,7 +20,8 @@
CDtest(object, print=TRUE, SWcrit=FALSE, ...)
-MOPtest(object, tau=0.10, size = 0.05, print=TRUE, ...)
+MOPtest(object, tau=0.10, size=0.05, print=TRUE,
+ estMethod = c("TSLS", "LIML"), simplified = TRUE, ...)
SYTables(object, print=TRUE, SWcrit=FALSE)
@@ -40,6 +41,14 @@
\item{tau}{The desired bias.}
\item{size}{The size of the test for the critical value.}
+
+ \item{estMethod}{Which method we want to test the weak instruments
+ for? The method determined the critical value associated with a given
+ relative bias.}
+
+\item{simplified}{Should we perform the simplified test? The test
+ produce more conservative critical values. The generalized test is
+ computationally more demanding.}
\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-07 17:01:28 UTC (rev 222)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-02-07 21:24:21 UTC (rev 223)
@@ -599,7 +599,6 @@
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$,
@@ -655,13 +654,12 @@
and $Z_2'v/\sqrt{n}$. This can be done by writing the model as a
system of equations with the same regressors and instruments. The
above `g` is the second equation, so we need to add the first in a
-list and repeat `h` twice in the list of instruments:
+list and put `h` in a list:
```{r}
dat <- as.data.frame(cbind(y=y,X2,Z2))
-g <- list(reformulate(colnames(Z2), "y", FALSE),
- g)
-h <- list(h,h)
+g <- list(reformulate(colnames(Z2), "y", FALSE), g)
+h <- list(h)
m <- sysMomentModel(g,h,data=dat,vcov="MDS")
b <- list(c(crossprod(Z2,y)/nrow(Z2)),
c(crossprod(Z2,X2)/nrow(Z2)))
@@ -685,7 +683,7 @@
see $\hat{W}_2$:
```{r}
-res <- momentfit:::getMOPW(mod2)
+res <- momentfit:::getMOPw(mod2)
res$w2
```
@@ -698,43 +696,42 @@
K_{eff} = \frac{[\tr(\hat W_2)]^2(1+2x)}{\tr(\hat{W}_2'\hat{W}_2)\max\eval(\hat W_2)}
\]
-In the simplified version, the parameter $x$ is simply set to
-$1/\tau$, where $\tau$ is the worst bias relative to the
-benchmark. For the generalized test, $x$ depends on
-$B_e(\hat{W},\hat{\omega})$ for $e$ = LIML or TSLS, which needs to be
-computed using a one dimentional optimizer.
+where $x=B_e(\hat{W},\hat{\Omega})/\tau$, where $\tau$ is the worst
+bias relative to the benchmark and $B_e(\hat{W},\hat{\omega})$ is some
+function. In the simplified version, the parameter $x$ is equal to
+$1/\tau$, so only $\hat W_2$ is needed. For the generalized test, $x$
+depends on $B_e(\hat{W},\hat{\omega})$ for $e$ = LIML or TSLS, which
+needs to be computed using a one dimentional optimizer. The function
+`getMOPx` returns the value of `x`. The main arguments are `w`, which
+is the output from `getMOPw`, `tau` that we explained above and the
+type of estimation we want to test the instruments for. As usual, LIML
+is less biased so it requires a smaller $F_{eff}$ to get the same
+relative bias as TSLS.
```{r}
-f <- function(beta, w, type)
-{
- s1 <- w$w1-2*beta*w$w12+beta^2*w$w2
- s2 <- w$w2
- s12 <- w$w12-beta*w$w2
- ts1 <- sum(diag(s1))
- ts2 <- sum(diag(s2))
- ts12 <- sum(diag(s12))
- omega <- w$omega
- if (type == "LIML")
- {
- lam <- eigen(2*s12-omega[1,2]*s1/omega[1,1])$value
- num <- ts12 - ts1*omega[1,2]/omega[1,1] - lam
- ne <- num/ts2
- } else {
- lam <- eigen(s12)$value
- num <- 1-2*lam/ts12
- ne <- num*ts12/ts2
- }
- ne <- max(abs(ne))
- BM <- sqrt(ts1/ts2)
- ne/BM
-}
-w <- momentfit:::getMOPW(mod2)
-f(1,w,"LIML")
+momentfit:::getMOPx(w=res, tau=0.10, type="LIML")
```
+By default, the `MOPtest` function computes the simplified test, which
+is the one obtained above. For the generalized test, we set the
+argument `simplified` to `FALSE`.
+```{r}
+MOPtest(mod2, simplified=FALSE, estMethod="LIML")
+```
+The test is the same as above, but the critical value is smaller,
+which is expected since the simplify test tends to have higher
+critical value, especially when the number of excluded instruments is
+small. We can compare the test with the one for TSLS:
+```{r}
+MOPtest(mod2, simplified=FALSE, estMethod="TSLS")
+```
+
+It seems that we are rejecting the weak instrument hypothesis for TSLS
+and not for LIML. This is surprising. We'll need to investigate.
+
## 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