[Gmm-commits] r138 - in pkg/gmm: . R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Apr 24 16:40:35 CEST 2019


Author: chaussep
Date: 2019-04-24 16:40:34 +0200 (Wed, 24 Apr 2019)
New Revision: 138

Added:
   pkg/gmm/man/getImpProb.Rd
Modified:
   pkg/gmm/NAMESPACE
   pkg/gmm/NEWS
   pkg/gmm/R/Methods.gel.R
   pkg/gmm/R/ategel.R
   pkg/gmm/R/gel.R
   pkg/gmm/man/ATEgel.Rd
   pkg/gmm/man/gel.Rd
   pkg/gmm/man/getLamb.Rd
Log:
Changed the algorithms for GEL-CUE and add getImpProb method; see NEWS

Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/NAMESPACE	2019-04-24 14:40:34 UTC (rev 138)
@@ -5,18 +5,28 @@
 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, 
-       residuals.gmm, fitted.gel, residuals.gel, plot.gmm, plot.gel,formula.gmm, formula.gel, charStable, specTest, 
-       specTest.gmm, specTest.gel, print.specTest, momentEstim.baseGmm.twoStep, momentEstim.baseGmm.twoStep.formula,
-       momentEstim.baseGmm.iterative.formula, momentEstim.baseGmm.iterative, momentEstim.baseGmm.cue.formula, 
-       momentEstim.baseGmm.cue, getModel.baseGmm, getModel.baseGel, getModel.constGmm, getModel.constGel,
-       FinRes.baseGmm.res, momentEstim.baseGel.mod, momentEstim.baseGel.modFormula,tsls,summary.tsls, print.summary.tsls,
-       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, getModel.ateGel, ATEgel,
-       summary.ategel, vcov.ategel, confint.ategel, marginal, marginal.ategel,checkConv)
+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, residuals.gmm,
+fitted.gel, residuals.gel, plot.gmm, plot.gel,formula.gmm,
+formula.gel, charStable, specTest, specTest.gmm, specTest.gel,
+print.specTest, momentEstim.baseGmm.twoStep,
+momentEstim.baseGmm.twoStep.formula,
+momentEstim.baseGmm.iterative.formula, momentEstim.baseGmm.iterative,
+momentEstim.baseGmm.cue.formula, momentEstim.baseGmm.cue,
+getModel.baseGmm, getModel.baseGel, getModel.constGmm,
+getModel.constGel, FinRes.baseGmm.res, momentEstim.baseGel.mod,
+momentEstim.baseGel.modFormula,tsls,summary.tsls, print.summary.tsls,
+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, getModel.ateGel, ATEgel,
+summary.ategel, vcov.ategel, confint.ategel, marginal,
+marginal.ategel,checkConv,getImpProb, getImpProb.gel)
 
 S3method(marginal, ategel)
 S3method(summary, gmm)
@@ -51,6 +61,7 @@
 S3method(formula, gel)
 S3method(specTest, gmm)
 S3method(specTest, gel)
+S3method(getImpProb, gel)
 S3method(print, specTest)
 S3method(print, confint)
 S3method(FinRes, baseGmm.res)

Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/NEWS	2019-04-24 14:40:34 UTC (rev 138)
@@ -8,6 +8,8 @@
   bwNeweyWest. The new version, takes it into account. It is not possible to obtain the same sandwich covariance matrix with gmm and OLS.
 o In previous versions, the weighting matrix and, as a result, the standard errors were based on the small sample adjusted vcovHAC for weakly dependent
   processes. The adjustment is not justified for GMM. The option adjust=FALSE is therefore added to vcovHAC in .myKernHAC()
+o Added a method getImpProb, which extract the implied probabilities. It allows  for negative probabilities as it often happens with CUE
+o lambda for CUE is solved analytically (it is minus the projection of gt on a vector of ones). The option RCue for optlam can be used to restrict the implied probabilities to be non-negative.
 
 Changes in version 1.6-1
 

Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/R/Methods.gel.R	2019-04-24 14:40:34 UTC (rev 138)
@@ -298,5 +298,32 @@
   }
 
 
