[Gmm-commits] r226 - in pkg/momentfit: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 16 23:44:27 CET 2024
Author: chaussep
Date: 2024-02-16 23:44:26 +0100 (Fri, 16 Feb 2024)
New Revision: 226
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 the Lewis-Mertens (2022) test for weak instruments
Modified: pkg/momentfit/NAMESPACE
===================================================================
--- pkg/momentfit/NAMESPACE 2024-02-09 22:15:30 UTC (rev 225)
+++ pkg/momentfit/NAMESPACE 2024-02-16 22:44:26 UTC (rev 226)
@@ -44,7 +44,7 @@
rhoET, rhoEL, rhoEEL, rhoHD, Wu_lam, EEL_lam, REEL_lam, getLambda,
solveGel, rhoETEL, rhoETHD, ETXX_lam, gelFit, evalGel, getImpProb,
evalGelObj, momFct, gel4, setCoef, lse, getK, kclassfit, CDtest,
- SYTables, SWtest, MOPtest)
+ SYTables, SWtest, MOPtest, LewMertest)
### S3 methods ###
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-02-09 22:15:30 UTC (rev 225)
+++ pkg/momentfit/R/weak.R 2024-02-16 22:44:26 UTC (rev 226)
@@ -565,3 +565,127 @@
}
+## Lewis and Mertens (2022)
+
+phiMat <- function(w2, k2, l2)
+{
+ phi <- matrix(NA, k2,k2)
+ if (length(w2) == 1)
+ return(matrix(w2,1,1))
+ for (i in 1:k2)
+ {
+ c0 <- 1+(i-1)*l2
+ c1 <- ncol(w2)
+ di <- diag(w2[,c0:c1])
+ tmp <- colSums(matrix(di, nrow=l2))
+ if (i == 1)
+ {
+ diag(phi) <- tmp
+ } else if (i==k2) {
+ phi[1,k2] <- phi[k2,1] <- tmp
+ } else {
+ j <- seq_len(length(tmp))
+ j <- cbind(j,j)
+ phi[-(1:(i-1)),][j] <- phi[,-(1:(i-1))][j] <- tmp
+ }
+ }
+ phi
+}
+
+LewMertest <- function(object, tau=0.10, size=0.05, print=TRUE,
+ estMethod = c("TSLS", "LIML"), simplified = TRUE,
+ digits = max(3L, getOption("digits") - 3L), ...)
+{
+ estMethod <- match.arg(estMethod)
+ if (!inherits(object, "linearModel"))
+ stop("object must be of class linearModel")
+ spec <- modelDims(object)
+ if (sum(spec$isEndo)<1)
+ {
+ warning("The model does not contain endogenous variables")
+ return(NA)
+ }
+ Z2 <- model.matrix(object, "excludedExo")
+ X1 <- model.matrix(object, "includedExo")
+ X2 <- model.matrix(object, "includedEndo")
+ y <- modelResponse(object)
+ if (!is.null(X1))
+ {
+ fitX1 <- lm.fit(X1, Z2)
+ 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="")
+ colnames(X2) <- paste("endo",1:ncol(X2), sep="")
+ b1 <- c(crossprod(Z2,y)/spec$n)
+ Pi2 <- crossprod(Z2,X2)/spec$n
+ g <- c(reformulate(colnames(Z2), "y", FALSE),
+ lapply(colnames(X2), function(ei)
+ reformulate(colnames(Z2), ei, 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%*%Pi2)
+ omega <- crossprod(v)/nrow(v)
+ w <- vcov(m, c(list(b1),lapply(1:ncol(Pi2), function(i) Pi2[,i])))
+ w2 <- w[(ncol(Z2)+1):ncol(w),(ncol(Z2)+1):ncol(w)]
+ phi <- phiMat(w2, ncol(X2), ncol(Z2))
+ if (dim(phi)[1] == 1)
+ {
+ gmin <- nrow(Z2)*c(crossprod(Pi2))/c(phi)
+ sqPhi <- 1/sqrt(phi)
+ } else {
+ svPhi <- svd(phi)
+ sqPhi <- svPhi$u%*%diag(1/sqrt(svPhi$d))%*%t(svPhi$v)
+ gmin <- c(nrow(Z2)*min(eigen(sqPhi%*%crossprod(Pi2)%*%sqPhi)$value))
+ }
+ lam <- 1/tau
+ cumul <- getCumul(w2, lam, sqPhi, ncol(Z2))
+ list(gmin=gmin, cumul=cumul,
+ crit=crit(size, cumul$k1, cumul$k2, cumul$k3))
+}
+
+getCumul <- function(w2, lam, phi, l2)
+{
+ if (length(phi)==1)
+ {
+ sig <- w2/c(phi)*l2
+ } else {
+ tmp <- kronecker(phi, diag(l2))
+ sig <- l2*(tmp%*%w2%*%tmp)
+ svSig <- svd(sig)
+ sig2 <- svSig$u%*%diag(svSig$d^2)%*%t(svSig$v)
+ sig3 <- svSig$u%*%diag(svSig$d^3)%*%t(svSig$v)
+ }
+ k1 <- l2*(1+lam)
+ if (length(sig) == 1)
+ {
+ k2 <- 2*(sig^2+2*lam*l2*sig)
+ k3 <- 8*(sig^3+3*lam*l2*sig^2)
+ } else {
+ svSig <- svd(sig)
+ sig2 <- svSig$u%*%diag(svSig$d^2)%*%t(svSig$v)
+ sig3 <- svSig$u%*%diag(svSig$d^3)%*%t(svSig$v)
+ tmp <- phiMat(sig2, dim(phi)[1], l2)
+ k2 <- 2*(max(eigen(tmp)$val)+2*lam*l2*max(svSig$d))
+ tmp <- phiMat(sig3, dim(phi)[1], l2)
+ k3 <- 8*(max(eigen(tmp)$val)+3*lam*l2*max(svSig$d)^2)
+ }
+ w <- k2/k3
+ v <- 8*k2*w^2
+ list(k1=k1, k2=k2, k3=k3, w=w, v=v)
+}
+
+crit <- function(alpha, k1, k2, k3)
+{
+ w <- k2/k3
+ v <- 8*k2*w^2
+ cr <- qchisq(1-alpha, v)
+ (cr-v)/(4*w)+k1
+}
+
+
+
Modified: pkg/momentfit/man/weakTest.Rd
===================================================================
--- pkg/momentfit/man/weakTest.Rd 2024-02-09 22:15:30 UTC (rev 225)
+++ pkg/momentfit/man/weakTest.Rd 2024-02-16 22:44:26 UTC (rev 226)
@@ -3,6 +3,7 @@
\alias{SWtest}
\alias{MOPtest}
\alias{SYTables}
+\alias{LewMertest}
\title{
Tests for weak Instruments
}
@@ -9,11 +10,12 @@
\description{
This is a collection of tests for weak instruments. It includes the
Cragg and Donald test for the reduced rank hypothesis, the conditional
-F-test of Sanderson and Windmeijer (2016) and the robust effective F
-test of Montiel Olea and Pflueger (2013). The critical values of Stock
-and Yogo (2005) for the null hypothesis of weak instruments are also
-provided.
-}
+F-test of Sanderson and Windmeijer (2016), the robust effective F
+test of Montiel Olea and Pflueger (2013) and the Lewis and Mertens
+(2022) robust test with multiple endogenous variables. The critical
+values of Stock and Yogo (2005) for the null hypothesis of weak
+instruments are also provided.}
+
\usage{
SWtest(object, j=1, print=TRUE, ...)
@@ -24,6 +26,10 @@
estMethod = c("TSLS", "LIML"), simplified = TRUE,
digits = max(3L, getOption("digits") - 3L), ...)
+LewMertest(object, tau=0.10, size=0.05, print=TRUE,
+ estMethod = c("TSLS", "LIML"), simplified = TRUE,
+ digits = max(3L, getOption("digits") - 3L), ...)
+
SYTables(object, print=TRUE, SWcrit=FALSE)
}
Modified: pkg/momentfit/vignettes/weak.Rmd
===================================================================
--- pkg/momentfit/vignettes/weak.Rmd 2024-02-09 22:15:30 UTC (rev 225)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-02-16 22:44:26 UTC (rev 226)
@@ -748,7 +748,7 @@
momentStrength(update(mod4, vcovOptions=list(type="HC0")))$strength
```
-## Testing for weak instruments: @lewis-mertens22
+## Testing for weak instruments: @lewis-mertens22 (very early stage)
The authors generalize the MOP test for models with multiple
endogenous variables. We first need to compute the matrix $\Phi =
@@ -774,30 +774,8 @@
t(A)%*%B%*%A
}
-phiMat <- function(w2, k2, l2)
-{
- psi <- matrix(NA, k2,k2)
- for (i in 1:k2)
- {
- c0 <- 1+(i-1)*l2
- c1 <- ncol(w2)
- di <- diag(w2[,c0:c1])
- tmp <- colSums(matrix(di, nrow=l2))
- if (i == 1)
- {
- diag(psi) <- tmp
- } else if (i==k2) {
- psi[1,k2] <- psi[k2,1] <- tmp
- } else {
- j <- seq_len(length(tmp))
- j <- cbind(j,j)
- psi[-(1:(i-1)),][j] <- psi[,-(1:(i-1))][j] <- tmp
- }
- }
- psi
-}
phiMatT(w2,k2,l2)
-phiMat(w2,k2,l2)
+momentfit:::phiMat(w2,k2,l2)
```
The statistic proposed by the authors is
@@ -810,56 +788,14 @@
because $k_2=1$ implies that $\Phi=\tr(W_2)$.
```{r}
-LewMertest <- function(object, tau=0.10, size=0.05, print=TRUE,
- estMethod = c("TSLS", "LIML"), simplified = TRUE,
- digits = max(3L, getOption("digits") - 3L), ...)
-{
- estMethod <- match.arg(estMethod)
- if (!inherits(object, "linearModel"))
- stop("object must be of class linearModel")
- spec <- modelDims(object)
- if (sum(spec$isEndo)<1)
- {
- warning("The model does not contain endogenous variables")
- return(NA)
- }
- Z2 <- model.matrix(object, "excludedExo")
- X1 <- model.matrix(object, "includedExo")
- X2 <- model.matrix(object, "includedEndo")
- y <- modelResponse(object)
- if (!is.null(X1))
- {
- fitX1 <- lm.fit(X1, Z2)
- 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="")
- colnames(X2) <- paste("endo",1:ncol(X2), sep="")
- b1 <- c(crossprod(Z2,y)/spec$n)
- Pi2 <- crossprod(Z2,X2)/spec$n
- g <- c(reformulate(colnames(Z2), "y", FALSE),
- lapply(colnames(X2), function(ei)
- reformulate(colnames(Z2), ei, 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%*%Pi2)
- omega <- crossprod(v)/nrow(v)
- w <- vcov(m, c(list(b1),lapply(1:ncol(Pi2), function(i) Pi2[,i])))
- w2 <- w[(ncol(Z2)+1):ncol(w),(ncol(Z2)+1):ncol(w)]
- phi <- phiMat(w2, ncol(X2), ncol(Z2))
- gmin <- eigen(solve(phi, nrow(Z2)*crossprod(Pi2)))$val[1]
- gmin
-}
-LewMertest(mod2)
-LewMertest(mod3)
+LewMertest(mod2)$gmin
+LewMertest(mod4)$gmin
```
+
+
## 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