[Gmm-commits] r80 - in pkg/gmm: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 22 21:53:59 CEST 2015
Author: chaussep
Date: 2015-09-22 21:53:58 +0200 (Tue, 22 Sep 2015)
New Revision: 80
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NEWS
pkg/gmm/R/gel.R
pkg/gmm/R/getModel.R
pkg/gmm/R/momentEstim.R
pkg/gmm/R/specTest.R
pkg/gmm/man/gel.Rd
pkg/gmm/man/getLamb.Rd
Log:
Added HD and ETHD to gel and fix LR test
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/DESCRIPTION 2015-09-22 19:53:58 UTC (rev 80)
@@ -1,6 +1,6 @@
Package: gmm
Version: 1.5-3
-Date: 2015-08-25
+Date: 2015-09-25
Title: Generalized Method of Moments and Generalized Empirical
Likelihood
Author: Pierre Chausse <pchausse at uwaterloo.ca>
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/NEWS 2015-09-22 19:53:58 UTC (rev 80)
@@ -1,6 +1,8 @@
Changes in version 1.5-3
o Fixed a typo in FinRes.R file. It was preventing to compute the proper vcov matrix for a very special case (fixed weights)
+o Added Helinger distance and Exponentially tilted Helinger estimator to gel()
+o Fixed the LR test of ETEL in gel()
Changes in version 1.5-2
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/R/gel.R 2015-09-22 19:53:58 UTC (rev 80)
@@ -11,7 +11,7 @@
# 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"), k = 1)
+.rho <- function(x, lamb, derive = 0, type = c("EL", "ET", "CUE", "HD", "ETEL", "ETHD"), k = 1)
{
type <- match.arg(type)
@@ -28,7 +28,21 @@
if (type == "ET")
rhomat <- -exp(gml)
if (type == "CUE")
- rhomat <- -gml -0.5*gml^2
+ rhomat <- -gml -0.5*gml^2
+ if (type == "HD")
+ rhomat <- -2/(1-gml/2)
+ if (type == "ETEL")
+ {
+ w <- -exp(gml)
+ w <- w/sum(w)
+ rhomat <- -log(w*NROW(gml))
+ }
+ if (type == "ETHD")
+ {
+ w <- -exp(gml)
+ w <- w/sum(w)
+ rhomat <- (sqrt(w)-1/sqrt(NROW(gml)))^2
+ }
}
if (derive==1)
{
@@ -37,23 +51,29 @@
if (type == "ET")
rhomat <- -exp(gml)
if (type == "CUE")
- rhomat <- -1 - gml
+ rhomat <- -1 - gml
+ if (type == "HD")
+ rhomat <- -1/((1-gml/2)^2)
+ if (any(type == c("ETEL","ETHD")))
+ rhomat <- NULL
}
if (derive==2)
{
if (type == "EL")
- rhomat <- -1/(1 - gml)^2
-
+ rhomat <- -1/(1 - gml)^2
if (type == "ET")
- rhomat <- -exp(gml)
-
+ rhomat <- -exp(gml)
if (type == "CUE")
- rhomat <- -rep(1,nrow(x))
+ rhomat <- -rep(1,nrow(x))
+ if (type == "HD")
+ rhomat <- -1/((1-gml/2)^3)
+ if (any(type == c("ETEL","ETHD")))
+ rhomat <- NULL
}
return(c(rhomat))
}
-.getCgelLam <- function(gt, l0, type = c('EL', 'ET', 'CUE'), method = c("nlminb", "optim", "constrOptim"), control=list(), k = 1, alpha = 0.01)
+.getCgelLam <- function(gt, l0, type = c('EL', 'ET', 'CUE', "HD"), method = c("nlminb", "optim", "constrOptim"), control=list(), k = 1, alpha = 0.01)
{
type <- match.arg(type)
method <- match.arg(method)
@@ -90,7 +110,7 @@
obj = mean(.rho(gt,res$par, derive=0,type=type,k=k))))
}
-getLamb <- function(gt, l0, type = c('EL', 'ET', 'CUE', "ETEL"), tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7,
+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())
{
method <- match.arg(method)
@@ -99,7 +119,7 @@
if (method == "iter")
{
- if (type == "ETEL")
+ if ((type == "ETEL") | (type == "ETHD"))
type = "ET"
for (i in 1:maxiterlam)
{
@@ -124,43 +144,34 @@
}
else
{
- fct <- function(l,X)
+ fct <- function(l,X,type,k)
{
r0 <- .rho(X,l,derive=0,type=type,k=k)
-mean(r0)
}
- Dfct <- function(l,X)
+ Dfct <- function(l,X,type,k)
{
r1 <- .rho(X,l,derive=1,type=type,k=k)
-colMeans(r1*X)
}
- DDfct <- function(l,X)
+ DDfct <- function(l,X,type,k)
{
r2 <- .rho(X,l,derive=2,type=type,k=k)
-t(X*r2)%*%X/nrow(X)
- }
- if (type == "ETEL")
- {
+ }
+ if ((type == "ETEL")|(type=="ETHD"))
type = "ET"
- ci <- -rep(1,nrow(gt))
- res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt)
- }
- else
- {
- if (method=="optim")
+ if (method=="optim")
{
- if (type != "EL")
- res <- optim(rep(0,ncol(gt)),fct,gr=Dfct,X=gt,method="BFGS",control=control)
- else
- {
+ if (type != "EL"){
+ res <- optim(rep(0,ncol(gt)),fct,gr=Dfct,X=gt,type=type,k=k,method="BFGS",control=control)
+ } else {
ci <- -rep(1,nrow(gt))
- res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt)
- }
+ res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt,type=type,k=k)
}
- else
- res <- nlminb(rep(0,ncol(gt)), fct, gradient = Dfct, hessian = DDfct, X = gt, control = control)
- }
-
+ } else {
+ res <- nlminb(rep(0,ncol(gt)), fct, gradient = Dfct, hessian = DDfct, X = gt, type=type, k=k, control = control)
+ }
l0 <- res$par
if (method == "optim" | method == "constrOptim")
conv <- list(convergence = res$convergence, counts = res$counts, message = res$message)
@@ -212,7 +223,7 @@
}
-gel <- function(g, x, tet0, gradv = NULL, smooth = FALSE, type = c("EL", "ET", "CUE", "ETEL"),
+gel <- function(g, x, tet0, gradv = NULL, smooth = FALSE, type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD"),
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"),
@@ -260,21 +271,25 @@
tol_obj = P$tol_obj, k = P$k1/P$k2, control = P$Lambdacontrol, method = "optim")
}
else
- lamb <- getLamb(gt, l0, type = P$type, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam,
+ lamb <- getLamb(gt, l0, type = P$type, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam,
tol_obj = P$tol_obj, k = P$k1/P$k2, control = P$Lambdacontrol, method = P$optlam)
}
else
{
- if (P$type=="ETEL")
- stop("CGEL not implemented for ETEL")
+ if ((P$type=="ETEL")|(P$type=="ETHD"))
+ stop("CGEL not implemented for ETEL nor for ETHD")
lamb <- try(.getCgelLam(gt, l0, type = P$type, method = "nlminb", control=P$Lambdacontrol,
k = P$k1/P$k2, alpha = P$CGEL),silent=TRUE)
if (class(lamb) == "try-error")
lamb <- try(.getCgelLam(gt, l0, type = P$type, method = "constrOptim", control=P$Lambdacontrol,
k = P$k1/P$k2, alpha = P$CGEL),silent=TRUE)
}
-
- obj <- mean(.rho(gt, lamb$lambda, type = P$typet, derive = 0, k = P$k1/P$k2))
+ 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))
assign("l0",lamb$lambda,envir=l0Env)
if(output == "obj")
return(obj)
@@ -282,4 +297,26 @@
return(list(obj = obj, lambda = lamb, gt = gt))
}
+.getImpProb <- function(gt, lambda, type, k1, k2)
+ {
+ if ((type == "ETEL")|(type=="ETHD"))
+ type <- "ET"
+ n <- NROW(gt)
+ pt <- -.rho(gt, lambda, type = type, derive = 1, k = k1/k2)/n
+ # Making sure pt>0
+ if (type=="CUE")
+ {
+ eps <- -length(pt)*min(min(pt),0)
+ pt <- (pt+eps/length(pt))/(1+eps)
+ }
+ ###################
+ conv_moment <- colSums(pt*gt)
+ conv_pt <- sum(as.numeric(pt))
+ pt <- pt/sum(pt)
+ attr(pt, "conv_moment") <- conv_moment
+ attr(pt, "conv_pt") <- conv_pt
+ pt
+ }
+
+
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/R/getModel.R 2015-09-22 19:53:58 UTC (rev 80)
@@ -223,16 +223,6 @@
{
P <- object
- if (P$type == "ETEL")
- {
- P$typel <- "ET"
- P$typet <- "EL"
- }
- else
- {
- P$typel <- P$type
- P$typet <- P$type
- }
if(P$optfct == "optim" | P$optfct == "nlminb")
P$k <- length(P$tet0)
else
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/R/momentEstim.R 2015-09-22 19:53:58 UTC (rev 80)
@@ -719,29 +719,18 @@
rlamb <- All$lambda
z <- list(coefficients = res$par, lambda = rlamb$lambda, conv_lambda = rlamb$conv, conv_par = res$convergence, dat=P$dat)
- rho1 <- .rho(gt, z$lambda, derive = 1, type = P$typel, k = P$k1/P$k2)
- z$foc_lambda <- crossprod(colMeans(rho1*gt))
z$type <- P$type
z$gt <- gt
- rhom <- .rho(z$gt, z$lambda, type = P$typet, k = P$k1/P$k2)
- z$pt <- -.rho(z$gt, z$lambda, type = P$typet, derive = 1, k = P$k1/P$k2)/n
- # Making sure pt>0
- if (P$type=="CUE")
- {
- eps <- -length(z$pt)*min(min(z$pt),0)
- z$pt <- (z$pt+eps/length(z$pt))/(1+eps)
- }
- ###################
+ pt <- .getImpProb(z$gt, z$lambda, P$type, P$k1, P$k2)
+ z$pt <- c(pt)
+ z$conv_moment <- attr(pt, "conv_moment")
+ z$conv_pt <- attr(pt, "conv_pt")
+ z$objective <- All$obj
- z$conv_moment <- colSums(z$pt*z$gt)
- z$conv_pt <- sum(as.numeric(z$pt))
- z$objective <- sum(rhom - .rho(1, 0, type = P$typet))/n
-
namex <- colnames(P$dat$x[, (P$dat$ny+1):(P$dat$ny+P$dat$k)])
nameh <- colnames(P$dat$x[, (P$dat$ny+P$dat$k+1):(P$dat$ny+P$dat$k+P$dat$nh)])
-
if (P$dat$ny > 1)
{
namey <- colnames(P$dat$x[, 1:P$dat$ny])
@@ -825,25 +814,15 @@
rlamb <- All$lambda
z <- list(coefficients = res$par, lambda = rlamb$lambda, conv_lambda = rlamb$conv, conv_par = res$convergence, dat=P$dat)
- rho1 <- .rho(gt, z$lambda, derive = 1, type = P$typel, k = P$k1/P$k2)
- z$foc_lambda <- crossprod(colMeans(rho1*gt))
z$type <- P$type
z$gt <- gt
- rhom <- .rho(z$gt, z$lambda, type = P$typet, k = P$k1/P$k2)
- z$pt <- -.rho(z$gt, z$lambda, type = P$typel, derive = 1, k = P$k1/P$k2)/n
+ pt <- .getImpProb(z$gt, z$lambda, P$type, P$k1, P$k2)
+ z$pt <- c(pt)
+ z$conv_moment <- attr(pt, "conv_moment")
+ z$conv_pt <- attr(pt, "conv_pt")
+ z$objective <- All$obj
-# making sure pt>0
- if (P$type=="CUE")
- {
- eps <- -length(z$pt)*min(min(z$pt),0)
- z$pt <- (z$pt+eps/length(pt))/(1+eps)
- }
-##################
- z$conv_moment <- colSums(as.numeric(z$pt)*z$gt)
- z$conv_pt <- sum(as.numeric(z$pt))
- z$objective <- sum(as.numeric(rhom) - .rho(1, 0, type = P$typet, k = P$k1/P$k2))/n
-
if(P$gradvf)
G <- P$gradv(z$coefficients, x)
else
Modified: pkg/gmm/R/specTest.R
===================================================================
--- pkg/gmm/R/specTest.R 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/R/specTest.R 2015-09-22 19:53:58 UTC (rev 80)
@@ -41,6 +41,8 @@
khat <- x$khat
gbar <- colMeans(x$gt)
LR_test <- 2*x$objective*n*x$k2/(x$bwVal*x$k1^2)
+ if (x$type == "ETHD")
+ LR_test <- LR_test*2
LM_test <- n*crossprod(x$lambda, crossprod(khat, x$lambda))/(x$bwVal^2)
J_test <- n*crossprod(gbar, solve(khat, gbar))/(x$k1^2)
test <- c(LR_test, LM_test, J_test)
Modified: pkg/gmm/man/gel.Rd
===================================================================
--- pkg/gmm/man/gel.Rd 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/man/gel.Rd 2015-09-22 19:53:58 UTC (rev 80)
@@ -8,7 +8,7 @@
Function to estimate a vector of parameters based on moment conditions using the GEL method as presented by Newey-Smith(2004) and Anatolyev(2005).
}
\usage{
-gel(g, x, tet0, gradv = NULL, smooth = FALSE, type = c("EL","ET","CUE","ETEL"),
+gel(g, x, tet0, gradv = NULL, smooth = FALSE, type = c("EL","ET","CUE","ETEL","HD","ETHD"),
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,
@@ -27,7 +27,11 @@
\item{smooth}{If set to TRUE, the moment function is smoothed as proposed by Kitamura(1997)}
-\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting, "CUE" for continuous updated estimator and "ETEL" for exponentially tilted empirical likelihood of Schennach(2007).}
+\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting,
+ "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).}
\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).}
@@ -123,7 +127,6 @@
}
-
\references{
Anatolyev, S. (2005), GMM, GEL, Serial Correlation, and Asymptotic Bias. \emph{Econometrica}, \bold{73}, 983-1002.
@@ -133,6 +136,10 @@
Kitamura, Yuichi (1997), Empirical Likelihood Methods With Weakly Dependent Processes.
\emph{The Annals of Statistics}, \bold{25}, 2084-2102.
+Kitamura, Y. and Otsu, T. and Evdokimov, K. (2013), Robustness,
+Infinitesimal Neighborhoods and Moment Restrictions.
+\emph{Econometrica}, \bold{81}, 1185-1201.
+
Newey, W.K. and Smith, R.J. (2004), Higher Order Properties of GMM and
Generalized Empirical Likelihood Estimators. \emph{Econometrica}, \bold{72}, 219-255.
Modified: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd 2015-08-25 20:26:33 UTC (rev 79)
+++ pkg/gmm/man/getLamb.Rd 2015-09-22 19:53:58 UTC (rev 80)
@@ -8,7 +8,7 @@
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"), tol_lam = 1e-7, maxiterlam = 100,
+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"), control = list())
}
\arguments{
@@ -16,7 +16,9 @@
\item{l0}{Vector of starting values for lambda}
-\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting and "CUE" for continuous updated estimator. See details for "ETEL".}
+\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".}
\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches \code{tol_lam} }
@@ -31,7 +33,15 @@
\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}.
+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.
}
\value{
More information about the Gmm-commits
mailing list