+getImpProb <- function(object, ...)
+    UseMethod("getImpProb")
 
+getImpProb.gel <- function(object, posProb=TRUE, normalize=TRUE,
+                           checkConv=FALSE, ...)
+    {
+        if (!normalize || (object$type == "CUE" && !posProb))
+            {
+                n <- NROW(object$gt)
+                pt <- -.rho(object$gt, object$lambda, type = object$type,
+                            derive = 1, k = object$k1/object$k2)/n
+                if (object$type == "CUE" && posProb)
+                    {
+                        eps <- -length(pt)*min(min(pt),0)
+                        pt <- (pt+eps/length(pt))/(1+eps)
+                    }
+                if (normalize)
+                    pt <- pt/sum(pt)
+            } else {
+                pt <- object$pt
+            }
+        if (checkConv)
+            attr(pt, "convergence") <- list(pt=sum(pt),
+                                            ptgt=colSums(pt*as.matrix(object$gt)))
+        pt
+    }
+
+
         

Modified: pkg/gmm/R/ategel.R
===================================================================
--- pkg/gmm/R/ategel.R	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/R/ategel.R	2019-04-24 14:40:34 UTC (rev 138)
@@ -4,7 +4,7 @@
                    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", "Wu"),
+                   optlam = c("nlminb", "optim", "iter", "Wu", "RCue"),
                    data=NULL, Lambdacontrol = list(),
                    model = TRUE, X = FALSE, Y = FALSE, ...)
     {
@@ -325,7 +325,7 @@
         coef
     }
 
-checkConv <- function(obj, tolConv=1e-4,verbose=TRUE)
+checkConv <- function(obj, tolConv=1e-4,verbose=TRUE, ...)
     {
         if (!any(class(obj)=="ategel"))
             stop("The function is for ategel objects produced by ATEgel()")
@@ -337,7 +337,7 @@
         nZ <- attr(obj$dat, "k")-1
         z <- dat[,3:(2+nZ),drop=FALSE]
         x <- dat[,-(1:(2+nZ)),drop=FALSE]
-        pt <- obj$pt
+        pt <- getImpProb(obj, ...)
         pt1 <- lapply(1:nZ, function(i) pt[z[,i]==1]/sum(pt[z[,i]==1]))
         pt0 <- pt[rowSums(z)==0]/sum(pt[rowSums(z)==0])
         m0 <- colSums(x[rowSums(z)==0,,drop=FALSE]*pt0)

Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/R/gel.R	2019-04-24 14:40:34 UTC (rev 138)
@@ -149,18 +149,56 @@
                         mean(.rho(gt,-M,derive=0,type="EL",k=K))))        
     }
 
