[Gmm-commits] r139 - in pkg/gmm: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 7 21:43:15 CEST 2019
Author: chaussep
Date: 2019-05-07 21:43:14 +0200 (Tue, 07 May 2019)
New Revision: 139
Added:
pkg/gmm/src/
pkg/gmm/src/Makevars
pkg/gmm/src/lambda_met.f
Modified:
pkg/gmm/NAMESPACE
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 an algorithm to restrict the implied probabilities for CUE to be non negative, converted Wu method in Fortran
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2019-04-24 14:40:34 UTC (rev 138)
+++ pkg/gmm/NAMESPACE 2019-05-07 19:43:14 UTC (rev 139)
@@ -1,3 +1,4 @@
+useDynLib(gmm)
import(stats)
importFrom(sandwich, estfun, bread, kernHAC, weightsAndrews, vcovHAC, bwAndrews, meatHC)
importFrom(methods, is)
Modified: pkg/gmm/R/ategel.R
===================================================================
--- pkg/gmm/R/ategel.R 2019-04-24 14:40:34 UTC (rev 138)
+++ pkg/gmm/R/ategel.R 2019-05-07 19:43:14 UTC (rev 139)
@@ -1,10 +1,10 @@
ATEgel <- function(g, balm, w=NULL, 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"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD","RCUE"),
tol_lam = 1e-9, tol_obj = 1e-9, tol_mom = 1e-9, maxiterlam = 100,
optfct = c("optim", "nlminb"),
- optlam = c("nlminb", "optim", "iter", "Wu", "RCue"),
+ optlam = c("nlminb", "optim", "iter", "Wu"),
data=NULL, Lambdacontrol = list(),
model = TRUE, X = FALSE, Y = FALSE, ...)
{
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2019-04-24 14:40:34 UTC (rev 138)
+++ pkg/gmm/R/gel.R 2019-05-07 19:43:14 UTC (rev 139)
@@ -11,12 +11,13 @@
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
-.rho <- function(x, lamb, derive = 0, type = c("EL", "ET", "CUE", "HD", "ETEL", "ETHD"),
- k = 1)
+.rho <- function(x, lamb, derive = 0, type = c("EL", "ET", "CUE", "HD", "ETEL", "ETHD",
+ "RCUE"), k = 1)
{
type <- match.arg(type)
- lamb <- matrix(lamb, ncol = 1)
- gml <- x%*%lamb*k
+ if (type == "RCUE")
+ type <- "CUE"
+ gml <- x%*%c(lamb)*k
if (derive == 0)
{
if (type == "EL")
@@ -111,7 +112,8 @@
}
return(list(lambda = res$par, convergence = conv,
- obj = mean(.rho(gt,res$par, derive=0,type=type,k=k))))
+ obj = mean(.rho(gt,res$par, derive=0,type=type,k=k)) -
+ .rho(1,0, derive=0,type=type,k=k)))
}
.Wu <- function(gt, tol_lam = 1e-8, maxiterlam = 50, K=1)
@@ -149,56 +151,104 @@
mean(.rho(gt,-M,derive=0,type="EL",k=K))))
}
-.CUE_lam <- function(gt, K=1, control=list(), numSol=FALSE)
+.Wu2 <- function(gt, tol_lam = 1e-8, maxiter = 50, K = 1)
{
- 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)
- }
+ gt <- as.matrix(gt)
+ res <- .Fortran("wu", as.double(gt), as.double(tol_lam),
+ as.integer(maxiter), as.integer(nrow(gt)),
+ as.integer(ncol(gt)), as.double(K),
+ conv=integer(1), obj=double(1),
+ lambda=double(ncol(gt)))
+ list(lambda=res$lambda, convergence=list(convergence=res$conv),
+ obj = res$obj)
+ }
+
+.CUE_lam <- function(gt, K=1)
+ {
+ q <- qr(gt)
+ n <- nrow(gt)
+ l0 <- -qr.coef(q, rep(1,n))
+ conv <- list(convergence=0)
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"),
+.CUE_lamPos <- function(gt, K=1, control=list())
+ {
+ getpt <- function(gt,lam)
+ {
+ gl <- c(gt%*%lam)
+ pt <- 1 + gl
+ pt/sum(pt)
+ }
+ maxit <- ifelse("maxit" %in% names(control),
+ control$maxit, 50)
+ res <-.CUE_lam(gt, K)
+ n <- nrow(gt)
+ i <- 1
+ pt <- getpt(gt, res$lambda)
+ w <- pt<0
+ while (TRUE)
+ {
+ gt2 <- gt[!w,]
+ n1 <- nrow(gt2)
+ if (n1 == n)
+ break
+ res <- try(.CUE_lam(gt2), silent=TRUE)
+ if (i > maxit)
+ return(list(lambda=rep(0,ncol(gt)), obj=0, pt=rep(1/n,n),
+ convergence=list(convergence=1)))
+ if (class(res) == "try-error")
+ return(list(lambda=rep(0,ncol(gt)), obj=0, pt=rep(1/n,n),
+ convergence=list(convergence=2)))
+ pt[!w] <- getpt(gt2, res$lambda)
+ pt[w] <- 0
+ if (all(pt>=0))
+ break
+ i <- i+1
+ w[!w] <- pt[!w]<0
+ }
+ n0 <- n-n1
+ res$obj <- res$obj*n1/n + n0/(2*n)
+ res
+ }
+
+.CUE_lamPos2 <- function(gt, K=1, control=list())
+ {
+ gt <- as.matrix(gt)
+ n <- nrow(gt)
+ q <- ncol(gt)
+ maxit <- ifelse("maxit" %in% names(control),
+ control$maxit, 50)
+ res <- try(.Fortran("lamcuep", as.double(gt),
+ as.integer(n), as.integer(q), as.double(K),
+ as.integer(maxit),conv=integer(1),
+ lam=double(q),pt=double(n),
+ obj=double(1)
+ ), silent=TRUE)
+ if (class(res) == "try-error")
+ return(list(lambda=rep(0,q), obj=0, pt=rep(1/n,n),
+ convergence=list(convergence=3)))
+ list(lambda=res$lam, obj=res$obj, pt=res$pt,
+ convergence=list(convergence=res$conv))
+ }
+
+getLamb <- function(gt, l0, type = c('EL', 'ET', 'CUE', "ETEL", "HD", "ETHD", "RCUE"),
tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7,
- k = 1, method = c("nlminb", "optim", "iter", "Wu", "RCue"),
+ k = 1, method = c("nlminb", "optim", "iter", "Wu"),
control=list())
{
method <- match.arg(method)
+ type <- match.arg(type)
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))
+ return(.Wu2(gt, tol_lam, maxiterlam, k))
if (type == "CUE")
- return(.CUE_lam(gt, k, control, method=="RCue"))
+ return(.CUE_lam(gt, k))
+ if (type == "RCUE")
+ return(.CUE_lamPos2(gt, k, control))
if (method == "iter")
{
if ((type == "ETEL") | (type == "ETHD"))
@@ -266,7 +316,8 @@
res$evaluations, message = res$message)
}
return(list(lambda = l0, convergence = conv, obj =
- mean(.rho(gt,l0,derive=0,type=type,k=k))))
+ mean(.rho(gt,l0,derive=0,type=type,k=k))-
+ .rho(1, 0, type = type, derive = 0, k = k)))
}
smoothG <- function (x, bw = bwAndrews, prewhite = 1, ar.method = "ols",
@@ -309,15 +360,14 @@
return(sx)
}
-
gel <- function(g, x, tet0 = NULL, gradv = NULL, smooth = FALSE,
- type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD","RCUE"),
kernel = c("Truncated", "Bartlett"), bw = bwAndrews,
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", "RCue"), data,
+ optlam = c("nlminb", "optim", "iter", "Wu"), data,
Lambdacontrol = list(),
model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", alpha = NULL,
eqConst = NULL, eqConstFullVcov = FALSE, onlyCoefficients=FALSE, ...)
@@ -353,12 +403,12 @@
}
evalGel <- function(g, x, tet0, gradv = NULL, smooth = FALSE,
- type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD","RCUE"),
kernel = c("Truncated", "Bartlett"), bw = bwAndrews,
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", "RCue"), data,
+ optlam = c("nlminb", "optim", "iter", "Wu"), data,
Lambdacontrol = list(),
model = TRUE, X = FALSE, Y = FALSE, alpha = NULL, ...)
{
@@ -385,7 +435,7 @@
z <- momentEstim(Model_info, ...)
class(z) <- "gel"
return(z)
- }
+ }
.thetf <- function(tet, P, output=c("obj","all"), l0Env)
{
@@ -427,14 +477,16 @@
control=P$Lambdacontrol,
k = P$k1/P$k2, alpha = P$CGEL),silent=TRUE)
}
- if (P$type != "ETHD")
- obj <- mean(.rho(gt, lamb$lambda, type = P$type, derive = 0,
- k = P$k1/P$k2) - .rho(1, 0, type = P$type,
- derive = 0, k = P$k1/P$k2))
- else
- obj <- sum(.rho(gt, lamb$lambda, type = P$type, derive = 0,
- k = P$k1/P$k2) - .rho(1, 0, type = P$type, derive = 0,
- k = P$k1/P$k2))
+ if (P$type == "ETEL")
+ obj <- mean(.rho(gt, lamb$lambda, type = P$type, derive = 0,
+ k = P$k1/P$k2) - .rho(1, 0, type = P$type,
+ derive = 0, k = P$k1/P$k2))
+ else if (P$type == "ETHD")
+ obj <- sum(.rho(gt, lamb$lambda, type = P$type, derive = 0,
+ k = P$k1/P$k2) - .rho(1, 0, type = P$type, derive = 0,
+ k = P$k1/P$k2))
+ else
+ obj <- lamb$obj
assign("l0",lamb$lambda,envir=l0Env)
if(output == "obj")
return(obj)
@@ -454,6 +506,8 @@
eps <- -length(pt)*min(min(pt),0)
pt <- (pt+eps/length(pt))/(1+eps)
}
+ if (type=="RCUE")
+ pt[pt<0] <- 0
###################
conv_moment <- colSums(pt*gt)
conv_pt <- sum(as.numeric(pt))
Modified: pkg/gmm/man/ATEgel.Rd
===================================================================
--- pkg/gmm/man/ATEgel.Rd 2019-04-24 14:40:34 UTC (rev 138)
+++ pkg/gmm/man/ATEgel.Rd 2019-05-07 19:43:14 UTC (rev 139)
@@ -12,10 +12,10 @@
\usage{
ATEgel(g, balm, w=NULL, 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"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD", "RCUE"),
tol_lam = 1e-9, tol_obj = 1e-9, tol_mom = 1e-9, maxiterlam = 100,
optfct = c("optim", "nlminb"),
- optlam = c("nlminb", "optim", "iter", "Wu", "RCue"), data=NULL,
+ optlam = c("nlminb", "optim", "iter", "Wu"), data=NULL,
Lambdacontrol = list(),
model = TRUE, X = FALSE, Y = FALSE, ...)
checkConv(obj, tolConv=1e-4, verbose=TRUE, ...)
@@ -64,31 +64,43 @@
"CUE" for continuous updated estimator, "ETEL" for exponentially
tilted empirical likelihood of Schennach(2007), "HD" for Hellinger
Distance of Kitamura-Otsu-Evdokimov (2013), and "ETHD" for the
- exponentially tilted Hellinger distance of Antoine-Dovonon (2015).}
+ exponentially tilted Hellinger distance of Antoine-Dovonon
+ (2015). "RCUE" is a restricted version of "CUE" in which the
+ probabilities are bounded below by zero. In that case, an analytical
+ Kuhn-Tucker method is used to find the solution.}
-\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches \code{tol_lamb} (see \code{\link{getLamb}}) }
+\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The
+ algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches
+ \code{tol_lamb} (see \code{\link{getLamb}}) }
-\item{maxiterlam}{The algorithm to compute \eqn{\lambda} stops if there is no convergence after "maxiterlam" iterations (see \code{\link{getLamb}}).}
+\item{maxiterlam}{The algorithm to compute \eqn{\lambda} stops if there
+ is no convergence after "maxiterlam" iterations (see
+ \code{\link{getLamb}}).}
-\item{tol_obj}{Tolerance for the gradiant of the objective function to compute \eqn{\lambda} (see \code{\link{getLamb}}).}
+\item{tol_obj}{Tolerance for the gradiant of the objective function to
+ compute \eqn{\lambda} (see \code{\link{getLamb}}).}
\item{optfct}{Algorithm used for the parameter estimates}
-\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{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"}. 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.}
+ \code{type="EL"}. The value of \code{optlam} is ignored for "CUE"
+ because in that case, the analytical solution exists.}
\item{data}{A data.frame or a matrix with column names (Optional). }
-\item{Lambdacontrol}{Controls for the optimization of the vector of Lagrange multipliers used by either \code{\link{optim}}, \code{\link{nlminb}} or \code{\link{constrOptim}}}
+\item{Lambdacontrol}{Controls for the optimization of the vector of
+ Lagrange multipliers used by either \code{\link{optim}},
+ \code{\link{nlminb}} or \code{\link{constrOptim}}}
-\item{model, X, Y}{logicals. If \code{TRUE} the corresponding components of the fit (the model frame, the model matrix, the response) are returned if g is a formula.}
+\item{model, X, Y}{logicals. If \code{TRUE} the corresponding
+ components of the fit (the model frame, the model matrix, the response)
+ are returned if g is a formula.}
\item{verbose}{If TRUE, a summary of the convergence is printed}
Modified: pkg/gmm/man/gel.Rd
===================================================================
--- pkg/gmm/man/gel.Rd 2019-04-24 14:40:34 UTC (rev 138)
+++ pkg/gmm/man/gel.Rd 2019-05-07 19:43:14 UTC (rev 139)
@@ -10,22 +10,22 @@
}
\usage{
gel(g, x, tet0 = NULL, gradv = NULL, smooth = FALSE,
- type = c("EL","ET","CUE","ETEL","HD","ETHD"),
+ type = c("EL","ET","CUE","ETEL","HD","ETHD","RCUE"),
kernel = c("Truncated", "Bartlett"), bw = bwAndrews,
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", "RCue"), 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, onlyCoefficients=FALSE, ...)
evalGel(g, x, tet0, gradv = NULL, smooth = FALSE,
- type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD"),
+ type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD","RCUE"),
kernel = c("Truncated", "Bartlett"), bw = bwAndrews,
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", "RCue"), data, Lambdacontrol = list(),
+ "iter", "Wu"), data, Lambdacontrol = list(),
model = TRUE, X = FALSE, Y = FALSE, alpha = NULL, ...)
}
\arguments{
@@ -46,7 +46,10 @@
"CUE" for continuous updated estimator, "ETEL" for exponentially
tilted empirical likelihood of Schennach(2007), "HD" for Hellinger
Distance of Kitamura-Otsu-Evdokimov (2013), and "ETHD" for the
- exponentially tilted Hellinger distance of Antoine-Dovonon (2015).}
+ exponentially tilted Hellinger distance of Antoine-Dovonon
+ (2015). "RCUE" is a restricted version of "CUE" in which the
+ probabilities are bounded below by zero. In that case, an analytical
+ Kuhn-Tucker method is used to find the solution.}
\item{kernel}{type of kernel used to compute the covariance matrix of the vector of sample moment conditions (see \code{\link{kernHAC}} for more details) and to smooth the moment conditions if "smooth" is set to TRUE. Only two types of kernel are available. The truncated implies a Bartlett kernel for the HAC matrix and the Bartlett implies a Parzen kernel (see Smith 2004).}
@@ -74,11 +77,8 @@
\item{optlam}{Algorithm used to solve for the lagrange multiplier in
\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.}
+ \code{type="EL"}. The value of \code{optlam} is ignored for "CUE"
+ because in that case, the analytical solution exists.}
\item{data}{A data.frame or a matrix with column names (Optional). }
Modified: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd 2019-04-24 14:40:34 UTC (rev 138)
+++ pkg/gmm/man/getLamb.Rd 2019-05-07 19:43:14 UTC (rev 139)
@@ -8,9 +8,9 @@
It computes the vector of Lagrange multipliers, which maximizes the GEL objective function, using an iterative Newton method.
}
\usage{
-getLamb(gt, l0, type = c("EL","ET","CUE", "ETEL", "HD","ETHD"),
+getLamb(gt, l0, type = c("EL","ET","CUE", "ETEL", "HD","ETHD","RCUE"),
tol_lam = 1e-7, maxiterlam = 100,
- tol_obj = 1e-7, k = 1, method = c("nlminb", "optim", "iter", "Wu","RCue"),
+ tol_obj = 1e-7, k = 1, method = c("nlminb", "optim", "iter", "Wu"),
control = list())
}
\arguments{
@@ -20,15 +20,26 @@
\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting,
"CUE" for continuous updated estimator, and "HD" for Hellinger
- Distance. See details for "ETEL" and "ETHD".}
+ Distance. See details for "ETEL" and "ETHD". "RCUE" is a restricted
+ version of "CUE" in which the probabilities are bounded below by
+ zero. In that case, an analytical Kuhn-Tucker method is used to find
+ the solution.}
-\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches \code{tol_lam} }
+\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The
+algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches
+\code{tol_lam} }
-\item{maxiterlam}{The algorithm stops if there is no convergence after "maxiterlam" iterations.}
+\item{maxiterlam}{The algorithm stops if there is no convergence after
+"maxiterlam" iterations.}
-\item{tol_obj}{Tolerance for the gradiant of the objective function. The algorithm returns a non-convergence message if \eqn{\max(|gradiant|)} does not reach \code{tol_obj}. It helps the \code{gel} algorithm to select the right space to look for \eqn{\theta}}
+\item{tol_obj}{Tolerance for the gradiant of the objective function. The
+algorithm returns a non-convergence message if \eqn{\max(|gradiant|)}
+does not reach \code{tol_obj}. It helps the \code{gel} algorithm to
+select the right space to look for \eqn{\theta}}
-\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{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
@@ -37,24 +48,20 @@
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 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.}
+ (2005). The value of \code{method} is ignored for "CUE" because in
+ that case, the analytical solution exists.}
-\item{control}{Controls to send to \code{\link{optim}}, \code{\link{nlminb}} or \code{\link{constrOptim}}}
-}
-\details{
-It solves the problem \eqn{\max_{\lambda} \frac{1}{n}\sum_{t=1}^n
- \rho(gt'\lambda)}. For the type "ETEL", it is only used by
+\item{control}{Controls to send to \code{\link{optim}},
+\code{\link{nlminb}} or \code{\link{constrOptim}}} } \details{ It solves
+the problem \eqn{\max_{\lambda} \frac{1}{n}\sum_{t=1}^n
+\rho(gt'\lambda)}. For the type "ETEL", it is only used by
\code{\link{gel}}. In that case \eqn{\lambda} is obtained by maximizing
\eqn{\frac{1}{n}\sum_{t=1}^n \rho(gt'\lambda)}, using
\eqn{\rho(v)=-\exp{v}} (so ET) and \eqn{\theta} by minimizing the same
equation but with \eqn{\rho(v)-\log{(1-v)}}. To avoid NA's,
\code{\link{constrOptim}} is used with the restriction \eqn{\lambda'g_t
- < 1}. The type "ETHD" is experimental and proposed by Antoine-Dovonon
-(2015). The paper is not yet available.
-}
+< 1}. The type "ETHD" is experimental and proposed by Antoine-Dovonon
+(2015). The paper is not yet available. }
\value{
lambda: A \eqn{q\times 1} vector of Lagrange multipliers which solve the system of equations given above.
Added: pkg/gmm/src/Makevars
===================================================================
--- pkg/gmm/src/Makevars (rev 0)
+++ pkg/gmm/src/Makevars 2019-05-07 19:43:14 UTC (rev 139)
@@ -0,0 +1 @@
+PKG_LIBS=$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -lgomp
\ No newline at end of file
Added: pkg/gmm/src/lambda_met.f
===================================================================
--- pkg/gmm/src/lambda_met.f (rev 0)
+++ pkg/gmm/src/lambda_met.f 2019-05-07 19:43:14 UTC (rev 139)
@@ -0,0 +1,131 @@
+ subroutine prep(gt, lam, n, q, d2)
+ integer n, q, i, info, ip(q)
+ double precision gt(n,q), lam(q), d2(q), dd(q,q)
+ double precision tmp(n), tmp2(n), tmp3(n,q)
+
+ call dgemv('n', n, q, 1.0d0, gt, n, lam, 1, 0.0d0, tmp, 1)
+ tmp = 1/(1+tmp)
+ call dgemv('t', n, q, 1.0d0, gt, n, tmp, 1, 0.0d0, d2, 1)
+ tmp2 = tmp**2
+ do i=1,q
+ tmp3(:,i) = -gt(:,i)*tmp2
+ end do
+ call dgemm('t','n',q,q,n,1.0d0, gt, n, tmp3, n, 0.0d0, dd, q)
+ call dgesv(q, 1, dd, q, ip, d2, q, info)
+ end
+
+
+ subroutine wu(gt, tol, maxit, n, q, k, conv, obj, lam)
+ integer n, q, i, conv, maxit
+ double precision obj, lam(q), tol, gt(n,q), k
+ double precision dif, d2(q), r, tmp(n), tmp2(q)
+
+ lam = 0.0d0
+ i = 1
+ dif = 1.0d0
+ do while (dif>tol .and. i <= maxit)
+ call prep(gt, lam, n, q, d2)
+ dif = maxval(abs(d2))
+ r = 1.0d0
+ do while (r>0)
+ r = 0.0d0
+ tmp2 = lam-d2
+ call dgemv('n', n, q, 1.0d0, gt, n, tmp2, 1, 0.0d0,
+ * tmp, 1)
+ if (minval(tmp)<=-1) then
+ r = r+1
+ end if
+ if (r>0) then
+ d2 = d2/2
+ end if
+ end do
+ lam = tmp2
+ i = i+1
+ end do
+ if (i>=maxit) then
+ lam = 0.0d0
+ conv = 1
+ else
+ lam = -lam
+ conv = 0
+ end if
+ obj = sum(log(1+tmp*k))/n
+ end
+
+ subroutine ols(x, y, n, m, lwork, nrhs, info, coef)
+ integer n, m, lwork, nrhs, info
+ double precision x(n,m), y(n,nrhs), work(lwork)
+ double precision coef(m,nrhs)
+ double precision xtmp(n,m), ytmp(n,nrhs)
+
+ xtmp = x
+ ytmp = y
+ call dgels('n', n, m, nrhs, xtmp, n, ytmp, n, work, -1, info)
+ lwork = min(m*n, int(work(1)))
+ if (info == 0) then
+ call dgels('n',n,m,nrhs,xtmp,n,ytmp,n,work,lwork,info)
+ coef = ytmp(1:m,:)
+ end if
+ end
+
+
+ subroutine lamcue(gt, n, q, k, lam, pt, obj)
+ integer n, q, lwork, info
+ double precision gt(n,q), lam(q), one(n), pt(n), obj, k
+
+ lwork = q*3
+ one = -1.0d0
+ call ols(gt, one, n, q, lwork, 1, info, lam)
+ call dgemv('n', n, q, 1.0d0, gt, n, lam, 1, 0.0d0,
+ * pt, 1)
+ pt = pt*k
+ obj = sum(-pt-(pt**2)/2)/n
+ pt = 1+pt
+ pt = pt/sum(pt)
+ end
+
+
+ subroutine lamcuep(gt, n, q, k, maxit, conv, lam, pt, obj)
+ integer n, q, i, maxit, n0, n1, conv, ind(n)
+ double precision gt(n,q), lam(q), pt(n), obj, pt0(n), k
+ integer wi(n), wni(n)
+ logical w(n)
+
+ call lamcue(gt, n, q, k, lam, pt, obj)
+ ind = (/ (i, i=1,n) /)
+ i = 1
+ conv = 0
+ w = (pt < 0)
+ do while (.not. all(pt >= 0))
+ n0 = count(w)
+ n1 = n-n0
+ wi(1:n1) = pack(ind, .not. w)
+ wni(1:n0) = pack(ind, w)
+ if (n1 < (q+1)) then
+ pt = 1.0d0/n
+ conv = 2
+ obj = 0.0d0
+ lam = 0.0d0
+ exit
+ end if
+ if (i > maxit) then
+ pt = 1.0d0/n
+ conv = 1
+ obj = 0.0d0
+ lam = 0.0d0
+ exit
+ end if
+ call lamcue(gt(wi(1:n1),:), n1, q, k, lam,
+ * pt0(1:n1), obj)
+ pt(wi(1:n1)) = pt0(1:n1)
+ pt(wni(1:n0)) = 0.0d0
+ i = i+1
+ w = w .or. (pt<0)
+ end do
+ if (conv == 0) then
+ obj = obj*n1/n + dble(n0)/(2*n)
+ end if
+ end
+
+
+
More information about the Gmm-commits
mailing list