[Gmm-commits] r137 - in pkg/gmm4: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 28 21:51:33 CEST 2018
Author: chaussep
Date: 2018-09-28 21:51:33 +0200 (Fri, 28 Sep 2018)
New Revision: 137
Modified:
pkg/gmm4/NAMESPACE
pkg/gmm4/R/gmmModel.R
pkg/gmm4/R/gmmModels-methods.R
pkg/gmm4/R/gmmfit-methods.R
pkg/gmm4/R/validity.R
pkg/gmm4/man/momentStrength-methods.Rd
pkg/gmm4/vignettes/empir.bib
pkg/gmm4/vignettes/gmmS4.Rnw
pkg/gmm4/vignettes/gmmS4.pdf
Log:
added clustered weighting matrix using meatCL of the sandwich package
Modified: pkg/gmm4/NAMESPACE
===================================================================
--- pkg/gmm4/NAMESPACE 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/NAMESPACE 2018-09-28 19:51:33 UTC (rev 137)
@@ -10,7 +10,7 @@
"fitted", "lm.fit", "pchisq", "pnorm", "printCoefmat", "anova",
"model.frame", "reformulate", "formula", "nlminb", "kernapply",
"constrOptim", "kernel")
-importFrom("sandwich", "vcovHAC", "estfun","kernHAC",
+importFrom("sandwich", "vcovHAC", "estfun","kernHAC","vcovCL", "meatCL",
"bread","bwAndrews","bwNeweyWest","weightsAndrews",
"weightsLumley", "vcovHC")
### S4 Methods and Classes
Modified: pkg/gmm4/R/gmmModel.R
===================================================================
--- pkg/gmm4/R/gmmModel.R 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/R/gmmModel.R 2018-09-28 19:51:33 UTC (rev 137)
@@ -26,7 +26,7 @@
names(option$bw) <- "Fixed"
} else if (type=="CL") {
option <- list(cluster=NULL, type="HC0", cadjust=TRUE,
- milti0=FALSE)
+ multi0=FALSE)
if (length(addO) > 0)
{
if (!all(names(addO) %in% names(option)))
Modified: pkg/gmm4/R/gmmModels-methods.R
===================================================================
--- pkg/gmm4/R/gmmModels-methods.R 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/R/gmmModels-methods.R 2018-09-28 19:51:33 UTC (rev 137)
@@ -31,8 +31,13 @@
cat("Number of moment conditions: ", d$q, "\n", sep="")
if (!inherits(x, "functionGmm"))
cat("Number of Endogenous Variables: ", sum(x at isEndo), "\n", sep="")
- cat("Sample size: ", d$n, "\n")})
+ if (!is.null(x at survOptions$weights) && x at survOptions$type == "frequency")
+ cat("Implied sample size (sum of weights): ", d$n, "\n")
+ else
+ cat("Sample size: ", d$n, "\n")})
+
+
setMethod("show", "gmmModels", function(object) print(object))
##### coef ########
@@ -405,9 +410,14 @@
sig <- sd(residuals(object, theta))
Z <- model.matrix(object, "instrument")
w <- sig^2*crossprod(Z)/nrow(Z)
+ } else if (object at vcov == "CL") {
+ gt <- evalMoment(object, theta)
+ class(gt) <- "gmmFct"
+ opt <- object at vcovOptions
+ opt$x <- gt
+ w <- do.call(meatCL, opt)
} else {
w <- vcovHAC(object, theta)
- attr(w, "inv") <- TRUE
}
w})
@@ -448,6 +458,13 @@
if (w$rank < object at q)
warning("The moment matrix is not full column rank")
type <- "qr"
+ } else if (object at vcov == "CL") {
+ gt <- evalMoment(object, theta)
+ class(gt) <- "gmmFct"
+ opt <- object at vcovOptions
+ opt$x <- gt
+ w <- chol(do.call(meatCL, opt))
+ type <- "chol"
} else {
w <- vcovHAC(object, theta)
wSpec <- attr(w,"Spec")
@@ -558,7 +575,7 @@
})
setMethod("momentStrength", signature("linearGmm"),
- function(object, theta, vcovType=c("OLS","HC","HAC")){
+ function(object, theta, vcovType=c("OLS","HC","HAC","CL")){
spec <- modelDims(object)
EndoVars <- !(spec$parNames %in% spec$momNames)
exoInst <- spec$momNames %in% spec$parNames
@@ -577,9 +594,11 @@
{
resu <- lm(X[,i ]~Z-1)
v <- switch(vcovType,
- OLS=vcov(resu),
- HC=vcovHC(resu,"HC1"),
- HAC=vcovHAC(resu))[!exoInst,!exoInst]
+ OLS=vcov(resu),
+ HC=vcovHC(resu,"HC1"),
+ HAC=vcovHAC(resu),
+ CL=do.call(vcovCL,c(object at vcovOptions, x=resu)))
+ v <- v[!exoInst,!exoInst]
b <- coef(resu)[!exoInst]
f <- b%*%solve(v, b)/df1
df2 <- resu$df.residual
@@ -676,12 +695,9 @@
x at modelF <- x at modelF[i,,drop=FALSE]
x at instF <- x at instF[i,,drop=FALSE]
if (!is.null(x at vcovOptions$cluster))
- {
- if (!is.null(dim(x at vcovOptions$cluster)))
- x at vcovOptions$cluster <- x at vcovOptions$cluster[i]
- else
- x at vcovOptions$cluster <- x at vcovOptions$cluster[i,,drop=FALSE]
- }
+ x at vcovOptions$cluster <- x at vcovOptions$cluster[i,,drop=FALSE]
+ if (!is.null(x at survOptions$weights))
+ x at survOptions$weights <- x at survOptions$weights[i]
x at n <- nrow(x at modelF)
x})
@@ -694,12 +710,9 @@
else
stop("X is not subsetable")
if (!is.null(x at vcovOptions$cluster))
- {
- if (!is.null(dim(x at vcovOptions$cluster)))
- x at vcovOptions$cluster <- x at vcovOptions$cluster[i]
- else
- x at vcovOptions$cluster <- x at vcovOptions$cluster[i,,drop=FALSE]
- }
+ x at vcovOptions$cluster <- x at vcovOptions$cluster[i,,drop=FALSE]
+ if (!is.null(x at survOptions$weights))
+ x at survOptions$weights <- x at survOptions$weights[i]
x at n <- nrow(x at X)
x})
@@ -707,12 +720,9 @@
function(x, i) {
x at modelF <- x at modelF[i,,drop=FALSE]
if (!is.null(x at vcovOptions$cluster))
- {
- if (!is.null(dim(x at vcovOptions$cluster)))
- x at vcovOptions$cluster <- x at vcovOptions$cluster[i]
- else
- x at vcovOptions$cluster <- x at vcovOptions$cluster[i,,drop=FALSE]
- }
+ x at vcovOptions$cluster <- x at vcovOptions$cluster[i,,drop=FALSE]
+ if (!is.null(x at survOptions$weights))
+ x at survOptions$weights <- x at survOptions$weights[i]
x at n <- nrow(x at modelF)
x})
Modified: pkg/gmm4/R/gmmfit-methods.R
===================================================================
--- pkg/gmm4/R/gmmfit-methods.R 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/R/gmmfit-methods.R 2018-09-28 19:51:33 UTC (rev 137)
@@ -187,7 +187,8 @@
vcovType <- switch(object at model@vcov,
HAC="HAC",
iid="OLS",
- MDS="HC")
+ MDS="HC",
+ CL="CL")
strength <- momentStrength(object at model, coef(object), vcovType)
dimnames(coef) <- list(names(par),
c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
Modified: pkg/gmm4/R/validity.R
===================================================================
--- pkg/gmm4/R/validity.R 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/R/validity.R 2018-09-28 19:51:33 UTC (rev 137)
@@ -1,11 +1,12 @@
############# Validity Methods ##################
#### All validity Methods for all objects are here
-.checkBased <- function(object)
+.checkBased <- function(object, n)
{
error <- character()
vcov <- c("HAC","MDS","iid","CL")
- chk <- try(.getVcovOptions(object at vcov, object at vcovOptions), silent=TRUE)
+ chk <- try(.getVcovOptions(object at vcov, NULL, object at vcovOptions), silent=TRUE)
+ chk2 <- try(.getSurvOptions(NULL, object at survOptions), silent=TRUE)
if (class(chk) == "try-error")
{
msg <- "Invalid vcovOptions"
@@ -17,23 +18,62 @@
error <- c(error, msg)
}
}
+ if (class(chk2) == "try-error")
+ {
+ msg <- "Invalid survOptions"
+ error <- c(error, msg)
+ } else {
+ if (!isTRUE(all.equal(chk2, object at survOptions)))
+ {
+ msg <- "Invalid vcovOptions"
+ error <- c(error, msg)
+ }
+ }
if ( !(object at vcov%in%vcov))
{
vcov <- paste(vcov, collapse=", ")
msg <- paste("vcov must be one of ",
vcov, sep="")
error <- c(error, msg)
- }
+ }
+ if (object at vcov == "CL")
+ {
+ if (!is.data.frame(object at vcovOptions$cluster))
+ {
+ msg <- "cluster must be a data.frame"
+ error <- c(error, msg)
+ } else {
+ if (nrow(object at vcovOptions$cluster) != n)
+ {
+ msg <- "Wrong number of observations in cluster"
+ error <- c(error, msg)
+ }
+ }
+ }
+ if (length(object at survOptions)>0)
+ {
+ if (!inherits(object at survOptions$weights, c("integer","numeric")))
+ {
+ msg <- "weights must be a numeric or integer vector"
+ error <- c(error, msg)
+ } else {
+ if (length(object at survOptions$weights) != n)
+ {
+ msg <- "Wrong number of observations in weights"
+ error <- c(error, msg)
+ }
+ }
+ }
}
.checkLinGmm <- function(object)
{
- error <- .checkBased(object)
ny <- NCOL(object at modelF[[1L]])
nM <- nrow(object at modelF)
nI <- nrow(object at instF)
k <- ncol(model.matrix(object))
q <- ncol(model.matrix(object, "instrument"))
+ error <- .checkBased(object, nM)
if (k != object at k)
{
msg <- "k does not correspond to the number of regressor in modelF"
@@ -121,17 +161,6 @@
msg <- "The Hypothesis matrix if not full rank"
error <- c(error, msg)
}
- #res <- try(.imposeRestrict(R,rhs, spec$originParNames), silent=TRUE)
- #if (any(class(res)=="try-error"))
- # {
- # msg <- "Fail to generate constraint specifications"
- # error <- c(error, msg)
- # }
- #if (!identical(res[-which(names(res)=="isEndo")], spec))
- # {
- # msg <- "cstSpec does not seem to correspond to cstLHS and cstRHS"
- # error <- c(error, msg)
- # }
if (length(error)==0)
TRUE
else
@@ -142,10 +171,10 @@
.checkNLinGmm <- function(object)
{
- error <- .checkBased(object)
nM <- nrow(object at modelF)
nI <- nrow(object at instF)
k <- length(object at theta0)
+ error <- .checkBased(object, nM)
if (k != object at k)
{
msg <- "k does not correspond to the number of resgressor in modelF"
@@ -206,34 +235,56 @@
setValidity("nonlinearGmm", .checkNLinGmm)
-.checkfGmm <- function(object)
+.checkformGmm <- function(object)
{
- error <- .checkBased(object)
+ n <- nrow(object at modelF)
+ error <- .checkBased(object, n)
k <- length(object at theta0)
if (k != object at k)
{
- msg <- "k does not correspond to the number of resgressor in modelF"
+ msg <- "k does not correspond to the number of regressor in modelF"
error <- c(error, msg)
}
- mom <- try(object at fct(object at theta0, object at X))
- if (!is.null(object at dfct))
+ varList <- c(as.list(object at theta0), as.list(object at modelF))
+ rhs <- sapply(object at fRHS, function(l)
+ class(try(eval(l, varList), silent=TRUE))=="try-error")
+ lhs <- sapply(object at fLHS, function(l)
+ class(try(eval(l, varList), silent=TRUE))=="try-error")
+ ql <- length(lhs)
+ qr <- length(rhs)
+ if (any(rhs))
{
- dmom <- try(object at dfct(object at theta0, object at X))
- if (any(class(dmom)=="try-error"))
- {
- msg <- paste("Cannot evaluate the provided derivatives of",
- "moments at theta0\n",
- attr(mom,"conditon"))
- error <- c(error, msg)
- } else {
- if (ncol(dmom) != object at k | nrow(dmom) != object at q)
- {
- msg <- paste("The dimention of the gradian of the ",
- "moments is not qxk",sep="")
- error <- c(error, msg)
- }
- }
+ msg <- paste("Some RHS's cannot be evaluated at theta0")
+ error <- c(error, msg)
}
+ if (any(lhs))
+ {
+ msg <- paste("Some LHS's cannot be evaluated at theta0")
+ error <- c(error, msg)
+ }
+ if (any(c(ql,qr) != length(object at momNames)))
+ {
+ msg <- "the length fRHS or fLHS does not match the length of momNames"
+ error <- c(error, msg)
+ }
+ if (k != length(object at parNames))
+ {
+ msg <- "the length of parNames is not equal to k"
+ error <- c(error, msg)
+ }
+ if (length(error)==0)
+ TRUE
+ else
+ error
+ }
+
+setValidity("formulaGmm", .checkformGmm)
+
+.checkfGmm <- function(object)
+ {
+ mom <- try(object at fct(object at theta0, object at X))
+ k <- length(object at theta0)
+ error <- character()
if (any(class(mom)=="try-error"))
{
msg <- paste("Cannot evaluate the moments at theta0\n",
@@ -242,6 +293,7 @@
} else {
q <- ncol(mom)
n <- nrow(mom)
+ error <- c(error, .checkBased(object, n))
if (q != object at q)
{
msg <- paste("q does not correspond to the number of ",
@@ -254,6 +306,29 @@
error <- c(error, msg)
}
}
+ if (k != object at k)
+ {
+ msg <- "k does not correspond to the number of resgressor in modelF"
+ error <- c(error, msg)
+ }
+ if (!is.null(object at dfct))
+ {
+ dmom <- try(object at dfct(object at theta0, object at X))
+ if (any(class(dmom)=="try-error"))
+ {
+ msg <- paste("Cannot evaluate the provided derivatives of",
+ "moments at theta0\n",
+ attr(mom,"conditon"))
+ error <- c(error, msg)
+ } else {
+ if (ncol(dmom) != object at k | nrow(dmom) != object at q)
+ {
+ msg <- paste("The dimention of the gradian of the ",
+ "moments is not qxk",sep="")
+ error <- c(error, msg)
+ }
+ }
+ }
if (q != length(object at momNames))
{
msg <- "the length of momNames is not equal to q"
Modified: pkg/gmm4/man/momentStrength-methods.Rd
===================================================================
--- pkg/gmm4/man/momentStrength-methods.Rd 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/man/momentStrength-methods.Rd 2018-09-28 19:51:33 UTC (rev 137)
@@ -13,7 +13,7 @@
}
\usage{
\S4method{momentStrength}{linearGmm}(object, theta,
-vcovType=c("OLS","HC","HAC"))
+vcovType=c("OLS","HC","HAC","CL"))
}
\arguments{
\item{object}{An object of class \code{"linearGmm"}}
@@ -25,7 +25,11 @@
the function \code{\link{vcovHC}} is used with "HC1", and
\code{\link{vcovHAC}} is used with the default setup is "HAC" is
chosen. In \code{summary} for \code{gmmfit} objects, it is adjusted
- to the type of covariance that is set in the object}
+ to the type of covariance that is set in the object. For type
+ \code{CL}, clustered covariance matrix is computed. The options are
+ the one included in the \code{vcovOptions} slot of the object (see
+ \code{\link{meatCL}}). The object must have be defined with clusters
+ for that to work. See \code{\link{gmmModel}}.}
}
\section{Methods}{
\describe{
Modified: pkg/gmm4/vignettes/empir.bib
===================================================================
--- pkg/gmm4/vignettes/empir.bib 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/vignettes/empir.bib 2018-09-28 19:51:33 UTC (rev 137)
@@ -201,6 +201,18 @@
url={http://www.jstatsoft.org/v16/i09/}
}
+ at article{berger-graham-zeileis17,
+ title = {Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in {R}},
+ author = {Berger, S. and Graham, N. and Zeileis, A.},
+ institution = {Working Papers in Economics and Statistics, Research Platform Empirical and Experimental Economics, Universit\"at Innsbruck},
+ year = {2017},
+ type = {Working Paper},
+ number = {2017-12},
+ month = {July},
+url = {http://EconPapers.RePEc.org/RePEc:inn:wpaper:2017-12}
+}
+ }
+
@article{chausse10,
title = {Computing Generalized Method of Moments and Generalized Empirical Likelihood with R},
author = {Pierre Chauss{\'e}},
Modified: pkg/gmm4/vignettes/gmmS4.Rnw
===================================================================
--- pkg/gmm4/vignettes/gmmS4.Rnw 2018-09-27 19:18:35 UTC (rev 136)
+++ pkg/gmm4/vignettes/gmmS4.Rnw 2018-09-28 19:51:33 UTC (rev 137)
@@ -93,7 +93,7 @@
\[
Y_i = X_i'\theta + \epsilon_i,
\]
-with the moment condition $\E[\epsilon_i(\theta)Z_i]=0$, where $X_i$ is $k\times 1$ and $Z_i$ is $q\times 1$ with $q\geq k$. We consider three possibilities for the asymptotic variance of $\sqrt{n}\bar{g}_i(\theta)$:
+with the moment condition $\E[\epsilon_i(\theta)Z_i]=0$, where $X_i$ is $k\times 1$ and $Z_i$ is $q\times 1$ with $q\geq k$. We consider four possibilities for the asymptotic variance of $\sqrt{n}\bar{g}_i(\theta)$:
\begin{enumerate}
\item[a)] ``iid'': Here we assume no autocorrelation and homoscedastic error with $\Var(\epsilon_i|Z_i)=\sigma^2$, which implies that the asymptotic variance $V$ is $\sigma^2\E[Z_iZ_i']$ and can be estimated by:
\[
@@ -110,6 +110,11 @@
\hat{V} = \sum_{i=-M}^{M} K_h(i)\hat{\Gamma}_i,
\]
where $K_h(i)$ is a kernel that depends on the bandwidth $h$, and $\hat{\Gamma}_i$ is an estimator of $\Gamma_i$.
+\item[d)] ``CL'': The sample is clustered. For one dimensional clusters, let $\bar{g}_i(\theta)$ be the sample mean of the moment function for cluster $i$ and $n_i$ be the number of observations in that cluster. Then, the clustered covariance matrix of the sample moment $\sqrt{n}\bar{g}(\theta)$ can be estimated as:
+ \[
+ \hat{V} = \frac{1}{n}\sum_{i=1}^{N_{cl}} n_i^2 \bar{g}_i(\hat{\theta}) \bar{g}_i'(\hat{\theta})
+ \]
+where $N_{cl}$ is the number of clusters. For higher dimensional clusters, like cities within provinces for example, we need to take into account that observations belong to more than one group. For a more detailed presentation with reference to recent developments, see \cite{berger-graham-zeileis17}.
\end{enumerate}
\item The nonlinear model:
\[
@@ -238,6 +243,42 @@
@
But the first approach may speed up estimation quite a bit for large dataset because it reduces the number of operations.
+Covariance matrix options can be modified using the argument ``vcovOptions''. For example, if we want a HAC matrix, several options such as the kernel and bandwidth can be modified. By default, the HAC is computed using the Quadratic Spectral kernel the optimal bandwidth of \cite{andrews91}. To modify the options, we proceed this way:
+
+<<>>=
+mod.hac <- gmmModel(y~x1+x2, ~x1+z2+z3, vcov="HAC",
+ vcovOptions=list(kernel="Bartlett", bw="NeweyWest"),
+ data=simData)
+mod.hac
+@
+
+See the help on ``vcovHAC'' of the sandwich package for more details on all possible parameters. For clustered covariance, we need to specify the clusters and some other options. Lets consider the following dataset:
+
+<<>>=
+data("InstInnovation", package = "sandwich")
+@
+
+We can use one-way clustering:
+
+<<>>=
+mod.cl1 <- gmmModel(sales~value, ~value, vcov="CL",
+ vcovOptions=list(cluster=~company),
+ data=InstInnovation)
+mod.cl1
+@
+
+or a two-way clustering:
+
+<<>>=
+mod.cl2 <- gmmModel(sales~value, ~value, vcov="CL",
+ vcovOptions=list(cluster=~company+year, multi0=TRUE),
+ data=InstInnovation)
+mod.cl2
+@
+
+The clustered covariance is computed using the ``meatCL'' function of the sandwich package. For more options, see its help file.
+
+
\subsection{Methods for gmmModels Classes} \label{sec:gmmmodels-methods}
\begin{itemize}
Modified: pkg/gmm4/vignettes/gmmS4.pdf
===================================================================
(Binary files differ)
More information about the Gmm-commits
mailing list