+.CUE_lam <- function(gt, K=1, control=list(), numSol=FALSE)
+    {
+        if (!numSol)
+            {
+                q <- qr(gt)
+                n <- nrow(gt)
+                l0 <- -qr.coef(q, rep(1,n))
+                conv <- list(convergence=0)
+            } else {
+		fct <- function(l,X,k)
+                    {
+			r0 <- .rho(X,l,derive=0,type="CUE",k=k)
+			-mean(r0)
+                    }
+		Dfct <- function(l,X,k)
+                    {
+			r1 <- .rho(X,l,derive=1,type="CUE",k=k)
+		        -colMeans(r1*X)
+                    }
+                ui <- gt
+                ci <- rep(-1,nrow(gt))
+                res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,ui,ci,
+                                   control=control,
+                                   X=gt,k=K)
+                l0 <- res$par
+                conv <- list(convergence = res$convergence, counts = res$counts,
+                             message = res$message)                
+            }        
+        list(lambda = l0, convergence = conv, obj =
+                 mean(.rho(gt,l0,derive=0,type="CUE",k=K)))
+    }
+
 getLamb <- function(gt, l0, type = c('EL', 'ET', 'CUE', "ETEL", "HD", "ETHD"),
                     tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7, 
-                    k = 1, method = c("nlminb", "optim", "iter", "Wu"), control=list())
+                    k = 1, method = c("nlminb", "optim", "iter", "Wu", "RCue"),
+                    control=list())
     {
 	method <- match.arg(method)
-	if(is.null(dim(gt)))
-            gt <- matrix(gt,ncol=1)
-
+        gt <- as.matrix(gt)
         if (method == "Wu" & type != "EL")
             stop("Wu (2005) method to compute Lambda is for EL only")
+        if (method == "RCue" & type != "CUE")
+            {
+                warning("RCue is for CUE only. optlam set to optim")
+                method <- "optim"
+            }
         if (method == "Wu")
             return(.Wu(gt, tol_lam, maxiterlam, k))
+        if (type == "CUE")
+            return(.CUE_lam(gt, k, control, method=="RCue"))
 	if (method == "iter")
             {
 		if ((type == "ETEL") | (type == "ETHD"))
@@ -279,7 +317,8 @@
                 tol_weights = 1e-7, tol_lam = 1e-9, tol_obj = 1e-9, 
 		tol_mom = 1e-9, maxiterlam = 100, constraint = FALSE,
                 optfct = c("optim", "optimize", "nlminb"), 
-                optlam = c("nlminb", "optim", "iter", "Wu"), data, Lambdacontrol = list(),
+                optlam = c("nlminb", "optim", "iter", "Wu", "RCue"), data,
+                Lambdacontrol = list(),
                 model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", alpha = NULL,
                 eqConst = NULL, eqConstFullVcov = FALSE, onlyCoefficients=FALSE, ...)
     {
@@ -303,7 +342,8 @@
                          optlam = optlam, model = model, X = X, Y = Y,
                          TypeGel = TypeGel, call = match.call(),
                          Lambdacontrol = Lambdacontrol, alpha = alpha, data = data,
-                         eqConst = eqConst, eqConstFullVcov = eqConstFullVcov, onlyCoefficients=onlyCoefficients)
+                         eqConst = eqConst, eqConstFullVcov = eqConstFullVcov,
+                         onlyCoefficients=onlyCoefficients)
 	class(all_args)<-TypeGel
 	Model_info<-getModel(all_args)
 	z <- momentEstim(Model_info, ...)
@@ -318,7 +358,7 @@
                     approx = c("AR(1)", "ARMA(1,1)"), prewhite = 1,
                     ar.method = "ols", tol_weights = 1e-7, tol_lam = 1e-9,
                     tol_obj = 1e-9, tol_mom = 1e-9, maxiterlam = 100,
-                    optlam = c("nlminb", "optim", "iter", "Wu"), data,
+                    optlam = c("nlminb", "optim", "iter", "Wu", "RCue"), data,
                     Lambdacontrol = list(),
                     model = TRUE, X = FALSE, Y = FALSE, alpha = NULL, ...)
     {
@@ -423,7 +463,6 @@
         pt
     }
 
-                
 .vcovGel <- function(gt, G, k1, k2, bw, pt=NULL,tol=1e-16)
     {
         q <- NCOL(gt)

Modified: pkg/gmm/man/ATEgel.Rd
===================================================================
--- pkg/gmm/man/ATEgel.Rd	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/man/ATEgel.Rd	2019-04-24 14:40:34 UTC (rev 138)
@@ -15,10 +15,10 @@
                    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", "Wu"), data=NULL,
+                   optlam = c("nlminb", "optim", "iter", "Wu", "RCue"), data=NULL,
                    Lambdacontrol = list(),
                    model = TRUE, X = FALSE, Y = FALSE, ...)
