[Gmm-commits] r228 - in pkg/momentfit: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 23 23:44:24 CET 2024
Author: chaussep
Date: 2024-02-23 23:44:23 +0100 (Fri, 23 Feb 2024)
New Revision: 228
Modified:
pkg/momentfit/NAMESPACE
pkg/momentfit/R/sysMomentModel-methods.R
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 and Mertens test for weak instruments -- from their matlab code
Modified: pkg/momentfit/NAMESPACE
===================================================================
--- pkg/momentfit/NAMESPACE 2024-02-21 19:47:21 UTC (rev 227)
+++ pkg/momentfit/NAMESPACE 2024-02-23 22:44:23 UTC (rev 228)
@@ -16,7 +16,8 @@
"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", "optimize")
+ "constrOptim", "kernel", "confint", "qnorm", "uniroot", "getCall", "qchisq", "optimize",
+ "integrate", "rnorm")
importFrom("sandwich", "vcovHAC", "estfun","kernHAC","vcovCL", "meatCL",
"bread","bwAndrews","bwNeweyWest","weightsAndrews",
"weightsLumley", "vcovHC", "meatHAC")
Modified: pkg/momentfit/R/sysMomentModel-methods.R
===================================================================
--- pkg/momentfit/R/sysMomentModel-methods.R 2024-02-21 19:47:21 UTC (rev 227)
+++ pkg/momentfit/R/sysMomentModel-methods.R 2024-02-23 22:44:23 UTC (rev 228)
@@ -758,9 +758,9 @@
}
weights <- weightsAndrews(x = gmat, bw = bw, kernel = options$kernel,
prewhite = options$prewhite, tol = options$tol)
- w <- vcovHAC(x = gmat, order.by = NULL, weights = weights,
- prewhite = options$prewhite, sandwich = FALSE,
- ar.method = options$ar.method)
+ w <- meatHAC(x = gmat, order.by = NULL, weights = weights,
+ prewhite = options$prewhite, diagnostics = FALSE,
+ adjust = options$adjust, ar.method = options$ar.method)
attr(w, "Spec") <- list(weights = weights, bw = bw, kernel = options$kernel)
w
})
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-02-21 19:47:21 UTC (rev 227)
+++ pkg/momentfit/R/weak.R 2024-02-23 22:44:23 UTC (rev 228)
@@ -365,6 +365,50 @@
# computing x for the generalized test
+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)
+ }
+ 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="")
+ b <- c(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 <- vcov(m, list(b1,b2))
+ 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)
+}
+
+
getMOPx <- function(w, tau, type = c("TSLS", "LIML"), e=0.0001, nP = 10000,
maxi = 1000)
{
@@ -521,50 +565,6 @@
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)
- }
- 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="")
- b <- c(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 <- vcov(m, list(b1,b2))
- 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)
-}
-
-
## Lewis and Mertens (2022)
phiMat <- function(w2, k2, l2)
@@ -593,11 +593,10 @@
}
LewMertest <- function(object, tau=0.10, size=0.05, print=TRUE,
- estMethod = c("TSLS", "LIML"), simplified = TRUE,
+ simplified = TRUE,
digits = max(3L, getOption("digits") - 3L),
- dfCorr = TRUE, ...)
+ npoints=10, ...)
{
- estMethod <- match.arg(estMethod)
if (!inherits(object, "linearModel"))
stop("object must be of class linearModel")
spec <- modelDims(object)
@@ -627,13 +626,14 @@
reformulate(colnames(Z2), ei, FALSE)))
h <- reformulate(colnames(Z2), NULL, FALSE)
dat <- as.data.frame(cbind(y=y,X2,Z2))
+ if (object at vcov == "HAC")
+ object at vcovOptions$adjust <- FALSE
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])))
- if (dfCorr)
- w <- w*nrow(Z2)/(nrow(Z2)-ncol(X1)-ncol(Z2))
+ w <- w*nrow(Z2)/(nrow(Z2)-ncol(X1)-ncol(Z2))
w2 <- w[(ncol(Z2)+1):ncol(w),(ncol(Z2)+1):ncol(w)]
phi <- phiMat(w2, ncol(X2), ncol(Z2))
if (dim(phi)[1] == 1)
@@ -645,50 +645,368 @@
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))
+ crit <- LewMerCrit(W=w, K=ncol(Z2), alpha=size, tau=tau,
+ points=npoints, generalized=!simplified)
+ if (!print)
+ return(list(test=gmin, crit=crit))
+ cat("Lewis and Mertens (2022) Test for Weak Instruments\n")
+ cat("***************************************************\n", sep="")
+ cat("Number of included Endogenous: ", ncol(X2), "\n", sep="")
+ cat("Number of excluded Exogenous: ", ncol(Z2), "\n", sep="")
+ cat("Statistics: ", formatC(gmin, digits=digits, ...), "\n", sep="")
+ cat("Simplified Critical value (", size*100, "%): ",
+ formatC(crit["Simplified"], digits=digits, ...), "\n", sep="")
+ if (!simplified)
+ cat("Generalized Critical value (", size*100, "%): ",
+ formatC(crit["Generalized"], digits=digits, ...), "\n", sep="")
+ invisible()
}
-getCumul <- function(w2, lam, phi, l2)
+
+
+nearestSPD <- function (A)
{
- if (length(phi)==1)
+ stopifnot(is.numeric(A), is.matrix(A))
+ eps <- .Machine$double.eps
+ m <- nrow(A)
+ n <- ncol(A)
+ if (m != n) {
+ stop("Argument 'A' must be a square matrix.")
+ }
+ else if (n == 1 && A <= 0)
+ return(as.matrix(eps))
+ B <- (A + t(A))/2
+ svdB <- svd(B)
+ H <- svdB$v %*% diag(svdB$d) %*% t(svdB$v)
+ Ahat <- (B + H)/2
+ Ahat <- (Ahat + t(Ahat))/2
+ k <- 0
+ not_pd <- TRUE
+ while (not_pd) {
+ k <- k + 1
+ try_R <- try(chol(Ahat), silent = TRUE)
+ if (inherits(try_R, "try-error")) {
+ mineig <- min(eigen(Ahat, symmetric = TRUE, only.values = TRUE)$values)
+ Ahat = Ahat + (-mineig * k^2 + eps(mineig)) * diag(1,
+ n)
+ }
+ else not_pd <- FALSE
+ }
+ Ahat
+}
+
+LewMerCrit <- function(model, W, K, alpha=0.05, tau=0.1, points=10,
+ generalized=FALSE)
+{
+ if (!missing(model))
{
- 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)
+ spec <- modelDims(model)
+ if (sum(spec$isEndo)<1)
+ {
+ warning("The model does not contain endogenous variables")
+ return(NA)
+ }
+ Z2 <- model.matrix(model, "excludedExo")
+ X1 <- model.matrix(model, "includedExo")
+ X2 <- model.matrix(model, "includedEndo")
+ y <- modelResponse(model)
+ 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))
+ if (model at vcov == "HAC")
+ model at vcovOptions$adjust <- FALSE
+ m <- sysMomentModel(g=g, list(h), data = dat, vcov=model at vcov,
+ vcovOptions = model 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])))
+ W <- W*nrow(Z2)/(nrow(Z2)-ncol(X1)-ncol(Z2))
+ K <- ncol(Z2)
}
- k1 <- l2*(1+lam)
- if (length(sig) == 1)
+ spKgen <- function(m,n)
{
- k2 <- 2*(sig^2+2*lam*l2*sig)
- k3 <- 8*(sig^3+3*lam*l2*sig^2)
+ I <- c(matrix(1:(m*n), n, m, byrow=TRUE))
+ K <- diag(m*n)[I,]
+ K
+ }
+ fMat <- function(A, FUN, eA=NULL)
+ {
+ if (is.null(eA))
+ eA <- eigen(A)
+ eA$vec%*%diag(FUN(eA$val))%*%t(eA$vec)
+ }
+
+ N <- nrow(W)/K-1
+ maxiter <- 1000
+ tol <- 1e-7
+ Bmax <- rep(NA,2)
+
+ RNK <- kronecker(diag(N), c(diag(K)))
+ RNN <- kronecker(diag(N), c(diag(N)))
+ RNpK <- kronecker(diag(N+1), c(diag(K)))
+ M2 <- RNK%*%t(RNK)/(1+N)-diag(N*K^2)
+ W1 <- W[1:K,1:K]
+ W2 <- W[-(1:K),-(1:K)]
+ W12 <- W[1:K,-(1:K)]
+ tmp <- crossprod(RNK, kronecker(W2, diag(K))%*%RNK)/K
+ tmp <- fMat(tmp, function(x) 1/sqrt(x))
+ tmp <- kronecker(tmp, diag(K))
+ tmp2 <- fMat(W2, sqrt)
+ S <- tmp%*%tmp2
+ Sigma <- S%*%t(S)
+ tmp2 <- crossprod(RNpK, kronecker(W,diag(K))%*%RNpK)
+ tmp2 <- fMat(tmp2, function(x) 1/sqrt(x))
+ Psi <- kronecker(tmp%*%t(rbind(W12,W2)),diag(K))%*%RNpK%*%tmp2
+
+ ## Simplified
+ if (N==1)
+ {
+ if (K>N+1)
+ {
+ Bmax[2] <- min(sqrt(2*(N+1)/K)*norm(M2%*%Psi, "2"),1);
+ } else {
+ Bmax[2] <- min(norm(Psi, "2"),1);
+ }
} 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)
+ if (K>N+1)
+ {
+ Bmax[2] <- min(sqrt(2*(N+1)/K)*norm(M2%*%Psi, "2"), norm(Psi, "2"));
+ } else {
+ Bmax[2] <- norm(Psi, "2");
+ }
+ }
+ ## Generalized
+ if (generalized)
+ {
+ if (K>N+1)
+ {
+ bm.iter <- numeric(points)
+ X1 <- kronecker(kronecker(diag(N),spKgen(K^2,N)),diag(N^2))%*%
+ kronecker(c(diag(N)),diag(K^2*N^2))%*%
+ kronecker(kronecker(diag(K),spKgen(K,N)),diag(N))%*%
+ (diag(N^2*K^2)+spKgen(N*K,N*K))
+ M1 <- crossprod(RNN, diag(N^3)+kronecker(spKgen(N,N),diag(N)))
+ M2PsiM2 <- M2%*%(Psi%*%t(Psi))%*%t(M2)
+ for (i in 1:points)
+ bm.iter[i] <- sqrt(-.optStiefel(M1,M2PsiM2,X1,N,K)$fval)
+ Bmax[1] <- max(bm.iter)
+ } else {
+ Bmax[1] <- Bmax[2]
+ }
}
- w <- k2/k3
- v <- 8*k2*w^2
- list(k1=k1, k2=k2, k3=k3, w=w, v=v)
+ cv <- rep(NA,2)
+ for (j in 1:2)
+ {
+ if (!is.na(Bmax[j]))
+ {
+ lmin <- Bmax[j]/tau
+ k <- numeric(3)
+ k[1] <- norm(crossprod(RNK,kronecker(Sigma,diag(K)))%*%RNK,"2")+K*lmin
+ k[2] <- 2*(norm(crossprod(RNK,kronecker(fMat(Sigma, function(x) x^2),diag(K)))%*%RNK,"2")+
+ 2*K*lmin*norm(Sigma,"2"))
+ k[3] <- 8*(norm(crossprod(RNK,kronecker(fMat(Sigma, function(x) x^3),diag(K)))%*%RNK,"2")+
+ 3*K*lmin*norm(Sigma,"2")^2)
+ cv[j] <- .getcritical(k, alpha, K)
+ }
+ }
+ names(cv) <- c("Generalized", "Simplified")
+ cv
}
-crit <- function(alpha, k1, k2, k3)
+.optStiefel <- function(M1,M2PsiM2,X1,N,K)
{
- w <- k2/k3
- v <- 8*k2*w^2
- cr <- qchisq(1-alpha, v)
- (cr-v)/(4*w)+k1
+ Fct <- function(x,M1,M2PsiM2,X1,N,K)
+ {
+ L0 <- t(x)
+ vecL0 <- as.matrix(c(L0))
+ QLL <- kronecker(kronecker(diag(N),L0),L0)
+ Mobj <- M1%*%QLL%*%M2PsiM2%*%t(QLL)%*%t(M1)/K
+ Mobj <- 0.5*(Mobj+t(Mobj))
+ Mobj <- nearestSPD(Mobj)
+ res <- eigen(Mobj)
+ ev <- res$vec[,which.max(res$val)]
+ fval <- -t(ev)%*%Mobj%*%ev
+ gradient = 2*kronecker(t(ev)%*%M1%*%QLL%*%M2PsiM2,t(ev)%*%M1)%*%
+ X1%*%kronecker(diag(N*K),vecL0)
+ gradient = -matrix(c(gradient),N,K)
+ gradient = t(gradient)
+ list(fct=c(fval), grad=gradient)
+ }
+ rustiefel <- function (m, R)
+ {
+ X <- matrix(rnorm(m * R), m, R)
+ tmp <- eigen(t(X) %*% X)
+ X %*% (tmp$vec %*% sqrt(diag(1/tmp$val, nrow = R)) %*% t(tmp$vec))
+ }
+ myQR <- function(XX,k)
+ {
+ res <- qr(XX)
+ Q <- qr.Q(res)
+ RR <- qr.R(res)
+ diagRR <- sign(diag(RR))
+ if (any(diagRR < 0))
+ Q <- Q%*%diag(diagRR)
+ list(Q=Q, R=RR)
+ }
+
+ iprod <- function(x,y)
+ {
+ sum(sum(Conj(x)*y))
+ }
+
+ X <- rustiefel(K, N)
+ opts <- list()
+ out <- list()
+
+ opts$tau <- 1e-3
+ opts$mxitr <- 1000
+ opts$record <- 0
+ xtol <- 1e-6
+ gtol <- 1e-6
+ ftol <- 1e-12
+ rhols <- 1e-4
+ STPEPS <- 1e-10
+ eta <- 0.1
+ gamma <- 0.85
+ retr <- 0
+ record <- opts$record
+ nt <- 5
+ crit <- numeric()
+ tiny <- 1e-13
+
+ FX <- Fct(X,M1,M2PsiM2,X1,N,K)
+ F <- FX$fct
+ G <- FX$grad
+ out$nfe = 1
+ GX = crossprod(G,X)
+
+ dtX = G - X%*%GX
+ nrmG <- norm(dtX, 'F')
+
+ Q <- 1; Cval <- F; tau <-opts$tau
+ for (itr in 1:opts$mxitr)
+ {
+ XP <- X
+ FP <- F
+ GP <- G
+ dtXP <- dtX
+ nls <- 1; deriv <- rhols*nrmG^2
+ while (TRUE)
+ {
+ X = myQR(XP - tau*dtX, N)$Q
+ if (norm(crossprod(X) - diag(N),'F') > tiny)
+ X <- myQR(X,N)
+ FX2 <- Fct(X,M1,M2PsiM2,X1,N,K)
+ F <- FX2$fct
+ G <- FX2$grad
+ out$nfe <- out$nfe + 1
+ if (F <= (Cval - tau*deriv) || nls >= 5)
+ break;
+ tau <- eta*tau
+ nls <- nls+1
+ }
+ GX <- crossprod(G,X)
+ dtX <- G - X%*%GX
+ nrmG <- norm(dtX, 'F')
+ S <- X - XP
+ XDiff <- norm(S,'F')/sqrt(K)
+ tau = opts$tau
+ FDiff = abs(FP-F)/(abs(FP)+1)
+ Y = dtX - dtXP
+ SY = abs(iprod(S,Y))
+ if ((itr%%2)==0)
+ {
+ tau <- norm(S,'F')^2/SY
+ } else {
+ tau <- SY/norm(Y,'F')^2
+ }
+ tau <- max(min(tau, 1e20), 1e-20)
+ crit <- rbind(crit, c(nrmG, XDiff, FDiff))
+ mcrit = colMeans(crit[(itr-min(nt,itr)+1):itr,,drop=FALSE])
+ if ((XDiff < xtol && FDiff < ftol ) ||
+ (nrmG < gtol) || all(mcrit[2:3] < 10*c(xtol, ftol)))
+ {
+ out$msg = "converge"
+ break
+ }
+ Qp <- Q
+ Q <- gamma*Qp + 1
+ Cval = (gamma*Qp*Cval + F)/Q
+
+ }
+ if (itr >= opts$mxitr)
+ out$msg = "exceed max iteration"
+
+ out$feasi = norm(crossprod(X)-diag(N),'F');
+ if (out$feasi > 1e-13)
+ {
+ X <- myQR(X,N);
+ FX <- Fct(X,M1,M2PsiM2,X1,N,K)
+ F <- FX$fct
+ G <- FX$grad
+ out$nfe <- out$nfe + 1;
+ out$feasi <- norm(crossprod(X)-diag(N),'F');
+ }
+
+ out$nrmG <- nrmG
+ out$fval <- F
+ out$itr <- itr
+ out
}
+.getcritical <- function(k, alpha, K)
+{
+ ome <- k[2]/k[3]
+ nu <- 8*k[2]*ome^2
+ cc <- qchisq(1-alpha,nu)
+ fun_phiz <- function(z)
+ ome*(1+(z-k[1])/(2*k[2]*ome))^(nu/2-1)*
+ exp(-nu/2*(1+(z-k[1])/(2*k[2]*ome)))*
+ nu^(nu/2-1)/(2^(nu/2-2))/gamma(nu/2)
+ G1fun <- function(q)
+ -1/2*(q-2*nu*(nu-2)/q+nu)+3*nu/2*((log(q/2))-psigamma(nu/2,0))
+ G2fun <- function(q)
+ 1/2*(q-nu*(nu-2)/q)-nu*((log(q/2))-psigamma(nu/2,0))
+ D1fun <- function(q)
+ (1+(q-k[1])*2*ome)/(2*k[2]*ome)*(1+(q-k[1])/(2*k[2]*ome))^(-1)*fun_phiz(q)
+ D2fun <- function(q)
+ fun_phiz(q)/k[2]*G1fun(nu+(q-k[1])*4*ome)
+ D3fun <- function(q)
+ G2fun(nu+(q-k[1])*4*ome)/k[3]*fun_phiz(q)
+
+ ID1fun <- function(xx) integrate(D1fun,xx,Inf)$value
+ ID2fun <- function(xx) integrate(D2fun,xx,Inf)$value
+ ID3fun <- function(xx) integrate(D3fun,xx,Inf)$value
+ kt_cond1 <- ID1fun(((cc-nu)/4/ome+k[1]))
+ kt_cond2 <- ID2fun(((cc-nu)/4/ome+k[1]))
+ kt_cond3 <- ID3fun(((cc-nu)/4/ome+k[1]))
+ kt_cond <- (kt_cond1>=0)&&(kt_cond2>=0)&&(kt_cond3>=0)
+ if (!kt_cond)
+ {
+ fun = function(x)
+ -((qchisq(1-alpha,8*x[2]*(x[2]/x[3])^2)-8*x[2]*(x[2]/x[3])^2)/4/(x[2]/x[3])+x[1])
+ k <- nlminb(k, fun, upper=k)$par
+ ome <- k[2]/k[3]
+ nu <- 8*k[2]*ome^2
+ cc <- qchisq(1-alpha,nu)
+ }
+ ((cc-nu)/4/ome+k[1])/K
+}
+
+
Modified: pkg/momentfit/man/weakTest.Rd
===================================================================
--- pkg/momentfit/man/weakTest.Rd 2024-02-21 19:47:21 UTC (rev 227)
+++ pkg/momentfit/man/weakTest.Rd 2024-02-23 22:44:23 UTC (rev 228)
@@ -27,9 +27,8 @@
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),
- dfCorr = TRUE, ...)
+ simplified = TRUE, digits = max(3L, getOption("digits") -
+ 3L), npoints=10, ...)
SYTables(object, print=TRUE, SWcrit=FALSE)
@@ -67,14 +66,17 @@
\item{digits}{The number of digits to print when \code{print} is set
to \code{TRUE}.}
-\item{dfCorr}{Should we correct the covariance matrix estimator for the
- degrees of freedom?}
-
\item{SWcrit}{If true, the critical values and the test are ajusted to
Sanderson and Windmeijer (2016). See details.}
-
+
+\item{npoints}{The number of minimizations to run in order to find the
+ global minimum, which is used for the Lewis and Mertens critical
+ value. It affects the results only when \code{simplified} is set to
+ \code{TRUE}.}
+
\item{...}{Arguments passed to \code{\link{formatC}} to print the
- statistic. For \code{MOPtest}, they are passed to \code{momentfit:::getMOPx}.
+ statistic. For \code{MOPtest}, they are passed to
+ \code{momentfit:::getMOPx}.
}
}
Modified: pkg/momentfit/vignettes/weak.Rmd
===================================================================
--- pkg/momentfit/vignettes/weak.Rmd 2024-02-21 19:47:21 UTC (rev 227)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-02-23 22:44:23 UTC (rev 228)
@@ -748,54 +748,27 @@
momentStrength(update(mod4, vcovOptions=list(type="HC0")))$strength
```
-## Testing for weak instruments: @lewis-mertens22 (very early stage)
+## Testing for weak instruments: @lewis-mertens22
The authors generalize the MOP test for models with multiple
-endogenous variables. We first need to compute the matrix $\Phi =
-(I_{k_2}\otimes \vect(I_{l_2}))'[W_2\otimes I_{l_2}](I_{k_2}\otimes
-\vect(I_{l_2}))$. This is clearly not something we want to write
-explicitely. It is unlikely to encounter very large $l_2$ and $k_2$,
-but we want to avoid creating large matrices just in case. If we write
-$W_2$ as a $k_2\times k_2$ block matrix, with $A_{ij}$ being the
-$l_2\times l_2$ matrix in the ith row and jth column, then
-$\Phi_{ij}=\tr(A_{ij})$. The following shows that the function
-`psiMat`, which uses this trace argument, produces the right matrix:
+endogenous variables. The statistic proposed by the authors is
-
-```{r}
-w2 <- matrix(1:144, 12, 12)
-w2[lower.tri(w2)] <- t(w2)[lower.tri(w2)]
-k2 <- 3
-l2 <- 4
-phiMatT <- function(w2, k2, l2)
-{
- A <- kronecker(diag(k2), c(diag(l2)))
- B <- kronecker(w2, diag(l2))
- t(A)%*%B%*%A
-}
-
-phiMatT(w2,k2,l2)
-momentfit:::phiMat(w2,k2,l2)
-```
-
-The statistic proposed by the authors is
-
-\[
-g_{min} = \min\eval[n\hat\Phi^{-1/2}\hat{\Pi}_2'\hat{\Pi}_2\hat\Phi^{-1/2}]
+\[ g_{min} =
+\min\eval[n\hat\Phi^{-1/2}\hat{\Pi}_2'\hat{\Pi}_2\hat\Phi^{-1/2}] \,,
\]
-It is clear that this is the MOP effective F test when $k_2=1$,
-because $k_2=1$ implies that $\Phi=\tr(W_2)$.
+where $\Phi = (I_{k_2}\otimes \vect(I_{l_2}))'[W_2\otimes
+I_{l_2}](I_{k_2}\otimes \vect(I_{l_2}))$. If we write $W_2$ as a
+$k_2\times k_2$ block matrix, with $A_{ij}$ being the $l_2\times l_2$
+matrix in the ith row and jth column, then $\Phi_{ij}=\tr(A_{ij})$. It
+is clear that this is the MOP effective F test when $k_2=1$, because
+$k_2=1$ implies that $\Phi=\tr(W_2)$.
```{r}
-LewMertest(mod2)$gmin
-LewMertest(mod4)$gmin
+LewMertest(mod3, simplified=FALSE)
```
-
-
-
## 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