[Gmm-commits] r107 - in pkg: gmm gmm/R gmm/man gmmExtra gmmExtra/R gmmExtra/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 13 23:26:24 CEST 2017
Author: chaussep
Date: 2017-06-13 23:26:23 +0200 (Tue, 13 Jun 2017)
New Revision: 107
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/man/gmm.Rd
pkg/gmmExtra/DESCRIPTION
pkg/gmmExtra/R/bootGmm.R
pkg/gmmExtra/man/bootGmm.Rd
pkg/gmmExtra/man/bootGmmMet.Rd
pkg/gmmExtra/man/bootJ.Rd
pkg/gmmExtra/man/linearHypothesis.Rd
Log:
Making gmmExtra work again
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmm/DESCRIPTION 2017-06-13 21:26:23 UTC (rev 107)
@@ -1,5 +1,5 @@
Package: gmm
-Version: 1.6-1
+Version: 1.6-2
Date: 2017-06-13
Title: Generalized Method of Moments and Generalized Empirical
Likelihood
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmm/R/getModel.R 2017-06-13 21:26:23 UTC (rev 107)
@@ -261,6 +261,7 @@
ar.method = object$ar.method,
approx = object$approx, tol = object$tol)
attr(object$x, "weight")$vcov <- object$vcov
+ attr(object$x, "mustar") <- object$mustar
object$g <- .momentFct
class(object) <- clname
return(object)
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmm/R/gmm.R 2017-06-13 21:26:23 UTC (rev 107)
@@ -15,7 +15,7 @@
kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen", "Tukey-Hanning"),crit=10e-7,bw = bwAndrews,
prewhite = 1, ar.method = "ols", approx="AR(1)",tol = 1e-7, itermax=100,optfct=c("optim","optimize","nlminb", "constrOptim"),
model=TRUE, X=FALSE, Y=FALSE, TypeGmm = "baseGmm", centeredVcov = TRUE, weightsMatrix = NULL, traceIter = FALSE, data, eqConst = NULL,
- eqConstFullVcov = FALSE, ...)
+ eqConstFullVcov = FALSE, mustar = NULL, ...)
{
type <- match.arg(type)
kernel <- match.arg(kernel)
@@ -37,7 +37,7 @@
crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method, approx = approx,
weightsMatrix = weightsMatrix, centeredVcov = centeredVcov, tol = tol, itermax = itermax,
optfct = optfct, model = model, X = X, Y = Y, call = match.call(), traceIter = traceIter,
- eqConst = eqConst, eqConstFullVcov = eqConstFullVcov)
+ eqConst = eqConst, eqConstFullVcov = eqConstFullVcov, mustar = mustar)
class(all_args)<-TypeGmm
Model_info<-getModel(all_args, ...)
z <- momentEstim(Model_info, ...)
@@ -52,7 +52,7 @@
"Tukey-Hanning"),crit=10e-7,bw = bwAndrews,
prewhite = FALSE, ar.method = "ols", approx="AR(1)",tol = 1e-7,
model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE,
- weightsMatrix = NULL, data)
+ weightsMatrix = NULL, data, mustar = NULL)
{
TypeGmm = "baseGmm"
type <- "eval"
@@ -74,7 +74,7 @@
centeredVcov = centeredVcov, tol = tol, itermax = 100,
model = model, X = X, Y = Y, call = match.call(),
traceIter = NULL, optfct="optim",
- eqConst = NULL, eqConstFullVcov = FALSE)
+ eqConst = NULL, eqConstFullVcov = FALSE, mustar = mustar)
class(all_args)<-TypeGmm
Model_info<-getModel(all_args)
class(Model_info) <- "baseGmm.eval"
@@ -261,6 +261,7 @@
}
includeExo <- which(colnames(xm)%in%colnames(hm))
inv <- attr(w, "inv")
+ mustar <- attr(dat, "mustar")
if (!is.null(type))
{
if(type=="2sls")
@@ -295,9 +296,7 @@
e2sls <- ym-xm%*%par
v2sls <- lm.fit(as.matrix(hm), e2sls)$fitted
value <- sum(v2sls^2)/sum(e2sls^2)
- }
- else
- {
+ } else {
par <- c(t(par))
g2sls <- g(par, dat)
w <- crossprod(g2sls)/n
@@ -305,44 +304,48 @@
value <- crossprod(gb, solve(w, gb))
}
}
- }
- else
- {
+ } else {
if (ny>1)
{
if (inv)
{
whx <- solve(w, (crossprod(hm, xm) %x% diag(ny)))
wvecyh <- solve(w, matrix(crossprod(ym, hm), ncol = 1))
- }
- else
- {
+ } else {
whx <- w%*% (crossprod(hm, xm) %x% diag(ny))
wvecyh <- w%*%matrix(crossprod(ym, hm), ncol = 1)
}
dg <- gradv(NULL, dat)
xx <- crossprod(dg, whx)
par <- solve(xx, crossprod(dg, wvecyh))
- }
- else
- {
+ } else {
if (nh>k)
{
if(inv)
xzwz <- crossprod(xm,hm)%*%solve(w,t(hm))
else
xzwz <- crossprod(xm,hm)%*%w%*%t(hm)
- par <- solve(xzwz%*%xm,xzwz%*%ym)
+ par <- solve(xzwz%*%xm,xzwz%*%ym)
+ } else {
+ par <- solve(crossprod(hm,xm),crossprod(hm,ym))
}
- else
- par <- solve(crossprod(hm,xm),crossprod(hm,ym)) }
+ }
gb <- matrix(colSums(g(par, dat))/n, ncol = 1)
if(inv)
value <- crossprod(gb, solve(w, gb))
else
value <- crossprod(gb, w%*%gb)
+ }
+ res <- list(par = par, value = value)
+ if (!is.null(mustar))
+ {
+ if (!is.null(type))
+ {
+ w <- crossprod(hm)/NROW(hm)
+ attr(w, "inv") <- TRUE
+ }
+ res <- .mustarLin(par, xm, hm, w, dat, mustar)
}
- res <- list(par = par, value = value)
if (!is.null(type))
{
if (type == "2sls")
@@ -357,6 +360,27 @@
return(res)
}
+.mustarLin <- function(par, xm, hm, w, dat, mustar)
+ {
+ if (ncol(xm) == ncol(hm))
+ {
+ par <- par-solve(crossprod(hm,xm),mustar)
+ } else {
+ hmxm <- crossprod(hm,xm)
+ if (attr(w, "inv"))
+ T1 <- solve(w, hmxm)
+ else
+ T1 <- w%*%hmxm
+ par <- par-solve(crossprod(hmxm, T1), crossprod(T1, mustar))
+ }
+ gb <- colSums(.momentFct(par, dat))/NROW(xm)
+ if(attr(w, "inv"))
+ value <- crossprod(gb, solve(w, gb))
+ else
+ value <- crossprod(gb, w%*%gb)
+ list(value=value, par=par)
+ }
+
.obj1 <- function(thet, x, w)
{
gt <- .momentFct(thet, x)
@@ -434,7 +458,14 @@
w <- attr(dat, "smooth")$w
gt <- smoothG(gt, bw = bw, weights = w)$smoothx
}
- return(as.matrix(gt))
+ gt <- as.matrix(gt)
+ if (!is.null(attr(dat, "mustar")))
+ {
+ if (length(attr(dat, "mustar")) != ncol(gt))
+ stop("dim of mustar must match the number of moment conditions")
+ gt <- sweep(gt, 2, attr(dat, "mustar"), "-")
+ }
+ return(gt)
}
.DmomentFct <- function(tet, dat, pt=NULL)
Modified: pkg/gmm/man/gmm.Rd
===================================================================
--- pkg/gmm/man/gmm.Rd 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmm/man/gmm.Rd 2017-06-13 21:26:23 UTC (rev 107)
@@ -17,12 +17,13 @@
tol = 1e-7, itermax=100,optfct=c("optim","optimize","nlminb", "constrOptim"),
model=TRUE, X=FALSE, Y=FALSE, TypeGmm = "baseGmm", centeredVcov = TRUE,
weightsMatrix = NULL, traceIter = FALSE, data, eqConst = NULL,
- eqConstFullVcov = FALSE, ...)
+ eqConstFullVcov = FALSE, mustar = NULL, ...)
evalGmm(g, x, t0, tetw=NULL, gradv=NULL, wmatrix = c("optimal","ident"),
vcov=c("HAC","iid","TrueFixed"), kernel=c("Quadratic Spectral","Truncated",
"Bartlett", "Parzen", "Tukey-Hanning"),crit=10e-7,bw = bwAndrews,
prewhite = FALSE, ar.method = "ols", approx="AR(1)",tol = 1e-7,
- model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, weightsMatrix = NULL, data)
+ model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, weightsMatrix = NULL,
+ data, mustar = NULL)
gmmWithConst(obj, which, value)
}
\arguments{
@@ -88,6 +89,13 @@
\item{eqConstFullVcov}{If FALSE, the constrained coefficients are assumed to be fixed and only the covariance of the unconstrained coefficients is computed. If TRUE, the covariance matrix of the full set of coefficients is computed.}
+\item{mustar}{If not null, it must be a vector with the number of
+ elements being equal to the number of moment conditions. In that case,
+ the vector is subtracted from the sample moment vector before
+ minimizing the objective function. It is useful to do a bootstrap
+ procedure.
+}
+
\item{...}{More options to give to \code{\link{optim}}.}
}
Modified: pkg/gmmExtra/DESCRIPTION
===================================================================
--- pkg/gmmExtra/DESCRIPTION 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmmExtra/DESCRIPTION 2017-06-13 21:26:23 UTC (rev 107)
@@ -1,6 +1,6 @@
Package: gmmExtra
-Version: 0.0-2
-Date: 2012-10-25
+Version: 0.0-3
+Date: 2017-06-13
Title: Extra tools for GMM estimation
Author: Pierre Chausse <pchausse at uwaterloo.ca>
Maintainer: Pierre Chausse <pchausse at uwaterloo.ca>
Modified: pkg/gmmExtra/R/bootGmm.R
===================================================================
--- pkg/gmmExtra/R/bootGmm.R 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmmExtra/R/bootGmm.R 2017-06-13 21:26:23 UTC (rev 107)
@@ -45,29 +45,31 @@
}
.getBlock <- function(x, l)
- {
- aX <- attributes(x)
- x <- as.matrix(x)
- n <- nrow(x)
- b <- floor(n/l)
- reDist <- n-b*l
- if (reDist != 0)
- {
- N <- sample(1:(n-l), reDist, replace=TRUE)
- i <- lapply(N, function(j) j:(j+l))
- i <- simplify2array(i)
- } else {
- i <- vector()
- }
- N <- sample(1:(n-l+1), (b-reDist), replace=TRUE)
- i2 <- lapply(N, function(j) j:(j+l-1))
- i2 <- simplify2array(i2)
- xF <- x[c(i,i2),]
- attributes(xF) <- aX
- attr(xF,"Blocks") <- c(i,i2)
- attr(xF,"l") <- l
- xF
- }
+ {
+ if (is.list(x))
+ n <- NROW(x[[1]])
+ aX <- attributes(x)
+ b <- floor(n/l)
+ reDist <- n - b * l
+ if (reDist != 0) {
+ N <- sample(1:(n - l), reDist, replace = TRUE)
+ i <- lapply(N, function(j) j:(j + l))
+ i <- simplify2array(i)
+ } else {
+ i <- vector()
+ }
+ N <- sample(1:(n - l + 1), (b - reDist), replace = TRUE)
+ i2 <- lapply(N, function(j) j:(j + l - 1))
+ i2 <- simplify2array(i2)
+ if (class(x) == "list")
+ xF <- lapply(1:length(x), function(j) as.matrix(x[[j]])[c(i,i2), ])
+ else
+ xF <- x[c(i, i2), ]
+ attributes(xF) <- aX
+ attr(xF, "Blocks") <- c(i, i2)
+ attr(xF, "l") <- l
+ xF
+ }
.getS <- function(gt, Blocks, l)
{
@@ -136,95 +138,70 @@
}
bootGmm <- function(obj, N, seed = NULL, niter=8, trace=TRUE)
- {
+ {
if (Sys.info()[["sysname"]] == "Windows")
- niter <- 8
- if(!is.null(seed))
- set.seed(seed)
- argCall <- obj$allArg
- argCall$call <- NULL
- argCall$t0 <- obj$coef
- bw <- attr(obj$w0, "Spec")$bw
- argCall$bw <- bw
- l <- attr(obj$w0, "Spec")$bw
- if (!is.null(l))
- {
- l <- floor(l)
- } else {
- l <- 1
- }
- l <- max(l,1)
+ {
+ warning("mc.cores set to 1 because of the Windows OS")
+ niter <- 1
+ }
+ if (!is.null(seed))
+ set.seed(seed)
+ if (is.null(obj$allArg$data))
+ stop("To run bootstrap GMM, gmm() must be called with non-null data argument")
+ l <- attr(obj$w0, "Spec")$bw
+ if (!is.null(l))
+ {
+ l <- floor(l)
+ } else {
+ l <- 1
+ }
+ l <- max(l, 1)
+ gtBlock <- na.omit(filter(obj$gt, rep(1/l, l)))
+ mustar <- colMeans(gtBlock)
+ Newdat <- list()
+ for (i in 1:N) Newdat[[i]] <- .getBlock(obj$allArg$data, l)
+ getAll <- function(i)
+ {
+ res0 <- update(obj, wmatrix="ident", bw=l, t0=obj$coef,data=Newdat[[i]],
+ mustar=mustar)
+ S <- .getS(res0$gt, attr(Newdat[[i]], "Blocks"),
+ attr(Newdat[[i]], "l"))
+ weightsMatrix <- try(solve(S), silent = TRUE)
+ if (class(weightsMatrix) == "try-error") {
+ stop("Singular S")
+ } else {
+ res1 <- update(obj, wmatrix="optimal", vcov="TrueFixed",
+ weightsMatrix=weightsMatrix,
+ mustar=mustar, t0=obj$coef, bw=l, data=Newdat[[i]])
+ res1$S <- S
+ }
+ res1
+ }
+ res <- .SimByBlock(getAll, niter, N, trace = trace)
+ chk <- sapply(1:N, function(i) !res[[i]]$bad)
+ res <- res[chk]
+ attr(res, "gmm") <- obj
+ class(res) <- "bootGmm"
+ return(res)
+ }
- gtBlock <- na.omit(filter(obj$gt,rep(1/l,l)))
- mustar <- colMeans(gtBlock)
-
- if(attr(obj$dat,"ModelType") == "linear")
- {
- argCall$x <- obj$dat
- attr(argCall$x, "oldg") <- obj$g
- dat <- obj$dat$x
- } else {
- dat <- argCall$x
- attr(argCall$x, "ModelType") <- "nonlinear"
- attr(argCall$x, "oldg") <- argCall$g
- }
- attr(argCall$x, "mu") <- mustar
-
- Newdat <- list()
- for (i in 1:N)
- Newdat[[i]] <- .getBlock(dat, l)
-
- getAll <- function(i, argCall, Newdat)
- {
- g2 <- function(theta, x)
- {
- mustar <- attr(x,"mu")
- gt <- attr(x,"oldg")(theta, x)
- sweep(gt, 2, mustar)
- }
- if (attr(argCall$x, "ModelType") == "nonlinear") {
- argCall$x <- Newdat[[i]]
- } else {
- argCall$x$x <- Newdat[[i]]}
-
- argCall$g <- g2
- argCall$wmatrix <- "ident"
-
- res0 <- do.call(gmm, argCall)
- argCall$wmatrix <- "optimal"
- argCall$vcov <- "TrueFixed"
- S <- .getS(res0$gt, attr(Newdat[[i]],"Blocks"), attr(Newdat[[i]],"l"))
- argCall$weightsMatrix <- try(solve(S),silent=TRUE)
- if (class(argCall$weightsMatrix) == "try-error")
- {
- stop("Singular S")
- } else {
- res1 <- do.call(gmm, argCall)
- res1$S <- S
- res1
- }
- }
- res <- .SimByBlock(getAll, niter, N, trace=trace, argCall=argCall, Newdat=Newdat)
- chk <- sapply(1:N, function(i) !res[[i]]$bad)
- res <- res[chk]
- attr(res,"gmm") <- obj
- class(res) <- "bootGmm"
- return(res)
- }
-
summary.bootGmm <- function(object, ...)
{
n <- length(object)
if (n == 0)
stop("The bootGmm object is empty")
coef <- sapply(1:n, function(i) object[[i]]$ans$coefficients)
- coef <- t(coef)
- conv <- sapply(1:n, function(i) object[[i]]$ans$algoInfo$convergence)
+ coef <- t(coef)
cat("Summary Statistics of Boostrap Estimates (N=",n,")\n")
cat("#######################################################\n")
- print(ans <- summary(coef))
- if (!is.null(conv[1]))
- cat("\nThe number of estimation that did not converge is ", sum(conv),"\n")
+ print(ans <- summary(coef))
+ if (attr(object[[1]]$ans$dat, "ModelType") != "linear")
+ {
+ conv <- sapply(1:n, function(i) object[[i]]$ans$algoInfo$convergence)
+ if (!is.null(conv[1]))
+ cat("\nThe number of estimation that did not converge is ",
+ sum(conv),"\n")
+ }
cat("The original estimates are\n")
print(round(attr(object,"gmm")$coefficients,4))
}
@@ -235,21 +212,29 @@
cat("#################\n")
cat("# GMM Bootstrap #\n")
cat("#################\n\n")
- cat("Original Call:\n ", paste(deparse(attr(x,"gmm")$call), sep = "\n", collapse = "\n"), "\n")
+ cat("Original Call:\n ", paste(deparse(attr(x,"gmm")$call),
+ sep = "\n", collapse = "\n"), "\n")
cat("Number of resamplings: ", n,"\n")
- conv <- sapply(1:n, function(i) x[[i]]$ans$algoInfo$convergence)
- cat("Number of successful estimations (convergence): ",n-sum(conv),"\n")
- coef <- t(sapply(1:n, function(i) x[[i]]$ans$coefficients))
- ansB <- rbind(colMeans(coef), apply(coef,2,sd))
- rownames(ansB) <- c("Mean","S-dev")
- cat("\nBootstrap results:\n")
- printCoefmat(ansB)
- res <- attr(x,"gmm")
- ans <- res$coefficients
- ans <- rbind(ans, summary(res)$coefficients[,2])
- cat("\nOriginal estimation:\n")
- rownames(ans) <- c("Estimates","S-dev")
- printCoefmat(ans)
+ if (n>0)
+ {
+ if (attr(x[[1]]$ans$dat, "ModelType") != "linear")
+ {
+ conv <- sapply(1:n, function(i) x[[i]]$ans$algoInfo$convergence)
+ cat("Number of successful estimations (convergence): ",
+ n-sum(conv),"\n")
+ }
+ coef <- t(sapply(1:n, function(i) x[[i]]$ans$coefficients))
+ ansB <- rbind(colMeans(coef), apply(coef,2,sd))
+ rownames(ansB) <- c("Mean","S-dev")
+ cat("\nBootstrap results:\n")
+ printCoefmat(ansB)
+ res <- attr(x,"gmm")
+ ans <- res$coefficients
+ ans <- rbind(ans, summary(res)$coefficients[,2])
+ cat("\nOriginal estimation:\n")
+ rownames(ans) <- c("Estimates","S-dev")
+ printCoefmat(ans)
+ }
}
Modified: pkg/gmmExtra/man/bootGmm.Rd
===================================================================
--- pkg/gmmExtra/man/bootGmm.Rd 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmmExtra/man/bootGmm.Rd 2017-06-13 21:26:23 UTC (rev 107)
@@ -49,7 +49,7 @@
set.seed(123)
d <- getdat(500)
-res <- gmm(d$Y~d$X-1,~d$Z-1,vcov="iid")
+res <- gmm(Y~X-1,~Z-1,vcov="iid", data=d)
# Just resampling 25 time to save time
resB <- bootGmm(res, 25, seed = 123, niter = 2, trace=TRUE)
}
Modified: pkg/gmmExtra/man/bootGmmMet.Rd
===================================================================
--- pkg/gmmExtra/man/bootGmmMet.Rd 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmmExtra/man/bootGmmMet.Rd 2017-06-13 21:26:23 UTC (rev 107)
@@ -53,7 +53,7 @@
set.seed(123)
d <- getdat(500)
-res <- gmm(d$Y~d$X-1,~d$Z-1)
+res <- gmm(Y~X-1,~Z-1, data=d)
# Just resampling 25 time to save time
resB <- bootGmm(res, 25, seed = 123, niter = 2, trace = TRUE)
resB
Modified: pkg/gmmExtra/man/bootJ.Rd
===================================================================
--- pkg/gmmExtra/man/bootJ.Rd 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmmExtra/man/bootJ.Rd 2017-06-13 21:26:23 UTC (rev 107)
@@ -44,7 +44,7 @@
set.seed(123)
d <- getdat(500)
-res <- gmm(d$Y~d$X-1,~d$Z-1)
+res <- gmm(Y~X-1,~Z-1, data=d)
resB <- bootGmm(res, 25, seed = 123, niter = 2, trace=TRUE)
J <- bootJ(resB)
J$test
Modified: pkg/gmmExtra/man/linearHypothesis.Rd
===================================================================
--- pkg/gmmExtra/man/linearHypothesis.Rd 2017-06-13 16:38:03 UTC (rev 106)
+++ pkg/gmmExtra/man/linearHypothesis.Rd 2017-06-13 21:26:23 UTC (rev 107)
@@ -47,9 +47,9 @@
set.seed(123)
d <- getdat(500)
-res <- gmm(d$Y~d$X-1,~d$Z-1)
+res <- gmm(Y~X-1,~Z-1, data=d)
# Just resampling 25 time to save time
resB <- bootGmm(res, 25, seed = 123, niter = 10, trace=TRUE)
-T <- linearHypothesis(resB, "d$X2 = 1")
+T <- linearHypothesis(resB, "X2 = 1")
T
}
More information about the Gmm-commits
mailing list