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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Jun 9 15:26:24 CEST 2017


Author: chaussep
Date: 2017-06-09 15:26:23 +0200 (Fri, 09 Jun 2017)
New Revision: 103

Modified:
   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:
Added the Wu (2005) method for EL

Modified: pkg/gmm/R/ategel.R
===================================================================
--- pkg/gmm/R/ategel.R	2017-05-26 19:49:59 UTC (rev 102)
+++ pkg/gmm/R/ategel.R	2017-06-09 13:26:23 UTC (rev 103)
@@ -3,7 +3,8 @@
                    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(),
+                   optlam = c("nlminb", "optim", "iter", "Wu"),
+                   data=NULL, Lambdacontrol = list(),
                    model = TRUE, X = FALSE, Y = FALSE, ...)
     {
         type <- match.arg(type)
@@ -28,16 +29,22 @@
 	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
+        res <- try(.vcovate(z), silent=TRUE)
+        if (any(class(res)=="try-error"))
+            {
+                warning("Could not compute the robust-to misspecification standard errors")
+            } else {
+                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")
+        #k <- attr(dat, "k")
+        k <- dat$k
         ZT <- c(x[,2:(k+1),drop=FALSE]%*%tet[1:k])
         if (is.null(attr(dat, "family")))
             e <- x[,1] - ZT
@@ -74,7 +81,8 @@
         if (is.null(pt))
             pt <- rep(1/nrow(dat$x), nrow(dat$x))
         x <- dat$x
-        k <- attr(dat, "k")
+        #k <- attr(dat, "k")
+        k <- dat$k
         q <- dat$nh*(k-1)+2*k-1
         Z <- x[,2:(k+1), drop=FALSE]
         ZT <- c(Z%*%tet[1:k])
@@ -99,7 +107,8 @@
 .DmomentFctATE2 <- function(tet, dat, pt=NULL)
     {
         G <- .DmomentFctATE(tet, dat, pt)
-        k <- attr(dat, "k")
+        #k <- attr(dat, "k")
+        k <- dat$k
         q <- dat$nh*(k-1)+2*k-1
         if (is.null(pt))
             pt <- rep(1/nrow(dat$x), nrow(dat$x))
@@ -261,6 +270,7 @@
         object$allArg$treat <- NULL
         object$allArg$popMom <- NULL
         object$allArg$momType <- NULL
+        object$allArg$family <- NULL
         object$allArg$x <- object$dat        
         return(confint.gel(object, parm, level, lambda, type, fact, corr, ...))
     }

Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R	2017-05-26 19:49:59 UTC (rev 102)
+++ pkg/gmm/R/gel.R	2017-06-09 13:26:23 UTC (rev 103)
@@ -114,14 +114,53 @@
                     obj = mean(.rho(gt,res$par, derive=0,type=type,k=k))))	
     }
 
