[Gmm-commits] r222 - in pkg/momentfit: R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 7 18:01:29 CET 2024
Author: chaussep
Date: 2024-02-07 18:01:28 +0100 (Wed, 07 Feb 2024)
New Revision: 222
Modified:
pkg/momentfit/R/momentModel-methods.R
pkg/momentfit/R/sysMomentModel-methods.R
pkg/momentfit/R/weak.R
pkg/momentfit/man/vcovHAC-methods.Rd
pkg/momentfit/vignettes/weak.Rmd
pkg/momentfit/vignettes/weak.pdf
Log:
working on weak instrument tests and added HAV and CL to system of equations
Modified: pkg/momentfit/R/momentModel-methods.R
===================================================================
--- pkg/momentfit/R/momentModel-methods.R 2024-02-06 15:02:14 UTC (rev 221)
+++ pkg/momentfit/R/momentModel-methods.R 2024-02-07 17:01:28 UTC (rev 222)
@@ -129,27 +129,25 @@
if (type == "regressors")
{
ti <- attr(object at modelF, "terms")
- mat <- as.matrix(model.matrix(ti, object at modelF)[,])
+ mat <- as.matrix(model.matrix(ti, object at modelF)[,,drop=FALSE])
} else if (type == "instruments"){
ti <- attr(object at instF, "terms")
- mat <- as.matrix(model.matrix(ti, object at instF)[,])
+ mat <- as.matrix(model.matrix(ti, object at instF)[,,drop=FALSE])
} else {
X <- model.matrix(object)
Z <- model.matrix(object, "instruments")
- z1 <- colnames(Z) %in% colnames(X)
+ x1n <- colnames(Z) %in% colnames(X)
endo <- modelDims(object)$isEndo
if (type == "excludedExo")
{
- if (all(!z1))
+ if (all(x1n))
{
- warning("No excluded Exogenous Variables")
return(NULL)
}
- mat <- Z[,!z1,drop=FALSE]
+ mat <- Z[,!x1n,drop=FALSE]
} else if (type == "includedExo") {
if (all(endo))
{
- warning("No included Exogenous Variables")
return(NULL)
}
mat <- X[,!endo,drop=FALSE]
@@ -156,7 +154,6 @@
} else {
if (all(!endo))
{
- warning("No included endogenous variables")
return(NULL)
}
mat <- X[,endo,drop=FALSE]
Modified: pkg/momentfit/R/sysMomentModel-methods.R
===================================================================
--- pkg/momentfit/R/sysMomentModel-methods.R 2024-02-06 15:02:14 UTC (rev 221)
+++ pkg/momentfit/R/sysMomentModel-methods.R 2024-02-07 17:01:28 UTC (rev 222)
@@ -499,8 +499,9 @@
### The following multiplies each block matrix of ZZ, when the dimensions
### are defined by dimr and dimc (dim row and dim column) by each element
### of Sigma. Sigma must be length(dimr) by length(dimc)
+### Just set the lowerTri to TRUE by default. The cost is negligible.
-.SigmaZZ <- function(ZZ, Sigma, dimr, dimc=NULL, lowerTri=FALSE, isSym=TRUE)
+.SigmaZZ <- function(ZZ, Sigma, dimr, dimc=NULL, lowerTri=TRUE, isSym=TRUE)
{
r1 <- 1
c1 <- 1
@@ -723,6 +724,47 @@
met(object, wObj, theta0, ...)
})
+## vcovHAC
+
+setMethod("vcovHAC", "sysModel",
+ function (x, theta) {
+ if (x at vcov != "HAC")
+ {
+ warning("Model set as ", x at vcov, ". The default HAC options are used")
+ x at vcov <- "HAC"
+ x at sSpec <- new("sSpec")
+ x at smooth <- FALSE
+ x at vcovOptions <- .getVcovOptions("HAC")
+ }
+ gmat <- evalMoment(x, theta)
+ gmat <- do.call(cbind, gmat)
+ if (x at centeredVcov)
+ gmat <- scale(gmat, scale = FALSE)
+ class(gmat) <- "momentFct"
+ options <- x at vcovOptions
+ if (is.character(options$bw))
+ {
+ if (options$bw == "Wilhelm")
+ warning("bw = Wilhelm is not available for system GMM")
+ obj <- gmat
+ bwFct <- get(paste("bw",options$bw,sep=""))
+ bwArgs <- options
+ bwArgs$bw <- NULL
+ bwArgs$tol <- NULL
+ bwArgs$x <- obj
+ bw <- do.call(bwFct, bwArgs)
+ } else {
+ bw <- options$bw
+ }
+ 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)
+ attr(w, "Spec") <- list(weights = weights, bw = bw, kernel = options$kernel)
+ w
+ })
+
## vcov
setMethod("vcov", signature("sysModel"),
@@ -730,26 +772,33 @@
spec <- modelDims(object)
q <- spec$q
if (object at vcov == "MDS")
+ {
+ gt <- evalMoment(object, theta)
+ gt <- do.call(cbind, gt)
+ if (object at centeredVcov)
+ gt <- scale(gt, scale=FALSE)
+ w <- crossprod(gt)/nrow(gt)
+ } else if (object at vcov == "iid") {
+ e <- residuals(object, theta)
+ Sigma <- crossprod(e)/nrow(e)
+ Z <- model.matrix(object, "instrument")
+ if (object at sameMom)
{
- gt <- evalMoment(object, theta)
- gt <- do.call(cbind, gt)
- if (object at centeredVcov)
- gt <- scale(gt, scale=FALSE)
- w <- crossprod(gt)/nrow(gt)
- } else if (object at vcov == "iid") {
- e <- residuals(object, theta)
- Sigma <- crossprod(e)/nrow(e)
- Z <- model.matrix(object, "instrument")
- if (object at sameMom)
- {
- w <- kronecker(Sigma, crossprod(Z[[1]])/nrow(e))
- } else {
- Z <- crossprod(do.call(cbind,Z))/nrow(e)
- w <- .SigmaZZ(Z, Sigma, q)
- }
+ w <- kronecker(Sigma, crossprod(Z[[1]])/nrow(e))
} else {
- stop("not yet implemented for HAC")
+ Z <- crossprod(do.call(cbind,Z))/nrow(e)
+ w <- .SigmaZZ(Z, Sigma, q)
}
+ } else if (object at vcov == "CL") {
+ gt <- evalMoment(object, theta)
+ gt <- do.call(cbind, gt)
+ class(gt) <- "momentFct"
+ opt <- object at vcovOptions
+ opt$x <- gt
+ w <- do.call(meatCL, opt)
+ } else {
+ w <- vcovHAC(object, theta)
+ }
wn <- paste(rep(spec$eqnNames, q), ".", do.call("c", spec$momNames),
sep = "")
dimnames(w) <- list(wn,wn)
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-02-06 15:02:14 UTC (rev 221)
+++ pkg/momentfit/R/weak.R 2024-02-07 17:01:28 UTC (rev 222)
@@ -27,7 +27,6 @@
## K-Class functions
-
getK <- function(object, alpha=1, returnRes=FALSE)
{
if (!inherits(object, "linearModel"))
@@ -38,12 +37,16 @@
{
k1 <- 1
} else {
- X <- model.matrix(object)
- endo <- cbind(model.response(object at modelF),
- X[,object at isEndo])
- X <- X[,!object at isEndo, drop=FALSE]
+ X2 <- model.matrix(object, "includedEndo")
+ endo <- cbind(modelResponse(object), X2)
+ X <- model.matrix(object, "includedExo")
Z <- model.matrix(object, "instruments")
- e1 <- lm.fit(X, endo)$residuals
+ if (!is.null(X))
+ {
+ e1 <- lm.fit(X, endo)$residuals
+ } else {
+ e1 <- endo
+ }
e2 <- lm.fit(Z, endo)$residuals
if (returnRes) Res <- e2[,-1,drop=FALSE]
e1 <- crossprod(e1)
@@ -199,12 +202,11 @@
stop("object must be of class linearModel")
spec <- modelDims(object)
Z <- model.matrix(object, "instrument")
- X <- model.matrix(object)
- z2n <- !(colnames(Z) %in% colnames(X))
- X1 <- X[,!spec$isEndo, drop=FALSE]
- X2 <- X[,spec$isEndo, drop=FALSE]
+ Z2 <- model.matrix(object, "excludedExo")
+ X1 <- model.matrix(object, "includedExo")
+ X2 <- model.matrix(object, "includedEndo")
+ z2n <- colnames(Z) %in% colnames(Z2)
secS <- lm.fit(Z, X2)
- Z2 <- Z[,z2n,drop=FALSE]
df <- ncol(Z2)
if (ncol(Z2)==1 & SWcrit)
{
@@ -218,7 +220,7 @@
else
Z[,z2n,drop=FALSE]%*%secS$coef[z2n]
e <- secS$residuals
- e2 <- lm.fit(X1, X2)$residuals
+ e2 <- if (!is.null(X1)) lm.fit(X1, X2)$residuals else X2
e <- crossprod(e)/(spec$n-spec$q)
e2 <- crossprod(e2, tmp)
test <- min(eigen(solve(e,e2))$val)/df
@@ -327,26 +329,24 @@
"Returning the F-test", sep=""))
return(CDtest(object, print))
}
- Z <- model.matrix(object, "instrument")
- X <- model.matrix(object)
- X2 <- X[, spec$isEndo, drop=FALSE]
- X1 <- X[, !spec$isEndo, drop=FALSE]
- if (ncol(X1))
+ Z2 <- model.matrix(object, "excludedExo")
+ X1 <- model.matrix(object, "includedExo")
+ X2 <- model.matrix(object, "includedEndo")
+ if (!is.null(X1))
{
- z2n <- !(colnames(Z) %in% colnames(X1))
- fitX1 <- lm.fit(X1, Z[,z2n])
- Z <- as.matrix(fitX1$residuals)
+ fitX1 <- lm.fit(X1, Z2)
+ Z2 <- fitX1$residuals
X2 <- qr.resid(fitX1$qr, X2)
}
Xj <- X2[,j]
Xjm <- X2[,-j,drop=FALSE]
- fsReg <- lm.fit(Z, Xjm)
+ fsReg <- lm.fit(Z2, Xjm)
Xjmhat <- as.matrix(fsReg$fitted)
fit <- lm.fit(Xjmhat, Xj)
e <- Xj-c(Xjm%*%fit$coefficients)
ehat <- qr.fitted(fsReg$qr, e)
sig <- sum((e-ehat)^2)/(spec$n-spec$q)
- test <- sum(ehat^2)/sig/(ncol(Z)-ncol(Xjm))
+ test <- sum(ehat^2)/sig/(ncol(Z2)-ncol(Xjm))
if (!print)
return(test)
cat("Sanderson and Windmeijer Test for Weak Instruments\n")
@@ -354,7 +354,7 @@
add1 <- "(-1 for the critical values)"
add2 <- paste("(-", ncol(X2)-1, " for the critical values)", sep="")
cat("Number of included Endogenous: ", ncol(X2), add1, "\n", sep="")
- cat("Number of excluded Exogenous: ", sum(z2n), add2, "\n", sep="")
+ cat("Number of excluded Exogenous: ", ncol(Z2), add2, "\n", sep="")
cat("The test is not robust to heteroskedasticity\n")
cat("Statistics: ", formatC(test, ...), "\n\n", sep="")
SYTables(object, TRUE, TRUE)
@@ -367,7 +367,7 @@
{
if (!inherits(object, "linearModel"))
stop("object must be of class linearModel")
- spec <- modelDims(object)
+ spec <- modelDims(object)
if (sum(spec$isEndo)<1)
{
warning("The model does not contain endogenous variables")
@@ -378,31 +378,33 @@
warning("The MOP test is defined for models with only one endogenous variable")
return(NA)
}
- Z <- model.matrix(object, "instrument")
- X <- model.matrix(object)
- X2 <- X[, spec$isEndo, drop=FALSE]
- X1 <- X[, !spec$isEndo, drop=FALSE]
- if (ncol(X1))
+ Z2 <- model.matrix(object, "excludedExo")
+ X1 <- model.matrix(object, "includedExo")
+ X2 <- model.matrix(object, "includedEndo")
+ y <- modelResponse(object)
+ if (!is.null(X1))
{
- z2n <- !(colnames(Z) %in% colnames(X1))
- fitX1 <- lm.fit(X1, Z[,z2n])
- Z <- as.matrix(fitX1$residuals)
+ fitX1 <- lm.fit(X1, Z2)
+ Z2 <- fitX1$residuals
X2 <- qr.resid(fitX1$qr, X2)
+ y <- qr.resid(fitX1$qr, y)
}
- Z <- qr.Q(qr(Z))
- colnames(Z) <- paste("Z", 1:ncol(Z), sep="")
- g <- reformulate(colnames(Z), colnames(X2), FALSE)
- h <- reformulate(colnames(Z), NULL, FALSE)
- dat <- data.frame(cbind(X2, Z))
+ 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 <- lm.fit(Z, X2)$coefficients
+ Pi <- c(crossprod(Z2,X2))/nrow(Z2)
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/spec$n
+ ## 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
Keff <- se^2*(1+2*x)/(se2+2*x*se*me)
crit <- qchisq(1-size, Keff, Keff*x)/Keff
@@ -437,28 +439,32 @@
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)
+ if (!is.null(X1))
+ {
+ 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))*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)]
+ w2 <- w[(ncol(Z2)+1):ncol(w), (ncol(Z2)+1):ncol(w)]
+ w12 <- w[(ncol(Z2)+1):ncol(w), 1:ncol(Z2)]
+ list(w1=w1,w2=w2,w12=w12,omega=omega)
}
Modified: pkg/momentfit/man/vcovHAC-methods.Rd
===================================================================
--- pkg/momentfit/man/vcovHAC-methods.Rd 2024-02-06 15:02:14 UTC (rev 221)
+++ pkg/momentfit/man/vcovHAC-methods.Rd 2024-02-07 17:01:28 UTC (rev 222)
@@ -2,9 +2,10 @@
\docType{methods}
\alias{vcovHAC-methods}
\alias{vcovHAC,momentModel-method}
+\alias{vcovHAC,sysModel-method}
\title{ ~~ Methods for Function \code{vcovHAC} in Package \pkg{sandwich} ~~}
\description{
- Methods to compute the HAC covariance matrix of the moment matrix ~~
+ Methods to compute the HAC covariance matrix of the moment objects ~~
}
\section{Methods}{
\describe{
@@ -11,6 +12,10 @@
\item{\code{signature(x = "momentModel")}}{
}
+
+\item{\code{signature(x = "sysModel")}}{
+}
+
}}
\examples{
Modified: pkg/momentfit/vignettes/weak.Rmd
===================================================================
--- pkg/momentfit/vignettes/weak.Rmd 2024-02-06 15:02:14 UTC (rev 221)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-02-07 17:01:28 UTC (rev 222)
@@ -15,6 +15,7 @@
%\VignetteEncoding{UTF-8}
header-includes:
- \newcommand{\E}{\mathrm{E}}
+- \newcommand{\eval}{\mathrm{eval}}
- \newcommand{\diag}{\mathrm{diag}}
- \newcommand{\Prob}{\mathrm{Pr}}
- \newcommand{\Var}{\mathrm{Var}}
@@ -21,6 +22,7 @@
- \newcommand{\Vect}{\mathrm{Vec}}
- \newcommand{\Cov}{\mathrm{Cov}}
- \newcommand{\Cor}{\mathrm{Cor}}
+- \newcommand{\tr}{\mathrm{tr}}
- \newcommand{\conP}{\overset{p}{\to}}
- \newcommand{\conD}{\overset{d}{\to}}
- \newcommand\R{ \mathbb{R} }
@@ -546,13 +548,16 @@
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 (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\,.$$
+ normalize the matrix $Z_2$ so that $Z_2'Z_2/n=I$ (we replace $Z$ by
+ $Q$ from its QR decomposition times $\sqrt{n}$) . 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$.
-- The test is $F_{eff} = (X_2'Z_2Z_2'X_2)/[n^2tr(\hat{W}_2)]$
+- The test is $F_{eff} = (X_2'Z_2Z_2'X_2)/[n\tr(\hat{W}_2)]$, where
+ $\tr(A)$ is the trace of $A$. Since
+ $\hat{\Pi}_2=(Z_2'Z_2)^{-1}Z_2'X_2=Z_2'X_2/n$, we can write the
+ test as $F_{eff} = n\hat{\Pi}_2'\hat{\Pi}_2/\tr(\hat{W}_2)$.
This is computed by the function `MOPtest`. For now, no critical
values are reported. Will be added soon.
@@ -595,8 +600,7 @@
\end{eqnarray*}
-For example, we can compute $W_2$ above
-as follows for `mod2`.
+For example, we can compute $W_2$ above as follows for `mod2`.
- We extract $Z_2$, $X_1$, $X_2$ and $y$,
@@ -610,7 +614,7 @@
```
- We project $X_1$ off $Z_2$, $X_2$ and $y$ and normalize $Z$ using
- its QR decomposition:
+ its QR decomposition times $\sqrt{n}$:
```{r}
fitX1 <- lm.fit(X1, Z2)
@@ -617,13 +621,15 @@
Z2 <- fitX1$residuals
X2 <- qr.resid(fitX1$qr, X2)
y <- qr.resid(fitX1$qr, y)
-Z2 <- qr.Q(qr(Z2))
+Z2 <- qr.Q(qr(Z2))*sqrt(nrow(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$:
+package. We just need to create a `linearModel` object with no
+endogenous variables. For $\hat{W}_2$, we regress $X_2$ on $Z_2$ and
+use $Z_2$ as instruments. We can set the argument `vcov` to `"MDS"` to
+obtain a robust to heteroskedasticity $\hat{W}_2$ (or to `"CL"` for
+clustered or `"HAC"` for serially correlated errors).
```{r}
colnames(Z2) = paste("Z", 1:ncol(Z2), sep="")
@@ -638,71 +644,97 @@
$Z_2'e/\sqrt{n}$:
```{r}
-b <- lm.fit(Z2,X2)$coef
+b <- crossprod(Z2,X2)/nrow(Z2)
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_{1}$, which is the asymptotic variance of $Z_2'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:
+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:
-\[
-\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
+dat <- as.data.frame(cbind(y=y,X2,Z2))
+g <- list(reformulate(colnames(Z2), "y", FALSE),
+ g)
+h <- list(h,h)
+m <- sysMomentModel(g,h,data=dat,vcov="MDS")
+b <- list(c(crossprod(Z2,y)/nrow(Z2)),
+ c(crossprod(Z2,X2)/nrow(Z2)))
+w <- vcov(m, b)
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$:
+Note that we need to adjust the sample size. The way the model `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.
+This is what the function `getMOPW` 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. Here we
+see $\hat{W}_2$:
-
```{r}
res <- momentfit:::getMOPW(mod2)
-res
+res$w2
```
+The function also returns the elements `w1`, `w12` and `omega`. The
+latter is $\hat\Omega=[\hat v, \hat e]'[\hat v, \hat e]/n$. The
+matrices $\hat W_1$ and $\hat W_{12}$ are needed to compute the
+effective degrees of freedom:
+\[
+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.
+
+```{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")
+```
+
+
+
+
## 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