[Gmm-commits] r218 - in pkg/momentfit: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 30 17:16:13 CET 2024
Author: chaussep
Date: 2024-01-30 17:16:12 +0100 (Tue, 30 Jan 2024)
New Revision: 218
Added:
pkg/momentfit/R/sysdata.rda
pkg/momentfit/man/CDtest.Rd
Modified:
pkg/momentfit/DESCRIPTION
pkg/momentfit/NAMESPACE
pkg/momentfit/R/weak.R
pkg/momentfit/man/kclassfit.Rd
pkg/momentfit/vignettes/weak.Rmd
pkg/momentfit/vignettes/weak.pdf
Log:
added the Cragg-Donald test and Stock-Yogo critical values
Modified: pkg/momentfit/DESCRIPTION
===================================================================
--- pkg/momentfit/DESCRIPTION 2024-01-25 22:40:48 UTC (rev 217)
+++ pkg/momentfit/DESCRIPTION 2024-01-30 16:16:12 UTC (rev 218)
@@ -5,7 +5,7 @@
Author: Pierre Chausse <pchausse at uwaterloo.ca>
Maintainer: Pierre Chausse <pchausse at uwaterloo.ca>
Description: Several classes for moment-based models are defined. The classes are defined for moment conditions derived from a single equation or a system of equations. The conditions can also be expressed as functions or formulas. Several methods are also offered to facilitate the development of different estimation techniques. The methods that are currently provided are the Generalized method of moments (Hansen 1982; <doi:10.2307/1912775>), for single equations and systems of equation, and the Generalized Empirical Likelihood (Smith 1997; <doi:10.1111/j.0013-0133.1997.174.x>, Kitamura 1997; <doi:10.1214/aos/1069362388>, Newey and Smith 2004; <doi:10.1111/j.1468-0262.2004.00482.x>, and Anatolyev 2005 <doi:10.1111/j.1468-0262.2005.00601.x>).Some work is being done to add tools to deal with weak and/or many instruments. This include K-Class estimators (LIML and Fuller), Anderson and Rubin statistics test, etc.
-Depends: R (>= 3.0.0), sandwich
+Depends: R (>= 3.5), sandwich
Imports: stats, methods, parallel
Suggests: lmtest, knitr, texreg, rmarkdown, ivmodel
Collate: 'allClasses.R' 'validity.R' 'momentData.R' 'momentModel-methods.R'
Modified: pkg/momentfit/NAMESPACE
===================================================================
--- pkg/momentfit/NAMESPACE 2024-01-25 22:40:48 UTC (rev 217)
+++ pkg/momentfit/NAMESPACE 2024-01-30 16:16:12 UTC (rev 218)
@@ -43,7 +43,8 @@
printRestrict, restModel, getRestrict, gmm4, sysMomentModel, ThreeSLS,
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)
+ evalGelObj, momFct, gel4, setCoef, lse, getK, kclassfit, CDtest,
+ SYTables)
### S3 methods ###
Added: pkg/momentfit/R/sysdata.rda
===================================================================
(Binary files differ)
Index: pkg/momentfit/R/sysdata.rda
===================================================================
--- pkg/momentfit/R/sysdata.rda 2024-01-25 22:40:48 UTC (rev 217)
+++ pkg/momentfit/R/sysdata.rda 2024-01-30 16:16:12 UTC (rev 218)
Property changes on: pkg/momentfit/R/sysdata.rda
___________________________________________________________________
Added: svn:mime-type
## -0,0 +1 ##
+application/octet-stream
\ No newline at end of property
Modified: pkg/momentfit/R/weak.R
===================================================================
--- pkg/momentfit/R/weak.R 2024-01-25 22:40:48 UTC (rev 217)
+++ pkg/momentfit/R/weak.R 2024-01-30 16:16:12 UTC (rev 218)
@@ -62,15 +62,23 @@
}
-kclassfit <- function(object, k, type=c("LIML", "Fuller"), alpha = 1)
+kclassfit <- function(object, k, type=c("LIML", "Fuller", "BTSLS"), alpha = 1)
{
if (!inherits(object, "linearModel"))
stop("object must be a model of class linearModel")
type <- match.arg(type)
+ spec <- modelDims(object)
if (missing(k))
{
- Kres <- getK(object, alpha, TRUE)
- k <- Kres$k[type]
+ if (type == "BTSLS")
+ {
+ Kres <- list()
+ k2 <- spec$q-sum(!spec$isEndo)
+ k <- spec$n/(spec$n-k2+2)
+ } else {
+ Kres <- getK(object, alpha, TRUE)
+ k <- Kres$k[type]
+ }
method <- type
} else {
method="K-Class"
@@ -93,7 +101,6 @@
return(tsls(object))
if (k==0)
return(lse(object))
- spec <- modelDims(object)
EndoVars <- !(spec$parNames %in% spec$momNames)
exoInst <- spec$momNames %in% spec$parNames
if (all(!EndoVars))
@@ -182,3 +189,91 @@
invisible()
})
setMethod("show","summaryKclass", function(object) print(object))
+
+
+## Stock and Yogo (2005)
+
+CDtest <- function(object, print=TRUE, ...)
+{
+ if (!inherits(object, "linearModel"))
+ 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]
+ secS <- lm.fit(Z, X2)
+ Z2 <- Z[,z2n,drop=FALSE]
+ tmp <- if (ncol(X2)>1)
+ Z[,z2n,drop=FALSE]%*%secS$coef[z2n,,drop=FALSE]
+ else
+ Z[,z2n,drop=FALSE]%*%secS$coef[z2n]
+ e <- secS$residuals
+ e2 <- lm.fit(X1, X2)$residuals
+ e <- crossprod(e)/(spec$n-spec$q)
+ e2 <- crossprod(e2, tmp)/ncol(Z2)
+ test <- min(eigen(solve(e,e2))$val)
+ if (!print)
+ return(test)
+ cat("Cragg and Donald Test for Weak Instruments\n",
+ "******************************************\n", sep="")
+ cat("Number of included Endogenous: ", ncol(X2), "\n", sep="")
+ cat("Number of excluded Exogenous: ", ncol(Z2), "\n", sep="")
+ cat("Statistics: ", formatC(test, ...), "\n\n", sep="")
+ SYTables(object)
+}
+
+SYTables <- function(object, print=TRUE)
+{
+ if (!inherits(object, "linearModel"))
+ stop("object must be of class linearModel")
+ s <- modelDims(object)
+ l <- sum(s$isEndo)
+ k <- s$q - (s$k - l)
+ t <- sizeTSLS
+ sizeTSLS <- t[attr(t, "excExo") == k, attr(t, "incEndo") == l]
+ names(sizeTSLS) <- paste("size=",
+ attr(t, "size")[attr(t, "incEndo") == l], sep="")
+ t <- biasFuller
+ biasFuller <- t[attr(t, "excExo") == k, attr(t, "incEndo") == l]
+ names(biasFuller) <- paste("bias=",
+ attr(t, "bias")[attr(t, "incEndo") == l], sep="")
+ t <- sizeLIML
+ sizeLIML <- t[attr(t, "excExo") == k, attr(t, "incEndo") == l]
+ names(sizeLIML) <- paste("size=",
+ attr(t, "size")[attr(t, "incEndo") == l], sep="")
+ if ((k-l)>=2)
+ {
+ t <- biasTSLS
+ biasTSLS <- t[attr(t, "excExo") == k, attr(t, "incEndo") == l]
+ names(biasTSLS) <- paste("bias=",
+ attr(t, "bias")[attr(t, "incEndo") == l], sep="")
+ } else {
+ biasTSLS <- NULL
+ }
+
+ crit <- list(biasTSLS=biasTSLS,
+ sizeTSLS=sizeTSLS,
+ biasFuller=biasFuller,
+ sizeLIML=sizeLIML)
+ if (!print)
+ return(crit)
+ nCrit <- c("Target relative bias for TSLS:\n",
+ "Target size for TSLS:\n",
+ "Target relative bias for Fuller-K:\n",
+ "Target size for LIML:\n")
+
+ cat("Stock and Yogo (2005) critical values\n")
+ cat("*************************************\n")
+ for (i in 1:4)
+ {
+ if (!is.null(crit[[i]]))
+ {
+ cat(nCrit[i])
+ print.default(format(crit[[i]]), print.gap = 2L, quote = FALSE)
+ cat("\n")
+ }
+ }
+ invisible()
+}
Added: pkg/momentfit/man/CDtest.Rd
===================================================================
--- pkg/momentfit/man/CDtest.Rd (rev 0)
+++ pkg/momentfit/man/CDtest.Rd 2024-01-30 16:16:12 UTC (rev 218)
@@ -0,0 +1,73 @@
+\name{CDtest}
+\alias{CDtest}
+\alias{SYTables}
+\title{
+Cragg and Donald test for Weak Instruments
+}
+\description{
+The Cragg and Donald test is based on the F-test from the first stage
+regression. For models with one endogenous variable, it is the usual
+F-test. For more than one endogenous variables, it is the smallest
+eigenvalue of the matrix representation of the F statistic. Along with
+the statistics, the critical values of Stock and Yogo (2005) for the
+null hypothesis of weak instruments are provided.
+}
+\usage{
+
+CDtest(object, print=TRUE, ...)
+
+SYTables(object, print=TRUE)
+
+}
+
+\arguments{
+ \item{object}{
+ A model of class \code{linearModel}
+ }
+ \item{print}{
+ If code{TRUE}, the result is printed and nothing is returned. See
+ below for details on what the functions return when it is set to
+ \code{FALSE}
+ }
+ \item{...}{Arguments passed to \code{\link{formatC}} to print the statistic.
+ }
+}
+
+\details{ The \code{CDtest} function computes the test for the model
+ object and runs \code{SYTables} to extract the critical values that are
+ specific to the model. }
+
+\value{
+ The function \code{CDtest} returns the Cragg anb Donald statistic when
+ \code{print} is set to \code{FALSE}.
+
+ The function \code{SYTable} returns a list with the following elements
+ when \code{print} is set to \code{FALSE}:
+
+ \item{biasTSLS, biasFuller}{Named vectors of critical values for TSLS
+ or Fuller (see \code{\link{kclassfit}}). These critical values are
+ used to achieve the maximum relative bias given by their names. The
+ values are defined only for models with a number of overidentifying
+ restrictions greater or equal to 2.}
+
+ \item{sizeTSLS, sizeLIML}{A named vector of critical values for TSLS
+ or LIML (see \code{\link{kclassfit}}). These critical values are used
+ to acheive a size that does not exceed the one given by their names.}
+
+ Both functions return \code{NULL} invisibly when \code{print} is set to
+ \code{TRUE}
+
+}
+\examples{
+
+data(simData)
+theta <- c(beta0=1,beta1=2)
+model1 <- momentModel(y~x1, ~z1+z2, data=simData)
+CDtest(model1)
+}
+
+\keyword{LIML}
+\keyword{Fuller}
+\keyword{Weak Instruments Test}
+\keyword{Two Stage Least Squares}
+
Modified: pkg/momentfit/man/kclassfit.Rd
===================================================================
--- pkg/momentfit/man/kclassfit.Rd 2024-01-25 22:40:48 UTC (rev 217)
+++ pkg/momentfit/man/kclassfit.Rd 2024-01-30 16:16:12 UTC (rev 218)
@@ -9,7 +9,7 @@
includes LIML, Fuller, TSLS and OLS.
}
\usage{
-kclassfit(object, k, type = c("LIML", "Fuller"), alpha = 1)
+kclassfit(object, k, type = c("LIML", "Fuller", "BTSLS"), alpha = 1)
getK(object, alpha = 1, returnRes = FALSE)
}
@@ -24,7 +24,7 @@
}
\item{type}{
Which value of \code{k} should we use to fit the model? Only used if
- \code{k} is missing.
+ \code{k} is missing.
}
\item{alpha}{A parameter for the Fuller method}
@@ -36,7 +36,9 @@
Let the model be \eqn{Y=X\beta+U}{Y=XB+U} and the matrix of instruments
be \eqn{Z}{Z}. The K-Class estimator is
\eqn{(X'(I-kM_z)X)^{-1}(X'(I-kM_z)Y)}. The function \code{getK} can be
- used to compute the value of \code{k} for both LIML and Fuller.
+ used to compute the value of \code{k} for both LIML and Fuller. When
+ \code{type="BTSLS"}, the bias-adjusted TSLS of Nagar (1959) is
+ computed.
}
\value{
Modified: pkg/momentfit/vignettes/weak.Rmd
===================================================================
--- pkg/momentfit/vignettes/weak.Rmd 2024-01-25 22:40:48 UTC (rev 217)
+++ pkg/momentfit/vignettes/weak.Rmd 2024-01-30 16:16:12 UTC (rev 218)
@@ -174,8 +174,9 @@
is equivalent to finding the smallest eigenvalue of
$(Y'M_zY)^{-1}Y'M_1Y$. An alternative to the LIML method was proposed
by @fuller77. The method is also a K-Class method with
-$\kappa_{ful}=\kappa_{liml}-\alpha/(n-q)$. The Fuller method happens
-to have better properties than LIML.
+$\kappa_{ful}=\kappa_{liml}-\alpha/(n-q)$, where $\alpha$ is
+parameter. It usually set to=1. The Fuller method happens to have
+better properties than LIML.
## Computing $\hat\kappa$
@@ -267,13 +268,14 @@
The function that computes the K-Class estimator is `kclassfit`. The
arguments are: `object`, the model object, `k`, the value of $\kappa$,
-`type`, the type of $\kappa$ to compute when `k` is missing (`"LIML"`
-or `"Fuller"`) and `alpha`, the parameter of the Fuller method (the
-default is 1). Note first that the estimator is a TSLS estimator when
-`k=1` and a LSE when it is equal to 0. The package already has a
-`tsls` method for `linearModel` objects, which is what `kclassfit`
-calls when `k=1`. For the LSE, a new method was created to facilitate
-the estimation of model objects by least squares. The method is `lse`:
+`type`, the type of $\kappa$ to compute when `k` is missing (`"LIML"`,
+`"Fuller"` or `"BTSLS"` for the biased corrected TSLE) and `alpha`,
+the parameter of the Fuller method (the default is 1). Note first that
+the estimator is a TSLS estimator when `k=1` and a LSE when it is
+equal to 0. The package already has a `tsls` method for `linearModel`
+objects, which is what `kclassfit` calls when `k=1`. For the LSE, a
+new method was created to facilitate the estimation of model objects
+by least squares. The method is `lse`:
```{r}
lse(mod2)
@@ -282,18 +284,22 @@
It is an object of class `lsefit` that contains the `lm` object from
the estimation. Therefore, the `kclassfit` function returns an object
of class `lsefit` when `k=0` and `tlsl` when `k=1`. For any other
-value, which includes LIML and Fuller, the function returns an object
-of class `kclassfit`. The object contains a `gmmfit` object, generated
-by the estimation of the artificially created just-identified model,
-the name of the method, the value of $\kappa$ and the original model.
+value, which includes LIML, Fuller and BTSLS ($\kappa=n/(n-l_2+2)$),
+the function returns an object of class `kclassfit`. The object
+contains a `gmmfit` object, generated by the estimation of the
+artificially created just-identified model, the name of the method,
+the value of $\kappa$ and the original model.
```{r}
(liml <- kclassfit(mod2))
(fuller <- kclassfit(mod2, type="Fuller"))
+(btsls <- kclassfit(mod2, type="BTSLS"))
```
-We see that the LIML and Fuller estimates I get are identical to the ones from the
-`ivmodel` package.
+Note that the biased-adjusted TSLS is just TSLS because the adjustment
+only affects the method when the number of excluded instruments is not
+equal to 2. We see in the following that the LIML and Fuller estimates
+I get are identical to the ones from the `ivmodel` package.
```{r}
print(mod$LIML$point.est,digits=10)
@@ -377,11 +383,49 @@
you should get identical point estimate and both standard errors are
equal to 0.0576098.
-# Weak Instruments (For later use)
+# Weak Instruments
+## Testing for weak instrument: Stock-Yogo (2005)
-## Data Generating Process
+This test and the critical values are for model with homoskedastic
+errors. The test is the smallest eigenvalue of the following:
+\[
+\hat\Sigma_e^{-1/2}
+\Big[
+X_2'M_1Z_2(Z_2'M_1Z_2)^{-1}Z_2'M_1X_2
+\Big]
+\hat\Sigma_e^{-1/2} = [M_1 X_2]'[Z_2\hat\Pi_2]
+\]
+
+where $\hat\Sigma_e = \hat{e}'\hat{e}/(n-l_2-k_1)$. The function
+`CDtest`, which stands for Cragg and Donald test, computes this
+statistic. We can see that it returns the F-statistics when the number
+of endogenous variables is equal to 1:
+
+```{r}
+CDtest(mod2, print=FALSE)
+momentStrength(mod2)
+```
+
+Here is what we obtain for `mod3`
+
+```{r}
+CDtest(mod3)
+```
+
+```{r}
+g <- reformulate(c("educ", Xname), "lwage")
+h <- reformulate(c(c("nearc4","nearc2","IQ","KWW"), Xname))
+mod5 <- momentModel(g, h, data=card.data, vcov="iid")
+CDtest(mod5)
+```
+
+
+
+
+## Data Generating Process (for later use)
+
The following function is used to generate dataset with $k$
instruments and different level of strength. The DGP is
@@ -439,3 +483,5 @@
# References
+
+
Modified: pkg/momentfit/vignettes/weak.pdf
===================================================================
(Binary files differ)
More information about the Gmm-commits
mailing list