+.Wu <- function(gt, tol_lam = 1e-8, maxiterlam = 50, K=1)
+    {
+        u <- as.matrix(gt)
+        n=length(nrow(u))
+        M=rep(0,ncol(u))
+        dif=1
+        tol=tol_lam
+        k=0
+        while(dif>tol & k<=maxiterlam)
+            {
+                D1=t(u)%*%((1/(1+u%*%M))*rep(1,n))
+                DD=-t(u)%*%(c((1/(1+u%*%M)^2))*u)
+                D2=solve(DD,D1,tol=1e-40)
+                dif=max(abs(D2))
+                rule=1
+                while(rule>0)
+                    {
+                        rule=0
+                        if(min(1+t(M-D2)%*%t(u))<=0) rule=rule+1
+                        if(rule>0) D2=D2/2
+                    }
+                M=M-D2
+                k=k+1
+            }
+        if(k>=maxiterlam)
+            {
+                M=rep(0,ncol(u))
+                conv = list(convergence=1)
+            } else {
+                conv = list(convergence=0)
+            }
+	return(list(lambda = c(-M), convergence = conv, obj =
+                        mean(.rho(gt,-M,derive=0,type="EL",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"), control=list())
+                    k = 1, method = c("nlminb", "optim", "iter", "Wu"), control=list())
     {
 	method <- match.arg(method)
 	if(is.null(dim(gt)))
             gt <- matrix(gt,ncol=1)
-        
+
+        if (method == "Wu" & type != "EL")
+            stop("Wu (2005) method to compute Lambda is for EL only")
+        if (method == "Wu")
+            return(.Wu(gt, tol_lam, maxiterlam, k))
 	if (method == "iter")
             {
 		if ((type == "ETEL") | (type == "ETHD"))
@@ -192,9 +231,11 @@
                         mean(.rho(gt,l0,derive=0,type=type,k=k))))
     }
 
-smoothG <- function (x, bw = bwAndrews, prewhite = 1, ar.method = "ols", weights = weightsAndrews,
-			kernel = c("Bartlett", "Parzen", "Truncated", "Tukey-Hanning"), approx = c("AR(1)", "ARMA(1,1)"),
-			tol = 1e-7) 
+smoothG <- function (x, bw = bwAndrews, prewhite = 1, ar.method = "ols",
+                     weights = weightsAndrews,
+                     kernel = c("Bartlett", "Parzen", "Truncated", "Tukey-Hanning"),
+                     approx = c("AR(1)", "ARMA(1,1)"),
+                     tol = 1e-7) 
     {
 	kernel <- match.arg(kernel)
 	approx <- match.arg(approx)
@@ -238,7 +279,7 @@
                 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"), data, Lambdacontrol = list(),
+                optlam = c("nlminb", "optim", "iter", "Wu"), data, Lambdacontrol = list(),
                 model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", alpha = NULL,
                 eqConst = NULL, eqConstFullVcov = FALSE, ...)
     {
@@ -277,7 +318,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"), data, Lambdacontrol = list(),
+                    optlam = c("nlminb", "optim", "iter", "Wu"), data,
+                    Lambdacontrol = list(),
                     model = TRUE, X = FALSE, Y = FALSE, alpha = NULL, ...)
     {
 	type <- match.arg(type)

Modified: pkg/gmm/man/ATEgel.Rd
===================================================================
--- pkg/gmm/man/ATEgel.Rd	2017-05-26 19:49:59 UTC (rev 102)
+++ pkg/gmm/man/ATEgel.Rd	2017-06-09 13:26:23 UTC (rev 103)
@@ -14,7 +14,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"), data=NULL,
+                   optlam = c("nlminb", "optim", "iter", "Wu"), data=NULL,
                    Lambdacontrol = list(),
                    model = TRUE, X = FALSE, Y = FALSE, ...)
 }
@@ -67,7 +67,8 @@
 
 \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}{The default is "iter" which solves for \eqn{\lambda} using the Newton iterative method \code{\link{getLamb}}. If set to "numeric", the algorithm \code{\link{optim}} is used to compute \eqn{\lambda} instead.}
+\item{optlam}{Algorithm used to solve for the lagrange multiplier in
+  \code{\link{getLamb}}. The algorithm Wu is only for \code{type="EL"}. }
 
 \item{data}{A data.frame or a matrix with column names (Optional). }
 
@@ -109,6 +110,10 @@
 Schennach, Susanne, M. (2007), Point Estimation with Exponentially Tilted Empirical Likelihood.
 \emph{Econometrica}, \bold{35}, 634-672.
 
+Wu, C. (2005), Algorithms and R codes for the pseudo empirical
+likelihood method in survey sampling.
+\emph{Survey Methodology}, \bold{31}(2), page 239.
+
 Chausse, P. (2010), Computing Generalized Method of Moments and Generalized Empirical Likelihood with R.
  \emph{Journal of Statistical Software}, \bold{34}(11), 1--35.
  URL \url{http://www.jstatsoft.org/v34/i11/}.

Modified: pkg/gmm/man/gel.Rd
===================================================================
--- pkg/gmm/man/gel.Rd	2017-05-26 19:49:59 UTC (rev 102)
+++ pkg/gmm/man/gel.Rd	2017-06-09 13:26:23 UTC (rev 103)
@@ -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"), data,
+    "nlminb"), optlam = c("nlminb", "optim", "iter", "Wu"), data,
     Lambdacontrol = list(), model = TRUE, X = FALSE, Y = FALSE,
     TypeGel = "baseGel", alpha = NULL, eqConst = NULL,
     eqConstFullVcov = FALSE, ...)