-checkConv(obj, tolConv=1e-4, verbose=TRUE)
+checkConv(obj, tolConv=1e-4, verbose=TRUE, ...)
 }
 \arguments{
 \item{g}{A formula as \code{y~z}, where code{y} is the response and
@@ -77,7 +77,12 @@
 \item{tol_mom}{It is the tolerance for the moment condition \eqn{\sum_{t=1}^n p_t g(\theta(x_t)=0}, where \eqn{p_t=\frac{1}{n}D\rho(<g_t,\lambda>)} is the implied probability. It adds a penalty if the solution diverges from its goal.}
 
 \item{optlam}{Algorithm used to solve for the lagrange multiplier in
-  \code{\link{getLamb}}. The algorithm Wu is only for \code{type="EL"}. }
+  \code{\link{getLamb}}. The algorithm Wu is only for
+  \code{type="EL"}. By default, the analytical solution is used for
+  lambda, which may result in negative implied probabilities. The "RCue"
+  method uses \code{\link{constrOptim}} as for "EL" to restrict implied
+  probabiltities to be positive. If the non-negativity constraint is
+  binding, the samples will not be properly balanced.}
 
 \item{data}{A data.frame or a matrix with column names (Optional). }
 
@@ -89,7 +94,9 @@
 
 \item{tolConv}{The tolerance for comparing moments between groups}
 
-\item{...}{More options to give to \code{\link{optim}} or \code{\link{nlminb}}.}
+\item{...}{More options to give to \code{\link{optim}} or
+  \code{\link{nlminb}}. In \code{checkConv}, they are options passed to
+  \code{\link{getImpProb}}.}
 
 }
 

Modified: pkg/gmm/man/gel.Rd
===================================================================
--- pkg/gmm/man/gel.Rd	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/man/gel.Rd	2019-04-24 14:40:34 UTC (rev 138)
@@ -15,7 +15,7 @@
     approx = c("AR(1)", "ARMA(1,1)"), prewhite = 1, ar.method = "ols", 
     tol_weights = 1e-7, tol_lam = 1e-9, tol_obj = 1e-9, tol_mom = 1e-9, 
     maxiterlam = 100, constraint = FALSE, optfct = c("optim", "optimize", 
-    "nlminb"), optlam = c("nlminb", "optim", "iter", "Wu"), data,
+    "nlminb"), optlam = c("nlminb", "optim", "iter", "Wu", "RCue"), data,
     Lambdacontrol = list(), model = TRUE, X = FALSE, Y = FALSE,
     TypeGel = "baseGel", alpha = NULL, eqConst = NULL,
     eqConstFullVcov = FALSE, onlyCoefficients=FALSE, ...)
