[Gmm-commits] r101 - in pkg/gmm: . R data man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 19 22:51:10 CEST 2017
Author: chaussep
Date: 2017-05-19 22:51:09 +0200 (Fri, 19 May 2017)
New Revision: 101
Added:
pkg/gmm/R/ategel.R
pkg/gmm/data/nsw.rda
pkg/gmm/man/ATEgel.Rd
pkg/gmm/man/marginal.Rd
pkg/gmm/man/nsw.Rd
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NAMESPACE
pkg/gmm/R/Methods.gel.R
pkg/gmm/R/gel.R
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/momentEstim.R
pkg/gmm/man/Finance.Rd
pkg/gmm/man/confint.Rd
pkg/gmm/man/getModel.Rd
pkg/gmm/man/summary.Rd
pkg/gmm/man/vcov.Rd
Log:
Added lots of stuff for ATE estimation
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/DESCRIPTION 2017-05-19 20:51:09 UTC (rev 101)
@@ -1,6 +1,6 @@
Package: gmm
Version: 1.6-1
-Date: 2016-05-17
+Date: 2017-05-19
Title: Generalized Method of Moments and Generalized Empirical
Likelihood
Author: Pierre Chausse <pchausse at uwaterloo.ca>
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/NAMESPACE 2017-05-19 20:51:09 UTC (rev 101)
@@ -3,6 +3,7 @@
importFrom(methods, is)
importFrom(graphics, abline, legend, lines, panel.smooth, par, plot, points)
importFrom(grDevices, dev.interactive, devAskNewPage, extendrange)
+importFrom(utils, tail)
export(gmm,summary.gmm,smoothG,getDat,summary.gel,getLamb,gel, estfun.gmmFct, estfun.gmm, estfun.gel, bread.gel, bread.gmm,
print.gmm,coef.gmm,vcov.gmm,print.summary.gmm, confint.gel, print.gel, print.summary.gel, vcov.gel, coef.gel, fitted.gmm,
@@ -14,12 +15,15 @@
KTest, print.gmmTests, gmmWithConst, estfun.tsls, model.matrix.tsls,vcov.tsls, bread.tsls, evalGmm, momentEstim.baseGmm.eval,
momentEstim.baseGel.eval, evalGel, confint.gmm, print.confint, sysGmm, getModel.sysGmm, momentEstim.sysGmm.twoStep.formula,
momentEstim.tsls.twoStep.formula, getModel.tsls, summary.sysGmm, print.sysGmm, print.summary.sysGmm,
- sur, threeSLS, randEffect, five, bwWilhelm)
-
+ sur, threeSLS, randEffect, five, bwWilhelm, getModel.ateGel, ATEgel,
+ summary.ategel, vcov.ategel, confint.ategel, marginal, marginal.ategel)
+
+S3method(marginal, ategel)
S3method(summary, gmm)
S3method(summary, sysGmm)
S3method(summary, tsls)
S3method(summary, gel)
+S3method(summary, ategel)
S3method(print, gmm)
S3method(print, sysGmm)
S3method(print, summary.gmm)
@@ -28,6 +32,7 @@
S3method(vcov, gmm)
S3method(vcov, tsls)
S3method(confint, gmm)
+S3method(confint, ategel)
S3method(fitted, gmm)
S3method(residuals, gmm)
S3method(plot, gmm)
@@ -38,6 +43,7 @@
S3method(print, gmmTests)
S3method(coef, gel)
S3method(vcov, gel)
+S3method(vcov, ategel)
S3method(confint, gel)
S3method(fitted, gel)
S3method(residuals, gel)
@@ -54,6 +60,7 @@
S3method(getModel, constGmm)
S3method(getModel, constGel)
S3method(getModel, tsls)
+S3method(getModel, ateGel)
S3method(momentEstim, baseGmm.twoStep)
S3method(momentEstim, baseGmm.eval)
S3method(momentEstim, baseGmm.twoStep.formula)
Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/R/Methods.gel.R 2017-05-19 20:51:09 UTC (rev 101)
@@ -299,3 +299,4 @@
+
Added: pkg/gmm/R/ategel.R
===================================================================
--- pkg/gmm/R/ategel.R (rev 0)
+++ pkg/gmm/R/ategel.R 2017-05-19 20:51:09 UTC (rev 101)
@@ -0,0 +1,294 @@
+ATEgel <- function(g, balm, y=NULL, treat=NULL, tet0=NULL,momType=c("bal","balSample","ATT"),
+ popMom = NULL, family=c("linear","logit", "probit"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD"),
+ tol_lam = 1e-9, tol_obj = 1e-9, tol_mom = 1e-9, maxiterlam = 100,
+ optfct = c("optim", "nlminb"),
+ optlam = c("nlminb", "optim", "iter"), data=NULL, Lambdacontrol = list(),
+ model = TRUE, X = FALSE, Y = FALSE, ...)
+ {
+ type <- match.arg(type)
+ optfct <- match.arg(optfct)
+ optlam <- match.arg(optlam)
+ momType <- match.arg(momType)
+ family <- match.arg(family)
+ TypeGel <- "ateGel"
+
+ all_args <- list(g = g, x = balm, y=y, treat=treat, tet0 = tet0, type = type,
+ tol_lam = tol_lam, tol_obj = tol_obj, tol_mom = tol_mom,
+ maxiterlam = maxiterlam, optfct = optfct,
+ optlam = optlam, model = model, X = X, Y = Y,
+ TypeGel = TypeGel, call = match.call(),
+ Lambdacontrol = Lambdacontrol, data = data,
+ constraint=FALSE, kernel = "Truncated", bw = bwAndrews,
+ approx = "AR(1)", prewhite = 1, ar.method = "ols",
+ tol_weights = 1e-7, alpha = NULL, eqConst = NULL,
+ eqConstFullVcov = FALSE, momType=momType, popMom=popMom,
+ family=family)
+ class(all_args)<-TypeGel
+ Model_info<-getModel(all_args)
+ z <- momentEstim(Model_info, ...)
+ class(z) <- c("ategel", "gel")
+ res <- .vcovate(z)
+ z$vcov_par <- res$vcov_par
+ z$vcov_lambda <- res$vcov_lambda
+ return(z)
+ }
+
+.momentFctATE <- function(tet, dat)
+ {
+ x <- dat$x
+ k <- attr(dat, "k")
+ ZT <- c(x[,2:(k+1),drop=FALSE]%*%tet[1:k])
+ if (is.null(attr(dat, "family")))
+ e <- x[,1] - ZT
+ else
+ e <- x[,1] - attr(dat, "family")$linkinv(ZT)
+ gt1 <- e * x[,2:(k+1)]
+ gt2 <- sweep(x[,3:(k+1),drop=FALSE],2,tet[(k+1):(2*k-1)],"-")
+ gt3 <- lapply(1:(k-1), function(i) gt2[,i]*x[,-(1:(k+1))])
+ gt3 <- do.call(cbind, gt3)
+ gt <- cbind(gt1,gt2,gt3)
+ if (is.null(attr(dat, "popMom")))
+ {
+ if (attr(dat, "momType") == "balSample")
+ {
+ momB <- scale(x[,-(1:(k+1))], scale=FALSE)
+ gt <- cbind(gt, momB)
+ }
+ if (attr(dat, "momType") == "ATT")
+ {
+ momB <- sweep(x[,-(1:3), drop=FALSE], 2,
+ colMeans(x[x[,3]==1,-(1:3), drop=FALSE]),
+ FUN="-")
+ gt <- cbind(gt, momB)
+ }
+ } else {
+ momB <- sweep(x[,-(1:(k+1)), drop=FALSE], 2, attr(dat, "popMom"), "-")
+ gt <- cbind(gt, momB)
+ }
+ return(as.matrix(gt))
+ }
+
+.DmomentFctATE <- function(tet, dat, pt=NULL)
+ {
+ if (is.null(pt))
+ pt <- rep(1/nrow(dat$x), nrow(dat$x))
+ x <- dat$x
+ k <- attr(dat, "k")
+ q <- dat$nh*(k-1)+2*k-1
+ Z <- x[,2:(k+1), drop=FALSE]
+ ZT <- c(Z%*%tet[1:k])
+ G <- matrix(0, q, 2*k-1)
+ if (is.null(attr(dat, "family")))
+ {
+ tau <- rep(1, nrow(x))
+ } else {
+ tau <- attr(dat, "family")$mu.eta(ZT)
+ }
+ G11 <- lapply(1:k, function(i) -colSums(pt*Z[,i]*tau*Z))
+ G[1:k, 1:k] <- do.call(rbind, G11)
+ G[(k+1):(2*k-1), (k+1):(2*k-1)] <- -sum(pt)*diag(k-1)
+ uK <- colSums(pt*x[,-(1:(k+1))])
+ G[(2*k):q, (k+1):(2*k-1)] <- -kronecker(diag(k-1), uK)
+ if (attr(dat, "momType") != "bal" | !is.null(attr(dat, "popMom")))
+ G <- rbind(G, matrix(0, dat$nh, 2*k-1))
+ return(G)
+ }
+
+
+.DmomentFctATE2 <- function(tet, dat, pt=NULL)
+ {
+ G <- .DmomentFctATE(tet, dat, pt)
+ k <- attr(dat, "k")
+ q <- dat$nh*(k-1)+2*k-1
+ if (is.null(pt))
+ pt <- rep(1/nrow(dat$x), nrow(dat$x))
+ if (attr(dat, "momType") != "bal" & is.null(attr(dat, "popMom")))
+ G <- cbind(G, rbind(matrix(0,q, dat$nh), -sum(pt)*diag(dat$nh)))
+ return(G)
+ }
+
+.psiGam <- function(object)
+ {
+ n <- nrow(object$dat$x)
+ nh <- object$dat$nh
+ lam <- object$lambda
+ q <- length(lam)
+ k <- attr(object$dat, "k")
+ theta <- object$coefficients
+ gt <- object$gt
+ rho1 <- .rho(x=gt, lamb=lam, derive=1, type=object$type)
+ rho2 <- .rho(x=gt, lamb=lam, derive=1, type=object$type)
+ Z <- object$dat$x[,2:(k+1)]
+ ZT <- c(Z%*%theta[1:k])
+ X <- object$dat$x[,-(1:(k+1))]
+ family <- attr(object$dat, "family")
+ momType <- attr(object$dat, "momType")
+ popMom <- attr(object$dat, "popMom")
+ if (is.null(family))
+ {
+ tau1 <- rep(1, n)
+ } else {
+ tau1 <- family$mu.eta(ZT)
+ tau2 <- family$mu.eta2(ZT, family)
+ }
+ lG1 <- sapply(1:k, function(i) -(tau1*Z[,i]*Z)%*%lam[1:k])
+ q2 <- nh*(k-1)+2*k-1
+ lamM <- matrix(lam[(2*k):q2], ncol=(k-1))
+ lG2 <- sapply(1:(k-1), function(i) -lam[k+i]-X%*%lamM[,i])
+ lG <- cbind(lG1, lG2)
+ G <- .DmomentFctATE2(theta, object$dat, rho1)
+ G22 <- crossprod(rho2*gt, gt)/n
+ if (momType == "bal" | !is.null(popMom))
+ {
+ Psi <- cbind(rho1*lG, rho1*gt)
+ G11 <- crossprod(rho2*lG, lG)/n
+ G12 <- t(G)/n + crossprod(rho2*lG, gt)/n
+ if (!is.null(family))
+ {
+ G12tmp <- lapply(1:k, function(i)
+ colSums(-rho1*tau2*Z[,i]*c(Z%*%lam[1:k])*Z))
+ G12.2 <- matrix(0, nrow(G12), ncol(G12))
+ G12.2[1:k,1:k] <- do.call(rbind, G12tmp)
+ G12 <- G12 + G12.2/n
+ }
+ Gamma <- rbind(cbind(G11, G12),
+ cbind(t(G12), G22))
+ addPar <- 0
+ } else {
+ lG <- cbind(lG, matrix(-tail(lam, nh), n, nh, byrow=TRUE))
+ G11 <- crossprod(rho2*lG, lG)/n
+ G12 <- t(G)/n + crossprod(rho2*lG, gt)/n
+ if (!is.null(family))
+ {
+ G12tmp <- lapply(1:k, function(i)
+ colSums(-rho1*tau2*Z[,i]*c(Z%*%lam[1:k])*Z))
+ G12.2 <- matrix(0, nrow(G12), ncol(G12))
+ G12.2[1:k,1:k] <- do.call(rbind, G12tmp)
+ G12 <- G12 + G12.2/n
+ }
+ if (momType == "balSample")
+ Xi <- rep(1,n)
+ else
+ Xi <- Z[,-1]
+ nj <- sum(Xi)
+ lam2 <- -sum(rho1)*tail(lam,nh)/nj
+ theta4 <- colSums(Xi*X)/nj
+ G13 <- rbind(matrix(0, 2*k-1, nh), -nj/n*diag(nh))
+ G23 <- matrix(0,q, nh)
+ G33 <- matrix(0, nh, nh)
+ Psi <- cbind(rho1*lG, rho1*gt,
+ Xi*sweep(X, 2, theta4, "-"))
+ Psi[,(2*k):(2*k+nh-1)] <- Psi[,(2*k):(2*k+nh-1)]-Xi%*%t(lam2)
+ Gamma <- rbind(cbind(G11, G12, G13),
+ cbind(t(G12), G22, G23),
+ cbind(t(G13), t(G23), G33))
+ addPar <- nh
+ }
+ list(Psi=Psi, Gamma=Gamma, k=length(theta), q=q, addPar=addPar, n=n)
+ }
+
+.vcovate <- function (object)
+ {
+ res <- .psiGam(object)
+ k <- res$k
+ q <- res$q
+ addPar <- res$addPar
+ qrPsi <- qr(res$Psi/sqrt(res$n))
+ piv <- sort.int(qrPsi$pivot, index.return=TRUE)$ix
+ R <- qr.R(qrPsi)[,piv]
+ T1 <- solve(res$Gamma, t(R))
+ V <- T1%*%t(T1)/res$n
+ allV <- list()
+ allV$vcov_par <- V[1:k, 1:k]
+ allV$vcov_lambda <- V[(k+addPar+1):(k+addPar+q), (k+addPar+1):(k+addPar+q)]
+ if (addPar > 0)
+ {
+ allV$vcov_Allpar <- V[1:(k+addPar), 1:(k+addPar)]
+ allV$vcov_Alllambda <- V[-(1:(k+addPar)), -(1:(k+addPar))]
+ }
+ allV
+ }
+
+
+vcov.ategel <- function(object, lambda = FALSE, robToMiss=TRUE, ...)
+ {
+ if (robToMiss)
+ {
+ return(vcov.gel(object, lambda))
+ } else {
+ object$lambda <- rep(0, length(object$lambda))
+ res <- .vcovate(object)
+ object$vcov_par <- res$vcov_par
+ object$vcov_lambda <- res$vcov_lambda
+ return(vcov.gel(object, lambda))
+ }
+ }
+
+summary.ategel <- function(object, robToMiss=TRUE, ...)
+ {
+ if (robToMiss)
+ {
+ ans <- summary.gel(object)
+ ans$typeDesc = paste(ans$typeDesc,
+ "\n(S.E. are robust to misspecification)", sep="")
+ } else {
+ object$vcov_par <- vcov(object, robToMiss=FALSE)
+ object$vcov_lambda <- vcov(object, TRUE, robToMiss=FALSE)
+ ans <- summary.gel(object)
+ ans$typeDesc = paste(ans$typeDesc,
+ "\n(S.E. are not robust to misspecification", sep="")
+ }
+ ans
+ }
+
+confint.ategel <- function (object, parm, level = 0.95, lambda = FALSE,
+ type = c("Wald", "invLR", "invLM", "invJ"), fact = 3,
+ corr = NULL, robToMiss=TRUE, ...)
+ {
+ type <- match.arg(type)
+ if (type=="Wald")
+ {
+ if (!robToMiss)
+ {
+ object$vcov_par <- vcov(object, robToMiss=FALSE)
+ object$vcov_lambda <- vcov(object, TRUE, robToMiss=FALSE)
+ }
+ return(confint.gel(object, parm, level, lambda, type, fact, corr, ...))
+ }
+ object$allArg$g <- .momentFctATE
+ object$allArg$y <- NULL
+ object$allArg$treat <- NULL
+ object$allArg$popMom <- NULL
+ object$allArg$momType <- NULL
+ object$allArg$x <- object$dat
+ return(confint.gel(object, parm, level, lambda, type, fact, corr, ...))
+ }
+
+
+marginal <- function(object, ...)
+ UseMethod("marginal")
+
+marginal.ategel <- function(object, ...)
+ {
+ family <- attr(object$dat, "family")
+ if (is.null(family))
+ return(summary(object)$coef)
+ k <- attr(object$dat, "k")
+ p0 <- family$linkinv(object$coef[1])
+ p1 <- family$linkinv(object$coef[1]+object$coef[2:k])
+ p01 <- family$mu.eta(object$coef[1])
+ p11 <- family$mu.eta(object$coef[1]+object$coef[2:k])
+ A <- cbind(p11-p01, p11)
+ V <- vcov(object, ...)[1:k,1:k]
+ sd0 <- p01*sqrt(V[1,1])
+ sdd <- sapply(1:(k-1), function(i)
+ sqrt(c(t(A[i,])%*%V[c(1,i+1), c(1, i+1)]%*%A[i,])))
+ coef <- cbind(c(p0,p1-p0), c(sd0, sdd))
+ coef <- cbind(coef, coef[,1]/coef[,2])
+ coef <- cbind(coef, 2*pnorm(-abs(coef[,3])))
+ colnames(coef) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
+ rownames(coef) <- c("Control",
+ paste("Treat", 1:(k-1) , " versus Control", sep=""))
+ coef
+ }
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/R/gel.R 2017-05-19 20:51:09 UTC (rev 101)
@@ -382,4 +382,30 @@
}
-
+.vcovGel <- function(gt, G, k1, k2, bw, pt=NULL,tol=1e-16)
+ {
+ q <- NCOL(gt)
+ n <- NROW(gt)
+ if (is.null(pt))
+ pt <- 1/n
+ G <- G/k1
+ gt <- gt*sqrt(pt*bw/k2)
+ qrGt <- qr(gt)
+ piv <- sort.int(qrGt$pivot, index.return=TRUE)$ix
+ R <- qr.R(qrGt)[,piv]
+ X <- forwardsolve(t(R), G)
+ Y <- forwardsolve(t(R), diag(q))
+ res <- lm.fit(X,Y)
+ u <- res$residuals
+ Sigma <- chol2inv(res$qr$qr)/n
+ diag(Sigma)[diag(Sigma)<0] <- tol
+ if (q==ncol(G))
+ {
+ SigmaLam <- matrix(0, q, q)
+ } else {
+ SigmaLam <- backsolve(R, u)/n*bw^2
+ diag(SigmaLam)[diag(SigmaLam)<0] <- tol
+ }
+ khat <- crossprod(R)
+ list(vcov_par=Sigma, vcov_lambda=SigmaLam,khat=khat)
+ }
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/R/getModel.R 2017-05-19 20:51:09 UTC (rev 101)
@@ -205,8 +205,10 @@
if (dat$ny > 1)
{
namey <- colnames(dat$x[,1:dat$ny, drop=FALSE])
- object$namesCoef <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
- object$namesgt <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny, dat$nh)), sep = "")
+ object$namesCoef <- paste(rep(namey, dat$k), "_",
+ rep(namex, rep(dat$ny, dat$k)), sep = "")
+ object$namesgt <- paste(rep(namey, dat$nh), "_",
+ rep(nameh, rep(dat$ny, dat$nh)), sep = "")
} else {
object$namesCoef <- namex
object$namesgt <- nameh
@@ -426,3 +428,109 @@
return(object)
}
+getModel.ateGel <- function(object, ...)
+ {
+ object$allArg <- c(object, list(...))
+ object$allArg$weights <- NULL
+ object$allArg$call <- NULL
+ if(is(object$g, "formula"))
+ {
+ dat <- getDat(object$g, object$x, data = object$data)
+ if (dat$ny>1 | dat$ny==0)
+ stop("You need one and only one dependent variable")
+ k <- dat$k
+ if (k>2 & object$momType=="ATT")
+ stop("Cannot compute ATT with multiple treatments")
+ if (attr(dat$mt, "intercept")!=1)
+ stop("An intercept is needed to compute the treatment effect")
+ if (!all(dat$x[,3:(k+1)] %in% c(0,1)))
+ stop("The treatment indicators can only take values 0 or 1")
+ if (colnames(dat$x)[k+2] == "(Intercept)")
+ {
+ dat$x <- dat$x[,-(k+2)]
+ dat$nh <- dat$nh-1
+ }
+ if (!is.null(object$popMom))
+ {
+ if (length(object$popMom)!=dat$nh)
+ stop("The dim. of the population moments does not match the dim. of the vector of covariates")
+ }
+ if (is.null(object$tet0))
+ {
+ tet0 <- lm(dat$x[,1]~dat$x[,3:(k+1)])$coef
+ tet0 <- c(tet0, colMeans(dat$x[,3:(k+1),drop=FALSE]))
+ names(tet0) <- NULL
+ object$tet0 <- tet0
+ } else {
+ if (length(object$tet0) != (2*k-1))
+ stop("tet0 should have a length equal to 2x(number of treatments)+1")
+ }
+ if (object$family != "linear")
+ {
+ if (any(!(dat$x[,1]%in%c(0,1))))
+ stop("For logit or probit, Y can only take 0s and 1s")
+ family <- binomial(link=object$family)
+ if (object$family == "logit")
+ family$mu.eta2 <- function(x, family) family$mu.eta(x)*(1-2*family$linkinv(x))
+ else
+ family$mu.eta2 <- function(x, family) -x*family$mu.eta(x)
+
+ } else {
+ family <- NULL
+ }
+ q <- dat$nh + 2*k+1
+ if (object$momType != "bal" | !is.null(object$popMom))
+ q <- q+dat$nh
+ object$gradv <- .DmomentFctATE
+ object$x <- dat
+ object$gradvf <- FALSE
+ object$gform<-object$g
+ namex <- colnames(dat$x[,2:(k+1)])
+ nameh <- colnames(dat$x[,(k+2):ncol(dat$x), drop=FALSE])
+ namesCoef <- c(namex, paste("TreatProb", 1:(k-1), sep=""))
+ namesgt <- paste(rep(paste("Treat", 1:(k-1), sep=""), rep(dat$nh, k-1)), "_", rep(nameh, k-1), sep="")
+ object$namesgt <- c(namesCoef,namesgt)
+ if (object$momType != "bal" | !is.null(object$popMom))
+ object$namesgt <- c(object$namesgt, paste("bal_", nameh, sep=""))
+ if (is.null(names(object$tet0)))
+ object$namesCoef <- namesCoef
+ else
+ object$namesCoef <- names(object$tet0)
+ attr(object$x,"ModelType") <- "linear"
+ attr(object$x, "k") <- k
+ attr(object$x, "q") <- q
+ attr(object$x, "n") <- nrow(dat$x)
+ attr(object$x, "momType") <- object$momType
+ attr(object$x, "popMom") <- object$popMom
+ attr(object$x, "family") <- family
+ } else {
+ stop("Not implemented yet for nonlinear regression")
+ }
+ if (object$momType == "ATT")
+ metname <- "ATT"
+ else
+ metname <- "ATE"
+ if (!is.null(object$popMom))
+ {
+ metname2 <- " with balancing based on population moments"
+ } else {
+ if (object$momType == "balSample")
+ metname2 <- " with balancing based on the sample moments"
+ else
+ metname2 <- " with unrestricted balancing"
+ }
+ metname3 <- paste("\nMethod: ", object$type, sep="")
+ if (!is.null(family))
+ metname3 <- paste(metname3, ", Family: Binomial with ", family$link, " link", sep="")
+ clname <- "baseGel.mod"
+ object$k1 <- 1
+ object$k2 <- 1
+ object$w <- kernel(1)
+ object$bwVal <- 1
+ object$g <- .momentFctATE
+ object$CGEL <- object$alpha
+ object$typeDesc <- paste(metname, metname2, metname3, sep="")
+ class(object) <- clname
+ return(object)
+ }
+
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/R/gmm.R 2017-05-19 20:51:09 UTC (rev 101)
@@ -46,10 +46,13 @@
z
}
-evalGmm <- function(g, x, t0, tetw=NULL, gradv=NULL, wmatrix = c("optimal","ident"), vcov=c("HAC","iid","TrueFixed"),
- kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen", "Tukey-Hanning"),crit=10e-7,bw = bwAndrews,
+evalGmm <- function(g, x, t0, tetw=NULL, gradv=NULL, wmatrix = c("optimal","ident"),
+ vcov=c("HAC","iid","TrueFixed"),
+ kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen",
+ "Tukey-Hanning"),crit=10e-7,bw = bwAndrews,
prewhite = FALSE, ar.method = "ols", approx="AR(1)",tol = 1e-7,
- model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, weightsMatrix = NULL, data)
+ model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE,
+ weightsMatrix = NULL, data)
{
TypeGmm = "baseGmm"
type <- "eval"
@@ -58,17 +61,19 @@
wmatrix <- match.arg(wmatrix)
if (is.null(tetw) & is.null(weightsMatrix))
stop("If the weighting matrix is not provided, you need to provide the vector of parameters tetw")
-
if(vcov=="TrueFixed" & is.null(weightsMatrix))
stop("TrueFixed vcov only for fixed weighting matrix")
if(!is.null(weightsMatrix))
wmatrix <- "optimal"
if(missing(data))
data<-NULL
- all_args<-list(data = data, g = g, x = x, t0 = t0, tetw = tetw, gradv = gradv, type = type, wmatrix = wmatrix, vcov = vcov, kernel = kernel,
- crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method, approx = approx,
- weightsMatrix = weightsMatrix, centeredVcov = centeredVcov, tol = tol, itermax = 100,
- optfct = NULL, model = model, X = X, Y = Y, call = match.call(), traceIter = NULL,
+ all_args<-list(data = data, g = g, x = x, t0 = t0, tetw = tetw, gradv = gradv,
+ type = type, wmatrix = wmatrix, vcov = vcov, kernel = kernel,
+ crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method,
+ approx = approx, weightsMatrix = weightsMatrix,
+ centeredVcov = centeredVcov, tol = tol, itermax = 100,
+ model = model, X = X, Y = Y, call = match.call(),
+ traceIter = NULL, optfct="optim",
eqConst = NULL, eqConstFullVcov = FALSE)
class(all_args)<-TypeGmm
Model_info<-getModel(all_args)
@@ -396,6 +401,7 @@
return(obj)
}
+
.momentFct <- function(tet, dat)
{
if (!is.null(attr(dat, "eqConst")))
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2017-04-28 14:44:25 UTC (rev 100)
+++ pkg/gmm/R/momentEstim.R 2017-05-19 20:51:09 UTC (rev 101)
@@ -779,28 +779,17 @@
G <- P$gradv(z$coefficients, P$x)
else
G <- P$gradv(z$coefficients, P$x, z$pt)
- khat <- crossprod(c(z$pt)*z$gt,z$gt)/(P$k2)*P$bwVa
- z$G <- G
- G <- G/P$k1
- kg <- solve(khat, G)
- z$vcov_par <- try(solve(crossprod(G, kg))/n, silent=TRUE)
- if (class(z$vcov_par) == "try-error")
+ allVcov <- try(.vcovGel(gt, G, P$k1, P$k2, P$bwVal, z$pt),
+ silent=TRUE)
+ if (class(allVcov) == "try-error")
{
- warning("Cannot compute the covariance matrix of the coefficients; matrix singular")
z$vcov_par <- matrix(NA, length(z$coefficients), length(z$coefficients))
- }
- if (dim(G)[1] == dim(G)[2])
- {
- z$vcov_lambda <- matrix(0, dim(G)[1], dim(G)[2])
+ z$vcov_lambda <- matrix(NA, length(z$lambda), length(z$lambda))
+ z$khat <- NULL
+ warning("Cannot compute the covariance matrices")
} else {
- z$vcov_lambda <- try(solve(khat, ( diag(ncol(khat)) - G %*% (z$vcov_par*n) %*% t(kg) ))/n*P$bwVal^2, silent=TRUE)
- if (class(z$vcov_lambda) == "try-error")
- {
- warning("Cannot compute the covariance matrix of the lambdas; matrix singular")
- z$vcov_lambda <- matrix(NA, length(z$lambda), length(z$lambda))
- }
- }
-
+ z <- c(z, allVcov)
+ }
z$weights <- P$w
z$bwVal <- P$bwVal
names(z$bwVal) <- "Bandwidth"
@@ -814,7 +803,6 @@
if(P$model) z$model <- P$x$mf
if(P$X) z$x <- as.matrix(P$x$x[,(P$x$ny+1):(P$x$ny+P$x$k)])
if(P$Y) z$y <- as.matrix(P$x$x[,1:P$x$ny])
- z$khat <- khat
class(z) <- paste(P$TypeGel, ".res", sep = "")
z$allArg <- P$allArg
return(z)
@@ -883,43 +871,31 @@
attr(x,"eqConst") <- NULL
z$specMod <- paste(z$specMod, "** Note: Covariance matrix computed for all coefficients based on restricted values \n Tests non-valid**\n\n")
}
-
if(P$gradvf)
G <- P$gradv(z$coefficients, x)
else
G <- P$gradv(z$coefficients, x, z$pt)
z$G <- G
- khat <- crossprod(c(z$pt)*z$gt, z$gt)/(P$k2)*P$bwVal
- G <- G/P$k1
- kg <- solve(khat, G)
- z$vcov_par <- try(solve(crossprod(G, kg))/n, silent=TRUE)
- if (class(z$vcov_par) == "try-error")
+ allVcov <- try(.vcovGel(gt, G, P$k1, P$k2, P$bwVal, z$pt),
+ silent=TRUE)
+ if (class(allVcov) == "try-error")
{
- warning("Cannot compute the covariance matrix of the coefficients; matrix singular")
z$vcov_par <- matrix(NA, length(z$coefficients), length(z$coefficients))
- }
- if (dim(G)[1] == dim(G)[2])
- {
- z$vcov_lambda <- matrix(0, dim(G)[1], dim(G)[2])
+ z$vcov_lambda <- matrix(NA, length(z$lambda), length(z$lambda))
+ z$khat <- NULL
+ warning("Cannot compute the covariance matrices")
} else {
- z$vcov_lambda <- try(solve(khat, ( diag(ncol(khat)) - G %*% (z$vcov_par*n) %*% t(kg) ))/n*P$bwVal^2, silent=TRUE)
- if (class(z$vcov_lambda) == "try-error")
- {
- warning("Cannot compute the covariance matrix of the lambdas; matrix singular")
- z$vcov_lambda <- matrix(NA, length(z$lambda), length(z$lambda))
- }
- }
+ z <- c(z, allVcov)
+ }
z$weights <- P$w
z$bwVal <- P$bwVal
names(z$bwVal) <- "Bandwidth"
-
dimnames(z$vcov_par) <- list(names(z$coefficients), names(z$coefficients))
dimnames(z$vcov_lambda) <- list(names(z$lambda), names(z$lambda))
if(P$X) z$x <- x
z$call <- P$call
z$k1 <- P$k1
z$k2 <- P$k2
- z$khat <- khat
z$CGEL <- P$CGEL
z$typeDesc <- P$typeDesc
z$allArg <- P$allArg
@@ -1108,17 +1084,9 @@
G <- P$gradv(z$coefficients, P$x)
else
G <- P$gradv(z$coefficients, P$x, z$pt)
- khat <- crossprod(c(z$pt)*z$gt,z$gt)/(P$k2)*P$bwVa
-
z$G <- G
- G <- G/P$k1
- kg <- solve(khat, G)
- z$vcov_par <- solve(crossprod(G, kg))/n
- if (dim(G)[1] == dim(G)[2])
- z$vcov_lambda <- matrix(0, dim(G)[1], dim(G)[2])
- else
- z$vcov_lambda <- solve(khat, ( diag(ncol(khat)) - G %*% (z$vcov_par*n) %*% t(kg) ))/n*P$bwVal^2
-
+ allVcov <- .vcovGel(gt, G, P$k1, P$k2, P$bwVal, z$pt)
+ z <- c(z, allVcov)
z$weights <- P$w
z$bwVal <- P$bwVal
names(z$bwVal) <- "Bandwidth"
@@ -1134,7 +1102,6 @@
if(P$X) z$x <- as.matrix(P$x$x[,(P$x$ny+1):(P$x$ny+P$x$k)])
if(P$Y) z$y <- as.matrix(P$x$x[,1:P$x$ny])
}
- z$khat <- khat
return(z)
}
Added: pkg/gmm/data/nsw.rda
===================================================================
(Binary files differ)
Property changes on: pkg/gmm/data/nsw.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/gmm/man/ATEgel.Rd
===================================================================
--- pkg/gmm/man/ATEgel.Rd (rev 0)
+++ pkg/gmm/man/ATEgel.Rd 2017-05-19 20:51:09 UTC (rev 101)
@@ -0,0 +1,135 @@
+\name{ATEgel}
+
+\alias{ATEgel}
+
+\title{ATE with Generalized Empirical Likelihood estimation}
+
+\description{
+Function to estimate the average treatment effect with the sample being
+balanced by GEL.
+}
+\usage{
+ATEgel(g, balm, y=NULL, treat=NULL, tet0=NULL,momType=c("bal","balSample","ATT"),
+ popMom = NULL, family=c("linear","logit", "probit"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD"),
+ tol_lam = 1e-9, tol_obj = 1e-9, tol_mom = 1e-9, maxiterlam = 100,
+ optfct = c("optim", "nlminb"),
+ optlam = c("nlminb", "optim", "iter"), data=NULL,
+ Lambdacontrol = list(),
+ model = TRUE, X = FALSE, Y = FALSE, ...)
+}
+\arguments{
+\item{g}{A formula as \code{y~z}, where code{y} is the response and
+ \code{z} the treatment indicator. If there is more than one
+ treatment, more indicators can be added or \code{z} can be set as a
+ factor. It can also be of the form
+ \code{g(theta, y, z)} for non-linear models. It is however, not implemented yet.}
+
+\item{balm}{A formula for the moments to be balanced between the treated
+ and control groups (see details)}
+
+\item{y}{The response variable when \code{g} is a function. Not
+ implemented yet}
+
+\item{treat}{The treatment indicator when \code{g} is a function. Not
+ implemented yet}
+
+\item{tet0}{A \eqn{3 \times 1} vector of starting values. If not
+ provided, they are obtained using an OLS regression}
+
+\item{momType}{How the moments of the covariates should be balanced. By
+ default, it is simply balanced without restriction. Alternatively,
+ moments can be set equal to the sample moments of the whole sample, or
+ to the sample moments of the treated group. The later will produce
+ the average treatment effect of the treated (ATT)}
+
+\item{popMom}{A vector of population moments to use for balancing. It
+ can be used of those moments are available from a census, for
+ example. When available, it greatly improves efficiency.}
+
+\item{family}{By default, the outcome is linearly related to the
+ treatment indicators. If the outcome is binary, it is possible to use
+ the estimating equations of either the logit or probit model.}
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 101
More information about the Gmm-commits
mailing list