@@ -25,7 +25,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"), data, Lambdacontrol = list(), model = TRUE, X = FALSE,
+        "iter", "Wu"), data, Lambdacontrol = list(), model = TRUE, X = FALSE,
         Y = FALSE, alpha = NULL, ...)
 }
 \arguments{
@@ -72,7 +72,8 @@
 
 \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}{The default is "iter" which solves for \eqn{\lambda} using the Newton iterative method \code{\link{getLamb}}. If set to "numeric", the algorithm \code{\link{optim}} is used to compute \eqn{\lambda} instead.}
+\item{optlam}{Algorithm used to solve for the lagrange multiplier in
+  \code{\link{getLamb}}. The algorithm Wu is only for \code{type="EL"}. }
 
 \item{data}{A data.frame or a matrix with column names (Optional). }
 
@@ -181,6 +182,10 @@
 Schennach, Susanne, M. (2007), Point Estimation with Exponentially Tilted Empirical Likelihood.
 \emph{Econometrica}, \bold{35}, 634-672.
 
+Wu, C. (2005), Algorithms and R codes for the pseudo empirical
+likelihood method in survey sampling.
+\emph{Survey Methodology}, \bold{31}(2), page 239.
+
 Zeileis A (2006), Object-oriented Computation of Sandwich Estimators.
 \emph{Journal of Statistical Software}, \bold{16}(9), 1--16.
 URL \url{http://www.jstatsoft.org/v16/i09/}.

Modified: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd	2017-05-26 19:49:59 UTC (rev 102)
+++ pkg/gmm/man/getLamb.Rd	2017-06-09 13:26:23 UTC (rev 103)
@@ -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"),
+	tol_obj = 1e-7, k = 1, method = c("nlminb", "optim", "iter", "Wu"),
         control = list())
 }
 \arguments{
@@ -30,7 +30,13 @@
 
 \item{k}{It represents the ratio k1/k2, where \eqn{k1=\int_{-\infty}^{\infty} k(s)ds} and \eqn{k2=\int_{-\infty}^{\infty} k(s)^2 ds}. See Smith(2004).}
 
-\item{method}{The iterative procedure uses a Newton method for solving the FOC. It i however recommended to use \code{optim} or \code{nlminb}. If type is set to "EL" and method to "optim", \code{\link{constrOptim}} is called to prevent \eqn{log(1-gt'\lambda)} from producing NA. The gradient and hessian is providewd to  \code{nlminb} which speed up the convergence. The latter is therefore the default value.}
+\item{method}{The iterative procedure uses a Newton method for solving
+  the FOC. It i however recommended to use \code{optim} or
+  \code{nlminb}. If type is set to "EL" and method to "optim",
+  \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).}
 
 \item{control}{Controls to send to \code{\link{optim}}, \code{\link{nlminb}} or \code{\link{constrOptim}}}
 }
@@ -57,6 +63,10 @@
 
 Smith, R.J. (2004), GEL Criteria for Moment Condition Models. \emph{Working paper, CEMMAP}.
 
+Wu, C. (2005), Algorithms and R codes for the pseudo empirical
+likelihood method in survey sampling.
+\emph{Survey Methodology}, \bold{31}(2), page 239.
+
 }
 
 \examples{



More information about the Gmm-commits mailing list