[Gmm-commits] r84 - in pkg/gmm: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 28 17:07:47 CET 2015
Author: chaussep
Date: 2015-10-28 17:07:46 +0100 (Wed, 28 Oct 2015)
New Revision: 84
Modified:
pkg/gmm/NEWS
pkg/gmm/R/FinRes.R
pkg/gmm/R/Methods.gel.R
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/momentEstim.R
pkg/gmm/R/specTest.R
Log:
Other change to meet whats in NEWS
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/NEWS 2015-10-28 16:07:46 UTC (rev 84)
@@ -14,6 +14,8 @@
o There is also an evalGel(). You fix the parameters and only the lambda are computed.
o Fixed NAMESPACE to avoid the notes given by CRAN.
o Modified getDat to allow for ~1 as instruments. Useful for inference on the mean using empirical likelihood.
+o The Lambdas and specTest from gel are printed even if the model is just identified. It is consitent with gmm which prints the j-test
+ for just identified models. The reason if to allow evalGel to generate the Lamdbas for testing purpose.
Changes in version 1.5-2
Modified: pkg/gmm/R/FinRes.R
===================================================================
--- pkg/gmm/R/FinRes.R 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/R/FinRes.R 2015-10-28 16:07:46 UTC (rev 84)
@@ -13,107 +13,103 @@
FinRes <- function(z, object, ...)
- {
+ {
# object is computed by the getModel method #
- UseMethod("FinRes")
- }
+ UseMethod("FinRes")
+ }
FinRes.baseGmm.res <- function(z, object, ...)
{
- P <- object
- x <- z$dat
- n <- ifelse(is.null(nrow(z$gt)),length(z$gt),nrow(z$gt))
-
- if (!is.null(attr(x,"eqConst")) & P$allArg$eqConstFullVcov)
- {
- eqConst <- attr(x,"eqConst")$eqConst
- coef <- rep(0,length(eqConst[,1])+length(z$coefficients))
- ncoef <- rep("",length(eqConst[,1])+length(z$coefficients))
- coef[-eqConst[,1]] <- z$coefficients
- ncoef[-eqConst[,1]] <- names(z$coefficients)
- coef[eqConst[,1]] <- eqConst[,2]
- ncoef[eqConst[,1]] <- rownames(eqConst)
- names(coef) <- ncoef
- z$coefficients <- coef
- if (!is.null(z$initTheta))
- {
- initTheta <- rep(0,length(z$coefficients))
- initTheta[eqConst[,1]] <- eqConst[,2]
- initTheta[-eqConst[,1]] <- z$initTheta
- z$initTheta <- initTheta
- }
- z$k <- z$k+nrow(eqConst)
- z$k2 <- z$k2+nrow(eqConst)
- attr(x, "eqConst") <- NULL
- z$specMod <- paste(z$specMod, "** Note: Covariance matrix computed for all coefficients based on restricted values **\n\n")
+ P <- object
+ x <- z$dat
+ n <- ifelse(is.null(nrow(z$gt)),length(z$gt),nrow(z$gt))
+
+ if (!is.null(attr(x,"eqConst")) & P$allArg$eqConstFullVcov)
+ {
+ eqConst <- attr(x,"eqConst")$eqConst
+ coef <- rep(0,length(eqConst[,1])+length(z$coefficients))
+ ncoef <- rep("",length(eqConst[,1])+length(z$coefficients))
+ coef[-eqConst[,1]] <- z$coefficients
+ ncoef[-eqConst[,1]] <- names(z$coefficients)
+ coef[eqConst[,1]] <- eqConst[,2]
+ ncoef[eqConst[,1]] <- rownames(eqConst)
+ names(coef) <- ncoef
+ z$coefficients <- coef
+ if (!is.null(z$initTheta))
+ {
+ initTheta <- rep(0,length(z$coefficients))
+ initTheta[eqConst[,1]] <- eqConst[,2]
+ initTheta[-eqConst[,1]] <- z$initTheta
+ z$initTheta <- initTheta
+ }
+ z$k <- z$k+nrow(eqConst)
+ z$k2 <- z$k2+nrow(eqConst)
+ attr(x, "eqConst") <- NULL
+ z$specMod <- paste(z$specMod, "** Note: Covariance matrix computed for all coefficients based on restricted values **\n\n")
+ }
+ z$G <- z$gradv(z$coefficients, x)
+ G <- z$G
+ v <- .weightFct(z$coefficient, x, P$vcov)
+ z$v <- v
+ if (P$vcov == "TrueFixed")
+ {
+ z$vcov=try(solve(crossprod(G, P$weightsMatrix) %*% G)/n, silent = TRUE)
+ if(class(z$vcov) == "try-error")
+ {
+ z$vcov <- matrix(Inf,length(z$coef),length(z$coef))
+ warning("The covariance matrix of the coefficients is singular")
+ }
+ } else if ( (is.null(P$weightsMatrix)) & (P$wmatrix != "ident") ) {
+ z$vcov <- try(solve(crossprod(G, solve(v, G)))/n, silent = TRUE)
+ if(class(z$vcov) == "try-error")
+ {
+ z$vcov <- matrix(Inf,length(z$coef),length(z$coef))
+ warning("The covariance matrix of the coefficients is singular")
+ }
+ } else {
+ if (is.null(P$weightsMatrix))
+ w <- diag(ncol(z$gt))
+ else
+ w <- P$weightsMatrix
+ if (dim(G)[1] == dim(G)[2]){
+ T1 <- try(solve(G), silent=TRUE)
+ } else {
+ T1 <- try(solve(t(G)%*%w%*%G,t(G)%*%w), silent = TRUE)
+ }
+ if(class(T1) == "try-error")
+ {
+ z$vcov <- matrix(Inf,length(z$coef),length(z$coef))
+ warning("The covariance matrix of the coefficients is singular")
+ } else {
+ z$vcov <- T1%*%v%*%t(T1)/n
+ }
+ }
+ dimnames(z$vcov) <- list(names(z$coefficients), names(z$coefficients))
+ z$call <- P$call
+
+ if(is.null(P$weightsMatrix))
+ {
+ if(P$wmatrix == "ident")
+ {
+ z$w <- diag(ncol(z$gt))
+ } else {
+ z$w <- try(solve(v), silent = TRUE)
+ if(class(z$w) == "try-error")
+ warning("The covariance matrix of the moment function is singular")
+ }
+ } else {
+ z$w <- P$weightsMatrix
+ }
+
+ z$weightsMatrix <- P$weightsMatrix
+ z$infVcov <- P$vcov
+ z$infWmatrix <- P$wmatrix
+ z$allArg <- P$allArg
+ z$met <- P$type
+ z$kernel <- P$kernel
+ z$coefficients <- c(z$coefficients)
+ class(z) <- "gmm"
+ return(z)
}
- z$G <- z$gradv(z$coefficients, x)
- G <- z$G
- v <- .weightFct(z$coefficient, x, P$vcov)
- z$v <- v
- if (P$vcov == "TrueFixed")
- {
- z$vcov=try(solve(crossprod(G, P$weightsMatrix) %*% G)/n, silent = TRUE)
- if(class(z$vcov) == "try-error")
- {
- z$vcov <- matrix(Inf,length(z$coef),length(z$coef))
- warning("The covariance matrix of the coefficients is singular")
- }
- }
- else if ( (is.null(P$weightsMatrix)) & (P$wmatrix != "ident") )
- {
- z$vcov <- try(solve(crossprod(G, solve(v, G)))/n, silent = TRUE)
- if(class(z$vcov) == "try-error")
- {
- z$vcov <- matrix(Inf,length(z$coef),length(z$coef))
- warning("The covariance matrix of the coefficients is singular")
- }
- }
- else
- {
- if (is.null(P$weightsMatrix))
- w <- diag(ncol(z$gt))
- else
- w <- P$weightsMatrix
- if (dim(G)[1] == dim(G)[2]){
- T1 <- try(solve(G), silent=TRUE)
- } else {
- T1 <- try(solve(t(G)%*%w%*%G,t(G)%*%w), silent = TRUE)
- }
- if(class(T1) == "try-error")
- {
- z$vcov <- matrix(Inf,length(z$coef),length(z$coef))
- warning("The covariance matrix of the coefficients is singular")
- }
- else
- z$vcov <- T1%*%v%*%t(T1)/n
- }
- dimnames(z$vcov) <- list(names(z$coefficients), names(z$coefficients))
- z$call <- P$call
-
- if(is.null(P$weightsMatrix))
- {
- if(P$wmatrix == "ident")
- z$w <- diag(ncol(z$gt))
- else
- {
- z$w <- try(solve(v), silent = TRUE)
- if(class(z$w) == "try-error")
- warning("The covariance matrix of the moment function is singular")
- }
- }
- else
- z$w <- P$weightsMatrix
- z$weightsMatrix <- P$weightsMatrix
- z$infVcov <- P$vcov
- z$infWmatrix <- P$wmatrix
- z$allArg <- P$allArg
- z$met <- P$type
- z$kernel <- P$kernel
- z$coefficients <- c(z$coefficients)
- class(z) <- "gmm"
- return(z)
- }
-
Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/R/Methods.gel.R 2015-10-28 16:07:46 UTC (rev 84)
@@ -75,14 +75,11 @@
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits),
print.gap = 2, quote = FALSE)
- if (length(x$coefficients)<length(x$lambda))
- {
- cat("\n")
- cat("Lambdas:\n")
- print.default(format(coef(x, lambda = TRUE), digits = digits),
- print.gap = 2, quote = FALSE)
- }
cat("\n")
+ 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")
if (length(x$coefficients)<length(x$lambda))
cat("Convergence code for Lambda: ", x$conv_lambda$convergence,"\n")
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/R/getModel.R 2015-10-28 16:07:46 UTC (rev 84)
@@ -68,10 +68,7 @@
{
object$gradv <- .DmomentFct
object$gradvf <- FALSE
- if (is.null(object$data))
- dat <- getDat(object$g, object$x)
- else
- dat <- getDat(object$g, object$x, object$data)
+ dat <- getDat(object$g, object$x, data = object$data)
if(is.null(object$weightsMatrix))
{
clname <- paste(class(object), ".", object$type, ".formula", sep = "")
@@ -81,11 +78,11 @@
}
object$x <- dat
object$gform<-object$g
- 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(dat$x[,(dat$ny+1):(dat$ny+dat$k), drop=FALSE])
+ nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh), drop=FALSE])
if (dat$ny > 1)
{
- namey <- colnames(dat$x[,1:dat$ny])
+ namey <- colnames(dat$x[,1:dat$ny, drop=FALSE])
object$namesCoef <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
object$namesgt <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny, dat$nh)), sep = "")
} else {
@@ -210,8 +207,8 @@
object$gradvf <- FALSE
object$x <- dat
object$gform<-object$g
- 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(dat$x[,(dat$ny+1):(dat$ny+dat$k), drop=FALSE])
+ nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh), drop=FALSE])
if (dat$ny > 1)
{
namey <- colnames(dat$x[,1:dat$ny])
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/R/gmm.R 2015-10-28 16:07:46 UTC (rev 84)
@@ -80,38 +80,39 @@
}
tsls <- function(g,x,data)
-{
-if(class(g) != "formula")
- stop("2SLS is for linear models expressed as formula only")
-ans <- gmm(g,x,data=data,vcov="iid")
-ans$met <- "Two Stage Least Squares"
-ans$call <- match.call()
-class(ans) <- c("tsls","gmm")
-return(ans)
-}
+ {
+ if(class(g) != "formula")
+ stop("2SLS is for linear models expressed as formula only")
+ ans <- gmm(g,x,data=data,vcov="iid")
+ ans$met <- "Two Stage Least Squares"
+ ans$call <- match.call()
+ class(ans) <- c("tsls","gmm")
+ return(ans)
+ }
.myKernHAC <- function(gmat, obj)
- {
+ {
+ gmat <- as.matrix(gmat)
if(obj$centeredVcov)
- gmat <- lm(gmat~1)
+ gmat <- lm(gmat~1)
else
- class(gmat) <- "gmmFct"
+ class(gmat) <- "gmmFct"
AllArg <- obj$WSpec$sandwich
AllArg$x <- gmat
if (is.function(AllArg$bw))
- {
+ {
bw <- AllArg$bw(gmat, order.by = AllArg$order.by, kernel = AllArg$kernel,
- prewhite = AllArg$prewhite, ar.method = AllArg$ar.method)
+ prewhite = AllArg$prewhite, ar.method = AllArg$ar.method)
AllArg$bw <- bw
- }
+ }
weights <- do.call(weightsAndrews,AllArg)
AllArg$sandwich <- FALSE
AllArg$weights <- weights
w <- do.call(vcovHAC, AllArg)
attr(w,"Spec") <- list(weights = weights, bw = AllArg$bw, kernel = AllArg$kernel)
w
- }
+ }
getDat <- function (formula, h, data)
{
@@ -127,13 +128,15 @@
y <- as.matrix(model.response(mf, "numeric"))
xt <- as.matrix(model.matrix(mt, mf, NULL))
n <- NROW(y)
-
if (inherits(h,'formula'))
- {
+ {
+ tmp <- as.character(formula)
+ h <- paste(tmp[2], "~", as.character(h)[2], sep="")
+ h <- as.formula(h)
mfh <- match.call(expand.dots = FALSE)
mh <- match(c("h", "data"), names(mfh), 0L)
mfh <- mfh[c(1L, mh)]
- mfh$formula <- mfh$h
+ mfh$formula <- h
mfh$h <- NULL
mfh$drop.unused.levels <- TRUE
mfh$na.action <- "na.pass"
@@ -141,19 +144,11 @@
mfh <- eval(mfh, parent.frame())
mth <- attr(mfh, "terms")
h <- as.matrix(model.matrix(mth, mfh, NULL))
- if (length(h) == 0)
- {
- h <- as.matrix(rep(1,n))
- colnames(h) <- "h.(Intercept)"
- }
}
else
- {
- if (!is.matrix(h))
- h <- cbind(rep(1,length(h)),h)
- else
- h <- cbind(rep(1,nrow(h)),h)
- h <- as.matrix(h)
+ {
+ h <- as.matrix(h)
+ h <- cbind(1,h)
if(is.null(colnames(h)))
colnames(h) <- c("h.(Intercept)",paste("h",1:(ncol(h)-1),sep=""))
else
@@ -166,9 +161,6 @@
ny <- ncol(y)
k <- ncol(xt)
nh <- ncol(h)
-
- if (nrow(y) != nrow(xt) | nrow(xt) != nrow(h) | nrow(y)!=nrow(h))
- stop("The number of observations of X, Y and H must be the same")
if (nh<k)
stop("The number of moment conditions must be at least equal to the number of coefficients to estimate")
if (is.null(colnames(y)))
@@ -378,7 +370,7 @@
w <- attr(dat, "smooth")$w
gt <- smoothG(gt, bw = bw, weights = w)$smoothx
}
- return(gt)
+ return(as.matrix(gt))
}
.DmomentFct <- function(tet, dat, pt=NULL)
@@ -409,7 +401,7 @@
if (!is.null(attr(dat, "eqConst")))
dgb <- dgb[,-attr(dat,"eqConst")$eqConst[,1], drop=FALSE]
}
- return(dgb)
+ return(as.matrix(dgb))
}
.weightFct <- function(tet, dat, type=c("HAC", "iid", "ident", "fct", "fixed"))
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/R/momentEstim.R 2015-10-28 16:07:46 UTC (rev 84)
@@ -674,14 +674,7 @@
G <- G/P$k1
kg <- solve(khat, G)
z$vcov_par <- solve(crossprod(G, kg))/n
- if (length(z$lambda) == length(z$coefficients))
- {
- z$vcov_lambda <- matrix(NA, rep(length(z$lambda), 2))
- z$lambda <- rep(NA, length(z$lambda))
- z$specMod <- paste(z$specMod, "\n Just identified model; no lambda nor specification test needed\n", sep="")
- } else {
- z$vcov_lambda <- solve(khat, ( diag(ncol(khat)) - G %*% (z$vcov_par*n) %*% t(kg) ))/n*P$bwVal^2
- }
+ z$vcov_lambda <- solve(khat, ( diag(ncol(khat)) - G %*% (z$vcov_par*n) %*% t(kg) ))/n*P$bwVal^2
z$weights <- P$w
z$bwVal <- P$bwVal
Modified: pkg/gmm/R/specTest.R
===================================================================
--- pkg/gmm/R/specTest.R 2015-10-27 21:37:12 UTC (rev 83)
+++ pkg/gmm/R/specTest.R 2015-10-28 16:07:46 UTC (rev 84)
@@ -49,6 +49,8 @@
df <- (ncol(x$gt) - length(x$coef))
ntest <- noquote(paste("Over-identifying restrictions tests: degrees of freedom is ", df, sep = ""))
vptest <- pchisq(test,df,lower.tail = FALSE)
+ if (df == 0)
+ vptest[] <- "###"
test <- cbind(test,vptest)
dimnames(test) <- list(c("LR test", "LM test", "J test"), c("statistics", "p-value"))
ans <- list(test = test, ntest = ntest)
More information about the Gmm-commits
mailing list