@@ -25,8 +25,8 @@
         approx = c("AR(1)", "ARMA(1,1)"), prewhite = 1,
         ar.method = "ols", tol_weights = 1e-7, tol_lam = 1e-9, tol_obj = 1e-9, 
         tol_mom = 1e-9, maxiterlam = 100, optlam = c("nlminb", "optim",
-        "iter", "Wu"), data, Lambdacontrol = list(), model = TRUE, X = FALSE,
-        Y = FALSE, alpha = NULL, ...)
+        "iter", "Wu", "RCue"), data, Lambdacontrol = list(),
+        model = TRUE, X = FALSE, Y = FALSE, alpha = NULL, ...)
 }
 \arguments{
 \item{g}{A function of the form \eqn{g(\theta,x)} and which returns a \eqn{n \times q} matrix with typical element \eqn{g_i(\theta,x_t)} for \eqn{i=1,...q} and \eqn{t=1,...,n}. This matrix is then used to build the q sample moment conditions. It can also be a formula if the model is linear (see details below).  }
@@ -73,7 +73,12 @@
 \item{tol_mom}{It is the tolerance for the moment condition \eqn{\sum_{t=1}^n p_t g(\theta(x_t)=0}, where \eqn{p_t=\frac{1}{n}D\rho(<g_t,\lambda>)} is the implied probability. It adds a penalty if the solution diverges from its goal.}
 
 \item{optlam}{Algorithm used to solve for the lagrange multiplier in
-  \code{\link{getLamb}}. The algorithm Wu is only for \code{type="EL"}. }
+  \code{\link{getLamb}}. The algorithm Wu is only for
+  \code{type="EL"}. The last method, is for "CUE". By default, the
+  analytical solution is used for lambda, which may result in negative
+  implied probabilities. The "RCue" method uses
+  \code{\link{constrOptim}} as for "EL" to restrict implied
+  probabiltities to be positive.}
 
 \item{data}{A data.frame or a matrix with column names (Optional). }
 

Added: pkg/gmm/man/getImpProb.Rd
===================================================================
--- pkg/gmm/man/getImpProb.Rd	                        (rev 0)
+++ pkg/gmm/man/getImpProb.Rd	2019-04-24 14:40:34 UTC (rev 138)
@@ -0,0 +1,57 @@
+\name{getImpProb}
+\alias{getImpProb}
+\alias{getImpProb.gel}
+
+\title{Implied Probabilities} \description{ It computes the implied
+probabilities from objects of class \code{gel} with additional options.
+} 
+
+\usage{ 
+\method{getImpProb}{gel}(object, posProb=TRUE, normalize=TRUE,
+                         checkConv=FALSE,...)
+}
+
+
+
+\arguments{ \item{object}{Object of class
+\code{gel}.}  \item{posProb}{Should the implied probabilities be
+transformed into positive probabilities?}  \item{normalize}{Should we
+normalize the probabilities so that they sum to one?}
+\item{checkConv}{Should we add the attribute convergence to check the
+  sum of the probabilities and the weighted sum of the moment conditions?}
+\item{...}{Additional arguments to pass to other methods}
+}
+
+\value{
+  A vector af implied probabilities.
+}
+
+\references{
+  Newey, W.K. and Smith, R.J. (2004), Higher Order Properties of GMM and
+  Generalized Empirical Likelihood Estimators. \emph{Econometrica},
+  \bold{72}, 219-255. 
+}
+
+
+
+\examples{
+
+#################
+n = 500
+phi<-c(.2,.7)
+thet <- 0
+sd <- .2
+x <- matrix(arima.sim(n=n,list(order=c(2,0,1),ar=phi,ma=thet,sd=sd)),ncol=1)
+y <- x[7:n]
+ym1 <- x[6:(n-1)]
+ym2 <- x[5:(n-2)]
+
+H <- cbind(x[4:(n-3)], x[3:(n-4)], x[2:(n-5)], x[1:(n-6)])
+g <- y ~ ym1 + ym2
+x <- H
+t0 <- c(0,.5,.5)
+
+res <- gel(g, x, t0)
+pt <- getImpProb(res)
+}
+

Modified: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd	2018-09-28 19:51:33 UTC (rev 137)
+++ pkg/gmm/man/getLamb.Rd	2019-04-24 14:40:34 UTC (rev 138)
@@ -10,7 +10,7 @@
 \usage{
 getLamb(gt, l0, type = c("EL","ET","CUE", "ETEL", "HD","ETHD"),
         tol_lam = 1e-7, maxiterlam = 100, 
-	tol_obj = 1e-7, k = 1, method = c("nlminb", "optim", "iter", "Wu"),
+	tol_obj = 1e-7, k = 1, method = c("nlminb", "optim", "iter", "Wu","RCue"),
         control = list())
 }
 \arguments{
@@ -36,7 +36,11 @@
   \code{\link{constrOptim}} is called to prevent \eqn{log(1-gt'\lambda)}
   from producing NA. The gradient and hessian is provided to
   \code{nlminb} which speed up the convergence. The latter is therefore
-  the default value. "Wu" is for "EL" only. It uses the algorithm of Wu (2005).}
+  the default value. "Wu" is for "EL" only. It uses the algorithm of Wu
+  (2005). The last method, is for "CUE". By default, the analytical
+  solution is used for lambda, which may result in negative implied
+  probabilities. The "RCue" method uses \code{\link{constrOptim}} as for
+  "EL" to restrict implied probabiltities to be positive.}
 
 \item{control}{Controls to send to \code{\link{optim}}, \code{\link{nlminb}} or \code{\link{constrOptim}}}
 }



More information about the Gmm-commits mailing list