[Gmm-commits] r44 - in pkg/gmm: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 1 00:00:48 CET 2011
Author: chaussep
Date: 2011-12-01 00:00:47 +0100 (Thu, 01 Dec 2011)
New Revision: 44
Removed:
pkg/gmm/man/rho.Rd
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NAMESPACE
pkg/gmm/NEWS
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/gel.Rd
pkg/gmm/man/getLamb.Rd
pkg/gmm/man/summary.Rd
Log:
Many change especially for computing lamda with GEL
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/DESCRIPTION 2011-11-30 23:00:47 UTC (rev 44)
@@ -1,6 +1,6 @@
Package: gmm
-Version: 1.3-9
-Date: 2011-07-19
+Version: 1.4-0
+Date: 2011-11-30
Title: Generalized Method of Moments and Generalized Empirical
Likelihood
Author: Pierre Chausse <pchausse at uwaterloo.ca>
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/NAMESPACE 2011-11-30 23:00:47 UTC (rev 44)
@@ -2,7 +2,7 @@
importFrom(sandwich, estfun, bread)
-export(gmm,summary.gmm,rho,smoothG,getDat,summary.gel,getLamb,gel, estfun.gmmFct, estfun.gmm, estfun.gel, bread.gel, bread.gmm,
+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,
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/NEWS 2011-11-30 23:00:47 UTC (rev 44)
@@ -1,4 +1,4 @@
-Changes in version 1.3-9
+Changes in version 1.4-0
o The method for GMM-CUE has been modified. Before, the weights for the kernel estimation of the weighting matrix
were flexible inside the optimizer which was making the algorithm long and unstable. It is now fixed using either the starting
@@ -11,6 +11,7 @@
o Something the GMM converge to weird values which is sometimes caused by bad first step estimates used to compute the weighting matrix.
Summary() prints the initial values to have more infortmation when convergence fails.
o There was a bug in specTest() for GEL. The degrees of freedom for the J-test were wrong. It is fixed.
+o For GEL type CUE the implied probabilities are computed according to Antoine, Bonnal, and Renault (2007) which solves the problem of negative probabilities.
Changes in version 1.3-8
Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/R/Methods.gel.R 2011-11-30 23:00:47 UTC (rev 44)
@@ -68,6 +68,9 @@
cat("Lambdas:\n")
print.default(format(coef(x, lambda = TRUE), digits = digits),
print.gap = 2, quote = FALSE)
+ cat("\n")
+ cat("Convergence code for the coefficients: ", x$conv_par,"\n")
+ cat("Convergence code for Lambda: ", x$conv_lambda$convergence,"\n")
invisible(x)
}
@@ -90,7 +93,7 @@
print.gap = 2, quote = FALSE)
cat("\nConvergence code for the coefficients: ", x$conv_par, "\n")
- cat("\nConvergence code for the lambdas: ", x$conv_lambda, "\n")
+ cat("\nConvergence code for the lambdas: ", x$conv_lambda$convergence, "\n")
invisible(x)
}
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/R/gel.R 2011-11-30 23:00:47 UTC (rev 44)
@@ -11,34 +11,18 @@
# 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"), drop = TRUE, k = 1)
+.rho <- function(x, lamb, derive = 0, type = c("EL", "ET", "CUE"), k = 1)
{
type <- match.arg(type)
lamb <- matrix(lamb, ncol = 1)
gml <- x%*%lamb*k
- ch <- 0
if (derive == 0)
{
if (type == "EL")
- {
- ch <- sum(gml >= 1)
- if (drop)
- {
- gml <- (gml < 1)*gml
- rhomat <- log(1 - gml)
- }
- else
- {
- if (ch > 0)
- rhomat <- NaN
- else
- rhomat <- log(1 - gml)
- }
- }
+ rhomat <- log(1 - gml)
if (type == "ET")
rhomat <- -exp(gml)
-
if (type == "CUE")
rhomat <- -gml -0.5*gml^2
}
@@ -46,10 +30,8 @@
{
if (type == "EL")
rhomat <- -1/(1 - gml)
-
if (type == "ET")
rhomat <- -exp(gml)
-
if (type == "CUE")
rhomat <- -1 - gml
}
@@ -62,81 +44,80 @@
rhomat <- -exp(gml)
if (type == "CUE")
- rhomat <- -1
+ rhomat <- -rep(1,nrow(x))
}
- rhom <-list(ch = ch, rhomat = rhomat)
- return(rhom)
+ return(c(rhomat))
}
-getLamb <- function(g, tet, x, type = c('EL', 'ET', 'CUE'), tol_lam = 1e-12, maxiterlam = 1000, tol_obj = 1e-7, k = 1)
+
+getLamb <- function(gt, type = c('EL', 'ET', 'CUE'), tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7, k = 1,
+ method = c("nlminb", "optim", "iter"), control=list())
{
- type <- match.arg(type)
- gt <- g(tet, x)
+ method <- match.arg(method)
+ if(is.null(dim(gt)))
+ gt <- matrix(gt,ncol=1)
+ l0 <- rep(0,ncol(gt))
- n <- nrow(gt)
- tol_cond=1e-12
- gb <- colMeans(gt)
- khat <- crossprod(gt)/n
- lamb0 <- -solve(khat,gb)
-
- conv_mes <- "Normal convergence"
- singular <-0
- crit <-1e30
- crit0 <- crit
- dcrit <- 10
- dgblam <- -10
- gblam0 <- NULL
-
- j <- 1
- while ((crit > tol_lam*( 1 + sqrt( crossprod(lamb0) ) ) ) & (j <= maxiterlam))
- {
- rho2 <- as.numeric(rho(gt, lamb0, derive = 2, type = type, k = k)$rhomat)
- rho1 <- as.numeric(rho(gt, lamb0, derive = 1, type = type, k = k)$rhomat)
- gblam <- colMeans(rho1*gt)
- klam <- crossprod(rho2*gt, gt)/n
- chklam <- sum(abs(klam))
- if (!is.null(gblam0))
- dgblam <- crossprod(gblam) - crossprod(gblam0)
-
- #if (is.na(chklam) | chklam == 0 | chklam == Inf | dgblam > 0 | dgblam == Inf | is.na(dgblam) | dcrit < 0)
- if (is.na(chklam) | chklam == 0 | chklam == Inf | dgblam == Inf | is.na(dgblam))
+ if (method == "iter")
+ {
+ for (i in 1:maxiterlam)
{
- lamb1 <- rep(0, length(lamb0))
- crit <- 0
- singular=2
- conv_mes <- "The algorithm produced singular system, NaN or Inf"
- }
- else
- {
- if (rcond(klam) > tol_cond)
+ r1 <- .rho(gt,l0,derive=1,type=type,k=k)
+ r2 <- .rho(gt,l0,derive=2,type=type,k=k)
+ F <- -colMeans(r1*gt)
+ J <- crossprod(r2*gt,gt)
+ if (sum(abs(F))<tol_obj)
{
- lamb1 <- lamb0 - solve(klam, gblam)
- crit <- sqrt(crossprod(lamb0 - lamb1))
- lamb0 <- lamb1
+ conv <- "Tolerance for the FOC reached"
+ break
}
- else
+ P <- solve(J,F)
+ if (sum(abs(P))<tol_lam)
{
- lamb1 <- rep(0 , length(lamb0))
- crit <- 0
- singular <- 2
- conv_mes <- "The algorithm produced singular system"
+ cov <- "Tolerance on lambda reached"
+ break
}
+ l0 <- l0 + P
+ cov <- "maxiterlam reached"
}
- gblam0 <- gblam
- j <- j + 1
- dcrit<- crit0 - crit
- crit0 <- crit
}
- z <- list("lambda" = lamb1, singular = singular, conv_mes = conv_mes)
- if (j > maxiterlam)
+ else
{
- singular <- 1
- conv_mes <- "No convergence after 'maxiterlam' iterations"
- z$singular <- singular
- stop("Maxiterlam reached.\n Increase it, try other starting values \n or use the option optlam=\"numeric\".")
+ fct <- function(l,X)
+ {
+ r0 <- .rho(X,l,derive=0,type=type,k=k)
+ -mean(r0)
+ }
+ Dfct <- function(l,X)
+ {
+ r1 <- .rho(X,l,derive=1,type=type,k=k)
+ -colMeans(r1*X)
+ }
+ DDfct <- function(l,X)
+ {
+ r2 <- .rho(X,l,derive=2,type=type,k=k)
+ -t(X*r2)%*%X/nrow(X)
+ }
+ if (method=="optim")
+ {
+ if (type != "EL")
+ res <- optim(rep(0,ncol(gt)),fct,gr=Dfct,X=gt,method="B",control=control)
+ else
+ {
+ ci <- -rep(1,nrow(gt))
+ res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt)
+ }
+ }
+ else
+ res <- nlminb(rep(0,ncol(gt)), fct, gradient = Dfct, hessian = DDfct, X = gt, control = control)
+
+ l0 <- res$par
+ if (method == "optim" | method == "constrOptim")
+ conv <- list(convergence = res$convergence, counts = res$counts, message = res$message)
+ if(method == "nlminb")
+ conv <- list(convergence = res$convergence, counts = res$evaluations, message = res$message)
}
- z$obj <- crossprod(gblam)
- return(z)
+ return(list(lambda = l0, convergence = conv, obj = mean(.rho(gt,l0,derive=0,type=type,k=k))))
}
smoothG <- function (x, bw = bwAndrews, prewhite = 1, ar.method = "ols", weights = weightsAndrews,
@@ -185,7 +166,7 @@
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("iter", "numeric"), model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", ...)
+ optlam = c("nlminb", "optim", "iter"), controlLam = list(), model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", ...)
{
type <- match.arg(type)
@@ -199,7 +180,8 @@
kernel = kernel, bw = bw, approx = approx, prewhite = prewhite, ar.method = ar.method,
tol_weights = tol_weights, tol_lam = tol_lam, tol_obj = tol_obj, tol_mom = tol_mom,
maxiterlam = maxiterlam, constraint = constraint, optfct = optfct, weights = weights,
- optlam = optlam, model = model, X = X, Y = Y, TypeGel = TypeGel, call = match.call())
+ optlam = optlam, model = model, X = X, Y = Y, TypeGel = TypeGel, call = match.call(),
+ controlLam = controlLam)
class(all_args)<-TypeGel
Model_info<-getModel(all_args)
@@ -212,50 +194,22 @@
.thetf <- function(tet, P)
{
- if(!is.null(P$gform))
- {
- dat <- P$dat
- x <- dat$x
- }
+ gt <- P$g(tet, P$dat)
+ if (P$optlam != "optim" & P$type == "EL")
+ {
+ lamb <- try(getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2,
+ control = P$controlLam, method = P$optlam)$lambda, silent = TRUE)
+ if(class(lamb) == "try-error")
+ lamb <- getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2,
+ control = P$controlLam, method = "optim")$lambda
+ }
else
- x <- P$x
+ lamb <- getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2,
+ control = P$controlLam, method = P$optlam)$lambda
- if (P$optlam == "iter")
- {
- lamblist <- getLamb(P$g, tet, x, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2)
- lamb <- lamblist$lambda
- gt <- P$g(tet, x)
- pt <- -rho(gt, lamb, type = P$typet, derive = 1, k = P$k1/P$k2)$rhomat/nrow(gt)
- checkmom <- sum(as.numeric(pt)*gt)
- if (lamblist$singular == 0)
- p <- sum(rho(gt, lamb, type = P$typet, k = P$k1/P$k2)$rhomat) + abs(checkmom)/P$tol_mom
- if (lamblist$singular == 1)
- p <- sum(rho(gt, lamb, type = P$typet, k = P$k1/P$k2)$rhomat) + abs(checkmom)/P$tol_mom + lamblist$obj/P$tol_mom
- if (lamblist$singular == 2)
- p <- 1e50*proc.time()[3]
- }
- else
- {
- gt <- P$g(tet, x)
- rhofct <- function(lamb)
- {
- rhof <- -sum(rho(gt, lamb, type = P$typel, k = P$k1/P$k2)$rhomat)
- return(rhof)
- }
- if (ncol(gt) > 1)
- rlamb <- optim(rep(0, ncol(gt)), rhofct, control = list(maxit = 1000,parscale=rep(.01,ncol(gt))),method="B")
- else
- {
- rlamb <- optimize(rhofct, c(-1,1))
- rlamb$par <- rlamb$minimum
- rlamb$value <- rlamb$objective
- }
- lamb <- rlamb$par
- pt <- -rho(gt, lamb, type = P$typet, derive = 1, k = P$k1/P$k2)$rhomat/nrow(gt)
- checkmom <- sum(as.numeric(pt)*gt)
- p <- -rlamb$value + (checkmom)^2/P$tol_mom + (sum(as.numeric(pt)) - 1)^2/P$tol_mom
- }
- return(p)
+ obj <- mean(.rho(gt, lamb, type = P$typet, derive = 0, k = P$k1/P$k2))
+
+ return(obj)
}
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/R/getModel.R 2011-11-30 23:00:47 UTC (rev 44)
@@ -120,10 +120,13 @@
{
clname <- paste(class(P), ".modFormula", sep = "")
dat <- getDat(P$g, P$x)
- x <- dat$x
- g <- function(tet, x, ny = dat$ny, nh = dat$nh, k = dat$k)
+ g <- function(tet, dat)
{
+ x <- dat$x
+ ny <- dat$ny
+ nh <- dat$nh
+ k <- dat$k
tet <- matrix(tet, ncol = k)
e <- x[,1:ny] - x[, (ny+1):(ny+k)]%*%t(tet)
gt <- e*x[, ny+k+1]
@@ -137,10 +140,18 @@
return(gt)
}
- gradv <- function(tet, x, ny = dat$ny, nh = dat$nh, k = dat$k)
+ gradv <- function(tet, dat, pt = NULL)
{
+ x <- dat$x
+ ny <- dat$ny
+ nh <- dat$nh
+ k <- dat$k
tet <- matrix(tet, ncol = k)
- dgb <- -(t(x[,(ny+k+1):(ny+k+nh)])%*%x[,(ny+1):(ny+k)])%x%diag(rep(1, ny))/nrow(x)
+ if (is.null(pt))
+ dgb <- -(t(x[,(ny+k+1):(ny+k+nh)])%*%x[,(ny+1):(ny+k)])%x%diag(rep(1, ny))/nrow(x)
+ else
+ dgb <- -(t(c(pt)*x[,(ny+k+1):(ny+k+nh)])%*%x[,(ny+1):(ny+k)])%x%diag(rep(1, ny))
+
return(dgb)
}
P$dat <- dat
@@ -150,7 +161,7 @@
}
else
{
- x <- P$x
+ P$dat <- P$x
clname <- paste(class(P), ".mod", sep = "")
P$gform <- NULL
if (!is.function(object$gradv))
@@ -179,9 +190,8 @@
P$k2 <- 2/3
}
P$g1 <- P$g
- rgmm <- gmm(P$g, x, P$tet0, wmatrix = "ident")
+ rgmm <- gmm(P$g, P$dat, P$tet0, wmatrix = "ident")
-
P$bwVal <- P$bw(centeredGt <- lm(P$g(rgmm$coefficients, x)~1), kernel = P$wkernel, prewhite = P$prewhite,
ar.method = P$ar.method, approx = P$approx)
P$w <- smoothG(residuals(centeredGt), bw = P$bwVal)$kern_weights
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/R/gmm.R 2011-11-30 23:00:47 UTC (rev 44)
@@ -212,7 +212,7 @@
return(obj)
}
-.Gf <- function(thet, x, g)
+.Gf <- function(thet, x, g, pt = NULL)
{
myenv <- new.env()
assign('x', x, envir = myenv)
@@ -220,7 +220,11 @@
barg <- function(x, thet)
{
gt <- g(thet, x)
- gbar <- as.vector(colMeans(gt))
+ if (is.null(pt))
+ gbar <- as.vector(colMeans(gt))
+ else
+ gbar <- as.vector(colSums(c(pt)*gt))
+
return(gbar)
}
G <- attr(numericDeriv(quote(barg(x, thet)), "thet", myenv), "gradient")
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/R/momentEstim.R 2011-11-30 23:00:47 UTC (rev 44)
@@ -695,10 +695,7 @@
{
P <- object
g <- P$g
- dat <- getDat(P$gform, P$x)
- x <- dat$x
- n <- nrow(x)
-
+ n <- nrow(P$dat$x)
if (!P$constraint)
{
if (P$optfct == "optim")
@@ -713,66 +710,68 @@
res$convergence <- "There is no convergence code for optimize"
}
}
-
if(P$constraint)
res <- constrOptim(P$tet0, .thetf, grad = NULL, P = P, ...)
+ gt <- P$g(res$par, P$dat)
+
- if (P$optlam == "iter")
- {
- rlamb <- getLamb(P$g, res$par, x, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2)
- z <- list(coefficients = res$par, lambda = rlamb$lam, conv_lambda = rlamb$conv_mes, conv_par = res$convergence)
- z$foc_lambda <- rlamb$obj
- }
+ if (P$optlam != "optim" & P$type == "EL")
+ {
+ rlamb <- try(getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2,
+ control = P$controlLam, method = P$optlam), silent = TRUE)
+ if(class(rlamb) == "try-error")
+ rlamb <- getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2,
+ control = P$controlLam, method = "optim")
+ }
+ else
+ rlamb <- getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2,
+ control = P$controlLam, method = P$optlam)
- if (P$optlam == "numeric")
- {
- gt <- P$g(res$par, x)
- rhofct <- function(lamb)
- {
- rhof <- -sum(rho(gt, lamb, type = P$typel, k = P$k1/P$k2)$rhomat)
- return(rhof)
- }
- rlamb <- optim(rep(0, ncol(gt)), rhofct, control = list(maxit = 1000,parscale=rep(.01,ncol(gt))), method="B")
- z <- list(coefficients = res$par, conv_par = res$convergence, lambda = rlamb$par)
- z$conv_lambda = paste("Lambda by optim. Conv. code = ", rlamb$convergence, sep = "")
- rho1 <- as.numeric(rho(gt, z$lambda, derive = 1, type = P$typel, k = P$k1/P$k2)$rhomat)
- z$foc_lambda <- crossprod(colMeans(rho1*gt))
- }
+
+ 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 <- P$g(z$coefficients, x)
- 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)$rhomat/n
- z$conv_moment <- colSums(as.numeric(z$pt)*z$gt)
+ 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)
+ }
+ ###################
+
+ z$conv_moment <- colSums(z$pt*z$gt)
z$conv_pt <- sum(as.numeric(z$pt))
- z$objective <- sum(as.numeric(rhom$rhomat) - rho(1, 0, type = P$typet)$rhomat)/n
+ z$objective <- sum(rhom - .rho(1, 0, type = P$typet))/n
- namex <- colnames(dat$x[, (dat$ny+1):(dat$ny+dat$k)])
- nameh <- colnames(dat$x[, (dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
+ 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 (dat$ny > 1)
+
+ if (P$dat$ny > 1)
{
- namey <- colnames(dat$x[, 1:dat$ny])
- names(z$coefficients) <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
- colnames(z$gt) <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny,dat$nh)), sep = "")
- names(z$lambda) <- paste("Lam(",rep(namey,dat$nh), "_", rep(nameh, rep(dat$ny,dat$nh)), ")", sep = "")
+ namey <- colnames(P$dat$x[, 1:P$dat$ny])
+ names(z$coefficients) <- paste(rep(namey, P$dat$k), "_", rep(namex, rep(P$dat$ny, P$dat$k)), sep = "")
+ colnames(z$gt) <- paste(rep(namey, P$dat$nh), "_", rep(nameh, rep(P$dat$ny,P$dat$nh)), sep = "")
+ names(z$lambda) <- paste("Lam(",rep(namey,P$dat$nh), "_", rep(nameh, rep(P$dat$ny,P$dat$nh)), ")", sep = "")
}
- if (dat$ny == 1)
+ if (P$dat$ny == 1)
{
names(z$coefficients) <- namex
colnames(z$gt) <- nameh
names(z$lambda) <- nameh
}
- if (P$type == "EL")
- {
- z$badrho <- rhom$ch
- names(z$badrho) <- "Number_of_bad_rho"
- }
-
- G <- P$gradv(z$coefficients, x)
-
- khat <- crossprod(z$gt)/(n*P$k2)*P$bwVal
+ G <- P$gradv(z$coefficients, P$dat, z$pt)
+ khat <- crossprod(c(z$pt)*z$gt,z$gt)/(P$k2)*P$bwVal
G <- G/P$k1
kg <- solve(khat, G)
@@ -783,19 +782,18 @@
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))
b <- z$coefficients
- y <- as.matrix(model.response(dat$mf, "numeric"))
- ny <- dat$ny
+ y <- as.matrix(model.response(P$dat$mf, "numeric"))
+ ny <- P$dat$ny
b <- t(matrix(b, nrow = ny))
- x <- as.matrix(model.matrix(dat$mt, dat$mf, NULL))
+ x <- as.matrix(model.matrix(P$dat$mt, P$dat$mf, NULL))
yhat <- x%*%b
z$fitted.values <- yhat
z$residuals <- y - yhat
- z$terms <- dat$mt
- if(P$model) z$model <- dat$mf
+ z$terms <- P$dat$mt
+ if(P$model) z$model <- P$dat$mf
if(P$X) z$x <- x
if(P$Y) z$y <- y
z$call <- P$call
@@ -830,48 +828,37 @@
if(P$constraint)
res <- constrOptim(P$tet0, .thetf, grad = NULL, P = P, ...)
+ gt <- P$g(res$par, x)
+ rlamb <- getLamb(gt, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj,
+ method = P$optlam, k = P$k1/P$k2, control = P$controlLam)
- if (P$optlam == "iter")
- {
- rlamb <- getLamb(P$g, res$par, x, type = P$typel, tol_lam = P$tol_lam, maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, k = P$k1/P$k2)
- z <- list(coefficients = res$par, lambda = rlamb$lam, conv_lambda = rlamb$conv_mes, conv_par = res$convergence)
- z$foc_lambda <- rlamb$obj
- }
+ 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))
- if (P$optlam == "numeric")
- {
- gt <- P$g(res$par, x)
- rhofct <- function(lamb)
- {
- rhof <- -sum(rho(gt, lamb, type = P$typel, k = P$k1/P$k2)$rhomat)
- return(rhof)
- }
- rlamb <- optim(rep(0, ncol(gt)), rhofct, control = list(maxit = 1000))
- z <- list(coefficients = res$par, conv_par = res$convergence, lambda = rlamb$par)
- z$conv_lambda = paste("Lambda by optim. Conv. code = ", rlamb$convergence, sep = "")
- rho1 <- as.numeric(rho(gt, z$lambda, derive = 1, type = P$typel, k = P$k1/P$k2)$rhomat)
- z$foc_lambda <- crossprod(colMeans(rho1*gt))
- }
z$type <- P$type
- z$gt <- P$g(z$coefficients, x)
- 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)$rhomat/n
+ 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(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$rhomat) - rho(1, 0, type = P$typet, k = P$k1/P$k2)$rhomat)/n
+ z$objective <- sum(as.numeric(rhom) - .rho(1, 0, type = P$typet, k = P$k1/P$k2))/n
- if (P$type == "EL")
- {
- z$badrho <- rhom$ch
- names(z$badrho) <- "Number_of_bad_rho"
- }
-
if(P$gradvf)
G <- P$gradv(z$coefficients, x)
else
- G <- P$gradv(z$coefficients, x, g = P$g)
+ G <- P$gradv(z$coefficients, x, g = P$g, pt)
- khat <- crossprod(z$gt)/(n*P$k2)*P$bwVal
+ khat <- crossprod(c(z$pt)*z$gt, z$gt)/(P$k2)*P$bwVal
G <- G/P$k1
kg <- solve(khat, G)
Modified: pkg/gmm/man/gel.Rd
===================================================================
--- pkg/gmm/man/gel.Rd 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/man/gel.Rd 2011-11-30 23:00:47 UTC (rev 44)
@@ -12,7 +12,7 @@
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("iter","numeric"), model = TRUE,
+ optfct = c("optim", "optimize", "nlminb"), optlam = c("nlminb", "optim", "iter"), controlLam = list(), model = TRUE,
X = FALSE, Y = FALSE, TypeGel = "baseGel",...)
}
\arguments{
@@ -54,6 +54,8 @@
\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{controlLam}{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{TypeGel}{The name of the class object created by the method \code{getModel}. It allows developers to extand the package and create other GEL methods.}
Modified: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/man/getLamb.Rd 2011-11-30 23:00:47 UTC (rev 44)
@@ -8,15 +8,12 @@
It computes the vector of Lagrange multipliers, which maximizes the GEL objective function, using an iterative Newton method.
}
\usage{
-getLamb(g, tet, x, type = c('EL','ET','CUE'), tol_lam = 1e-12, maxiterlam = 1000, tol_obj = 1e-7, k = 1)
+getLamb(gt, type = c('EL','ET','CUE'), tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7, k = 1, method = c("nlminb", "optim", "iter"),
+ control = list())
}
\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.}
+\item{gt}{A \eqn{n \times q} matrix with typical element \eqn{g_i(\theta,x_t)}}
-\item{tet}{A \eqn{k \times 1} vector of parameters at which the function \eqn{g(\theta,x)} has to be evaluated}
-
-\item{x}{The matrix or vector of data from which the function \eqn{g(\theta,x)} is computed.}
-
\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting and "CUE" for continuous updated estimator.}
\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches \code{tol_lam} }
@@ -26,15 +23,18 @@
\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{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{control}{Controls to send to \code{\link{optim}}, \code{\link{nlminb}} or \code{\link{constrOptim}}}
}
\details{
-It solves the problem \eqn{\frac{1}{n}\sum_{t=1}^n D\rho(<g(\theta,x_t),\lambda>)g(\theta,x_t)=0}.
+It solves the problem \eqn{\max_{\lambda} - \frac{1}{n}\sum_{t=1}^n \rho(gt'\lambda)}.
}
\value{
lambda: A \eqn{q\times 1} vector of Lagrange multipliers which solve the system of equations given above.
-singular: 0 for a normal solution, 1 if the algorithm does not converge and 2 if the algorithm produces a singular system, NaN or Inf values.
-\code{conv_mes}: A message with details about the convergence.
+\code{conv}: Details on the type of convergence.
}
\references{
@@ -58,7 +58,8 @@
thet <- 0.2
sd <- .2
x <- matrix(arima.sim(n = n, list(order = c(2, 0, 1), ar = phi, ma = thet, sd = sd)), ncol = 1)
-getLamb(g, c(0, phi), x, type = "EL")
+gt <- g(c(0,phi),x)
+getLamb(gt, type = "EL",method="optim")
}
Deleted: pkg/gmm/man/rho.Rd
===================================================================
--- pkg/gmm/man/rho.Rd 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/man/rho.Rd 2011-11-30 23:00:47 UTC (rev 44)
@@ -1,48 +0,0 @@
-\name{rho}
-
-\alias{rho}
-
-\title{Objective function of Generalized Empirical Likelihood (GEL)}
-
-\description{
- It computes the objective function of GEL, its first and second analytical derivatives. It is called by \code{\link{gel}.}
-}
-\usage{
-rho(x, lamb, derive = 0, type = c("EL", "ET", "CUE"), drop = TRUE, k = 1)
-}
-\arguments{
-\item{x}{A \eqn{n\times q} matrix with typical element \eqn{(t,i)}, \eqn{g_i(\theta,x_t)}}
-
-\item{lamb}{A \eqn{q \times 1} vector of lagrange multipliers}
-
-\item{derive}{0 for the objective function, 1 for the first derivative with respect to \eqn{\lambda} and 2 for the second derivative with respect to \eqn{\lambda}.}
-
-\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting and "CUE" for continuous updated estimator.}
-
-\item{drop}{Because the solution may not be in the domain of \eqn{\rho(v)} \eqn{\forall t} in small sample, we can drop those observations to avoid the return of NaN}
-
-\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).}
-}
-\details{
-The objective function is \eqn{\frac{1}{n}\sum_{t=1}^n \rho(<g(\theta,x_t),\lambda>)}, where \eqn{\rho(v)=\log{(1-v)}} for empirical likelihood, \eqn{-e^v} for exponential tilting and \eqn{(-v-0.5v^2)} for continuous updated estimator.
-}
-
-\value{
-'rho' returns a scalar if "derive=0", a \eqn{q\time 1} vector if "derive" = 1 and a \eqn{q\times q} matrix if derive = 2.
- }
-
-\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.
-
-Hansen, L.P. and Heaton, J. and Yaron, A.(1996),
-Finit-Sample Properties of Some Alternative GMM Estimators.
-\emph{Journal of Business and Economic Statistics}, \bold{14}
-262-280.
-
-Smith, R.J. (2004), GEL Criteria for Moment Condition Models. \emph{Working paper, CEMMAP}.
-
-
-}
-
-
Modified: pkg/gmm/man/summary.Rd
===================================================================
--- pkg/gmm/man/summary.Rd 2011-11-30 15:59:56 UTC (rev 43)
+++ pkg/gmm/man/summary.Rd 2011-11-30 23:00:47 UTC (rev 44)
@@ -82,7 +82,7 @@
# GEL #
-t0 <- c(0,.5,.5)
+t0 <- res$coef
res <- gel(g, x, t0)
summary(res)
More information about the Gmm-commits
mailing list