[Gmm-commits] r45 - in pkg/gmm: . R data inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 1 23:29:17 CET 2011
Author: chaussep
Date: 2011-12-01 23:29:16 +0100 (Thu, 01 Dec 2011)
New Revision: 45
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NEWS
pkg/gmm/R/Methods.gel.R
pkg/gmm/R/Methods.gmm.R
pkg/gmm/R/gel.R
pkg/gmm/R/getModel.R
pkg/gmm/R/momentEstim.R
pkg/gmm/data/Finance.rda
pkg/gmm/inst/doc/gmm_with_R.rnw
pkg/gmm/man/gel.Rd
pkg/gmm/man/getLamb.Rd
Log:
Many changes, still not stable
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/DESCRIPTION 2011-12-01 22:29:16 UTC (rev 45)
@@ -12,7 +12,7 @@
methods that belong to the Generalized Empirical Likelihood
(GEL) family of estimators, as presented by Smith(1997),
Kitamura(1997), Newey-Smith(2004) and Anatolyev(2005).
-Depends: R (>= 2.0.0), sandwich
+Depends: R (>= 2.10.0), sandwich
Suggests: mvtnorm, car, fBasics, MASS, timeDate, timeSeries
Imports: stats
License: GPL (>= 2)
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/NEWS 2011-12-01 22:29:16 UTC (rev 45)
@@ -11,7 +11,20 @@
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.
+o For GEL type CUE the implied probabilities are computed according to Antoine, Bonnal, and Renault (2007) which solves the
+ problem of negative probabilities.
+o The data file Finance has been resaved which implies that the package depends now on R version 2.10.0 and higher
+o The function rho() had been changed to .rho() because it is not useful outside the gel() function
+o The function getLamb() has been modified. It is now more efficient. The default for the gel() option optlam is now nlminb.
+ The gradient and the hessian is provided which makes it much faster to solve for lambda. The argument of the function also changed.
+ Instead of providing the function g() and the vector of coefficients, we provide the matrix gt. I may be useful to call the function
+ sometimes to solve for lambda and it is easier that way.
+o The choice "iter" from the option optlam is kept for cases in which the optimizer fails to solve for lambda. It is often the case
+ when type="ETEL" is chosen which tends to produce NA's.
+o For GEL of type "EL" and optfct="optim", the algorith constrOptim() is launched to make sure lambda'gt is always greater than 1
+ and avoid NA's when computing log(1-lambda'gt).
+o Sometimes, problems happen in GMM estimation because of the bad first step estimates used to compute the weighting matrix.
+ The first step estimates are usually computed using the identity matrix. The vector is now printed for better control.
Changes in version 1.3-8
Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/R/Methods.gel.R 2011-12-01 22:29:16 UTC (rev 45)
@@ -133,9 +133,7 @@
ans$conv_pt <- z$conv_pt
ans$conv_moment <- cbind(z$conv_moment)
ans$conv_lambda <- z$conv_lambda
- names(ans$conv_par) <- "Convergence_code_theta"
names(ans$conv_pt) <- "Sum_of_pt"
- names(ans$conv_lambda) <- "Convergence_code_for_lambda"
dimnames(ans$conv_moment) <- list(names(z$gt), "Sample_moment_with_pt")
class(ans) <- "summary.gel"
ans
Modified: pkg/gmm/R/Methods.gmm.R
===================================================================
--- pkg/gmm/R/Methods.gmm.R 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/R/Methods.gmm.R 2011-12-01 22:29:16 UTC (rev 45)
@@ -30,7 +30,10 @@
if(z$met=="cue")
ans$cue <- object$cue
if (!is.null(object$initTheta))
+ {
ans$initTheta <- object$initTheta
+ names(ans$initTheta) <- names(z$coefficients)
+ }
class(ans) <- "summary.gmm"
ans
}
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/R/gel.R 2011-12-01 22:29:16 UTC (rev 45)
@@ -20,7 +20,11 @@
if (derive == 0)
{
if (type == "EL")
+ {
+ if (any(gml>=1))
+ stop("Computation of Lambda fails because NAs produced by log(1-gt*l)")
rhomat <- log(1 - gml)
+ }
if (type == "ET")
rhomat <- -exp(gml)
if (type == "CUE")
@@ -50,8 +54,8 @@
}
-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())
+getLamb <- function(gt, type = c('EL', 'ET', 'CUE', "ETEL"), tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7,
+ k = 1, method = c("nlminb", "optim", "iter"), control=list())
{
method <- match.arg(method)
if(is.null(dim(gt)))
@@ -68,17 +72,17 @@
J <- crossprod(r2*gt,gt)
if (sum(abs(F))<tol_obj)
{
- conv <- "Tolerance for the FOC reached"
+ conv <- list(convergence="Tolerance for the FOC reached")
break
}
P <- solve(J,F)
if (sum(abs(P))<tol_lam)
{
- cov <- "Tolerance on lambda reached"
+ conv <- list(convergence="Tolerance on lambda reached")
break
}
l0 <- l0 + P
- cov <- "maxiterlam reached"
+ conv <- list(convergence="maxiterlam reached")
}
}
else
@@ -98,18 +102,27 @@
r2 <- .rho(X,l,derive=2,type=type,k=k)
-t(X*r2)%*%X/nrow(X)
}
- if (method=="optim")
+ if (type == "ETEL")
{
- if (type != "EL")
- res <- optim(rep(0,ncol(gt)),fct,gr=Dfct,X=gt,method="B",control=control)
+ 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 (type != "EL")
+ res <- optim(rep(0,ncol(gt)),fct,gr=Dfct,X=gt,method="BFGS",control=control)
+ else
+ {
+ ci <- -rep(1,nrow(gt))
+ res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt)
+ }
+ }
else
- {
- ci <- -rep(1,nrow(gt))
- res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt)
- }
+ 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, control = control)
l0 <- res$par
if (method == "optim" | method == "constrOptim")
@@ -166,7 +179,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("nlminb", "optim", "iter"), controlLam = list(), model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", ...)
+ optlam = c("nlminb", "optim", "iter"), Lambdacontrol = list(), model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", ...)
{
type <- match.arg(type)
@@ -181,7 +194,7 @@
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(),
- controlLam = controlLam)
+ Lambdacontrol = Lambdacontrol)
class(all_args)<-TypeGel
Model_info<-getModel(all_args)
@@ -192,24 +205,28 @@
}
- .thetf <- function(tet, P)
+ .thetf <- function(tet, P, output=c("obj","all"))
{
+ output <- match.arg(output)
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)
+ lamb <- try(getLamb(gt, 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), 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
+ lamb <- getLamb(gt, 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 = "optim")
}
else
- 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
+ lamb <- getLamb(gt, 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)
- obj <- mean(.rho(gt, lamb, type = P$typet, derive = 0, k = P$k1/P$k2))
-
- return(obj)
+ obj <- mean(.rho(gt, lamb$lambda, type = P$typet, derive = 0, k = P$k1/P$k2))
+ if(output == "obj")
+ return(obj)
+ else
+ return(list(obj = obj, lambda = lamb, gt = gt))
}
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/R/getModel.R 2011-12-01 22:29:16 UTC (rev 45)
@@ -100,6 +100,7 @@
getModel.baseGel <- function(object, ...)
{
+
P <- object
if (P$type == "ETEL")
{
@@ -190,9 +191,9 @@
P$k2 <- 2/3
}
P$g1 <- P$g
+
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,
+ P$bwVal <- P$bw(centeredGt <- lm(P$g(rgmm$coefficients, P$dat)~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/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/R/momentEstim.R 2011-12-01 22:29:16 UTC (rev 45)
@@ -713,22 +713,10 @@
if(P$constraint)
res <- constrOptim(P$tet0, .thetf, grad = NULL, P = P, ...)
- gt <- P$g(res$par, P$dat)
-
+ All <- .thetf(res$par, P, "all")
+ gt <- All$gt
+ rlamb <- All$lambda
- 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)
-
-
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))
@@ -737,9 +725,9 @@
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
+print("hello")
-
# Making sure pt>0
if (P$type=="CUE")
{
@@ -828,9 +816,9 @@
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)
+ All <- .thetf(res$par, P, "all")
+ gt <- All$gt
+ 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)
@@ -848,15 +836,14 @@
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
- G <- P$gradv(z$coefficients, x, g = P$g, pt)
+ G <- P$gradv(z$coefficients, x, g = P$g, z$pt)
+
khat <- crossprod(c(z$pt)*z$gt, z$gt)/(P$k2)*P$bwVal
G <- G/P$k1
Modified: pkg/gmm/data/Finance.rda
===================================================================
(Binary files differ)
Modified: pkg/gmm/inst/doc/gmm_with_R.rnw
===================================================================
--- pkg/gmm/inst/doc/gmm_with_R.rnw 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/inst/doc/gmm_with_R.rnw 2011-12-01 22:29:16 UTC (rev 45)
@@ -723,9 +723,11 @@
@
A fourth method is available which is called the exponentially tilted empirical likelihood (ETEL) and was proposed by \cite{schennach07}. However, it does not belong to the family of GEL estimators. It solves the problem of misspecified models. In such models there may not exist any pseudo value to which $\hat{\theta}$ converges as the sample size increases. ETEL uses the $\rho()$ of ET to solve for $\lambda$ and the $\rho()$ of EL to solve for $\theta$. \cite{schennach07} shows that ETEL shares the same asymptotic properties of EL without having to impose restrictions on the domain of $\rho(v)$ when solving for $\lambda$.
<<>>=
-res <- gel(g1,x1,tet0,type="ETEL")
-coef(res)
+res_etel <- gel(g1,x1,c(mu=1,sig=1),type="ETEL")
+coef(res_etel)
@
+The type ETEL is experimental for now. Although it is supposed to be more stable because no restrictions are required to solve for $\lambda$, once we substitute $\lambda(\theta)$ in the EL objective function to estimate $\theta$, we still need to restrict $\lambda'g_t$ to avoid having NA's. The solution used in gel() is to obtain $\lambda(\theta)$ with \code{constrOptim} with the restriction $\lambda'gt > 1$ even if it is not required by ET ($\rho(v)=-\exp{(v)}$). It is however sensitive to starting values. That's the reason why we used different ones above.
+
\subsection{Estimating the AR coefficients of an ARMA process}
Because the moment conditions are weakly dependent, we need to set the option \code{smooth=TRUE}. Before going to the estimation procedure, we need to understand the relationship between the smoothing kernel and the HAC estimator. The reason why we need to smooth the moment function is that GEL estimates the covariance matrix of $\bar{g}(\theta,x_t)$, as if we had iid observations, using the expression $(1/T)\sum_{t=1}^T(g_tg_t')$. We can show that substituting $g_t$ by $g^w_t$ in this expression results in an HAC estimator. However, the relationship between the smoothing kernel and the kernel that appears in the HAC estimator is not obvious. For example, we can show that if the smoothing kernel is Truncated, then the kernel in the HAC estimator is the Bartlett. Let us consider the truncated kernel with a bandwidth of 2. This implies that $w(s)=1/5$ for $|s|\leq 2$ and 0 otherwise. Then, the expression for the covariance matrix becomes:
Modified: pkg/gmm/man/gel.Rd
===================================================================
--- pkg/gmm/man/gel.Rd 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/man/gel.Rd 2011-12-01 22:29:16 UTC (rev 45)
@@ -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("nlminb", "optim", "iter"), controlLam = list(), model = TRUE,
+ optfct = c("optim", "optimize", "nlminb"), optlam = c("nlminb", "optim", "iter"), Lambdacontrol = list(), model = TRUE,
X = FALSE, Y = FALSE, TypeGel = "baseGel",...)
}
\arguments{
@@ -54,7 +54,7 @@
\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{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.}
Modified: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd 2011-11-30 23:00:47 UTC (rev 44)
+++ pkg/gmm/man/getLamb.Rd 2011-12-01 22:29:16 UTC (rev 45)
@@ -8,13 +8,13 @@
It computes the vector of Lagrange multipliers, which maximizes the GEL objective function, using an iterative Newton method.
}
\usage{
-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())
+getLamb(gt, type = c("EL","ET","CUE", "ETEL"), tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7, k = 1,
+ method = c("nlminb", "optim", "iter"), control = list())
}
\arguments{
\item{gt}{A \eqn{n \times q} matrix with typical element \eqn{g_i(\theta,x_t)}}
-\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting and "CUE" for continuous updated estimator.}
+\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting and "CUE" for continuous updated estimator. See details for "ETEL".}
\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches \code{tol_lam} }
@@ -29,7 +29,7 @@
\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)}.
+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{\ling{constrOptim}} is used with the restriction \eqn{\lambda'g_t < 1}.
}
\value{
More information about the Gmm-commits
mailing list