[Gmm-commits] r86 - in pkg/gmm: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 29 19:48:29 CET 2015
Author: chaussep
Date: 2015-10-29 19:48:28 +0100 (Thu, 29 Oct 2015)
New Revision: 86
Modified:
pkg/gmm/NAMESPACE
pkg/gmm/NEWS
pkg/gmm/R/Methods.gel.R
pkg/gmm/R/Methods.gmm.R
pkg/gmm/R/getModel.R
pkg/gmm/R/momentEstim.R
pkg/gmm/R/specTest.R
pkg/gmm/man/confint.Rd
Log:
See NEWS
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2015-10-28 19:28:30 UTC (rev 85)
+++ pkg/gmm/NAMESPACE 2015-10-29 18:48:28 UTC (rev 86)
@@ -12,7 +12,7 @@
momentEstim.baseGmm.cue, getModel.baseGmm, getModel.baseGel, getModel.constGmm, getModel.constGel,
FinRes.baseGmm.res, momentEstim.baseGel.mod, momentEstim.baseGel.modFormula,tsls,summary.tsls, print.summary.tsls,
KTest, print.gmmTests, gmmWithConst, estfun.tsls, model.matrix.tsls,vcov.tsls, bread.tsls, evalGmm, momentEstim.baseGmm.eval,
- momentEstim.baseGel.eval, evalGel, confint.gmm)
+ momentEstim.baseGel.eval, evalGel, confint.gmm, print.confint)
S3method(summary, gmm)
S3method(summary, tsls)
@@ -41,6 +41,7 @@
S3method(specTest, gmm)
S3method(specTest, gel)
S3method(print, specTest)
+S3method(print, confint)
S3method(FinRes, baseGmm.res)
S3method(getModel, baseGmm)
S3method(getModel, baseGel)
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2015-10-28 19:28:30 UTC (rev 85)
+++ pkg/gmm/NEWS 2015-10-29 18:48:28 UTC (rev 86)
@@ -13,9 +13,12 @@
for the HAC estimator are generated. It is possible to give a different vector os parameters for the moment function and the weighting matrix.
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 Modified getDat to allow for ~1 as instruments. Useful for inference on the mean using empirical likelihood with gel(x~1,~1).
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.
+ for just identified models. The reason if to allow evalGel to generate the Lamdbas for testing purpose.
+o confint.gxx() has been modified. It now create an object and a print method has been created as well.
+o confint.gel includes the possibility of building a confidence interval using the inversion of one of the three specification tests.
+ It is therefore possible de construct an Empirical Likelihood confidence interval as described in Owen 2001. See manual.
Changes in version 1.5-2
Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R 2015-10-28 19:28:30 UTC (rev 85)
+++ pkg/gmm/R/Methods.gel.R 2015-10-29 18:48:28 UTC (rev 86)
@@ -12,43 +12,133 @@
# http://www.r-project.org/Licenses/
-confint.gel <- function(object, parm, level = 0.95, lambda = FALSE, ...)
+.invTest <- function(object, which, level = 0.95, fact = 3, type=c("LR", "LM", "J"), corr=NULL)
{
+ type <- match.arg(type)
+ argCall <- object$allArg
+ if (length(which) > 1)
+ stop("tests are inverted only for one parameter")
+ if (is.character(which))
+ {
+ which <- which(which == names(object$coefficients))
+ if (length(which) == 0)
+ stop("Wrong parameter names")
+ }
+ wtest <- which(type==c("LR", "LM", "J"))
+ if (object$df > 0)
+ {
+ test0 <- specTest(object)$test[wtest,]
+ if (test0[2] < (1-level))
+ stop("The model is rejected at the point estimate")
+ test0 <- test0[1]
+ } else {
+ test0 <- 0
+ }
+ if (length(object$coefficients) == 1)
+ {
+ argCall$optfct <- NULL
+ argCall$constraint <- NULL
+ argCall$eqConst <- NULL
+ argCall$eqConstFullVcov <- NULL
+ fctCall <- "evalGel"
+ } else {
+ fctCall <- "gel"
+ argCall$eqConst <- which
+ }
+ if (length(object$coefficients) == 2)
+ {
+ sdcoef <- sqrt(diag(vcov(object))[-which])
+ coef <- object$coefficients[-which]
+ upper <- coef + fact*sdcoef
+ lower <- coef - fact*sdcoef
+ argCall$method <- "Brent"
+ argCall$lower <- lower
+ argCall$upper <- upper
+ }
+ sdcoef <- sqrt(diag(vcov(object))[which])
+ coef <- object$coefficients[which]
+ int1 <- c(coef, coef + fact*sdcoef)
+ int2 <- c(coef - fact*sdcoef, coef)
+ fct <- function(coef, which, wtest, level, test0, corr=NULL)
+ {
+ argCall$tet0 <- object$coefficients
+ argCall$tet0[which] <- coef
+ obj <- do.call(get(fctCall), argCall)
+ test <- as.numeric(specTest(obj)$test[wtest,1]) - test0
+ if (is.null(corr))
+ level - pchisq(test, 1)
+ else
+ level - pchisq(test/corr, 1)
+ }
+ res1 <- uniroot(fct, int1, which = which, wtest=wtest, level=level,
+ test0=test0, corr=corr)
+ res2 <- uniroot(fct, int2, which = which, wtest=wtest, level=level,
+ test0=test0, corr=corr)
+ sort(c(res1$root, res2$root))
+ }
+
+
+confint.gel <- function(object, parm, level = 0.95, lambda = FALSE,
+ type = c("Wald", "invLR", "invLM", "invJ"),
+ fact = 3, corr = NULL, ...)
+ {
+ type <- match.arg(type)
z <- object
n <- nrow(z$gt)
-
- se_par <- sqrt(diag(z$vcov_par))
- par <- z$coefficients
- tval <- par/se_par
-
- se_parl <- sqrt(diag(z$vcov_lambda))
- lamb <- z$lambda
-
- zs <- qnorm((1 - level)/2, lower.tail=FALSE)
- ch <- zs*se_par
-
- if(!lambda)
+ if (type == "Wald")
{
- ans <- cbind(par-ch, par+ch)
- dimnames(ans) <- list(names(par), c((1 - level)/2, 0.5+level/2))
- }
- if(lambda)
- {
- if (length(z$coefficients) == length(z$lambda))
+ ntest <- "Direct Wald type confidence interval"
+ se_par <- sqrt(diag(z$vcov_par))
+ par <- z$coefficients
+ tval <- par/se_par
+ se_parl <- sqrt(diag(z$vcov_lambda))
+ lamb <- z$lambda
+ zs <- qnorm((1 - level)/2, lower.tail=FALSE)
+ ch <- zs*se_par
+ if(!lambda)
{
- cat("\nNo confidence intervals for lambda when the model is just identified.\n")
- return(NULL)
- } else {
- chl <- zs*se_parl
- ans <- cbind(lamb - chl, lamb + chl)
- dimnames(ans) <- list(names(lamb), c((1 - level)/2, 0.5 + level/2))
+ ans <- cbind(par-ch, par+ch)
+ dimnames(ans) <- list(names(par), c((1 - level)/2, 0.5+level/2))
}
- }
- if(!missing(parm))
- ans <- ans[parm,]
+ if(lambda)
+ {
+ if (length(z$coefficients) == length(z$lambda))
+ {
+ cat("\nNo confidence intervals for lambda when the model is just identified.\n")
+ return(NULL)
+ } else {
+ chl <- zs*se_parl
+ ans <- cbind(lamb - chl, lamb + chl)
+ dimnames(ans) <- list(names(lamb), c((1 - level)/2, 0.5 + level/2))
+ }
+ }
+ if(!missing(parm))
+ ans <- ans[parm,]
+ } else {
+ if(missing(parm))
+ parm <- names(object$coefficients)
+ type <- strsplit(type, "v")[[1]][2]
+ ntest <- paste("Confidence interval based on the inversion of the ", type, " test", sep="")
+ ans <- lapply(parm, function(w) .invTest(object, w, level = level, fact = fact, type=type, corr=corr))
+ ans <- do.call(rbind, ans)
+ if (!is.character(parm))
+ parm <- names(object$coefficients)[parm]
+ dimnames(ans) <- list(parm, c((1 - level)/2, 0.5+level/2))
+ }
+ ans <- list(test=ans,ntest=ntest)
+ class(ans) <- "confint"
ans
}
+print.confint <- function(x, digits = 5, ...)
+ {
+ cat("\n", x$ntest, sep="")
+ cat("\n#######################################\n")
+ print.default(format(x$test, digits = digits),
+ print.gap = 2, quote = FALSE)
+ invisible(x)
+ }
+
coef.gel <- function(object, lambda = FALSE, ...)
{
if(!lambda)
Modified: pkg/gmm/R/Methods.gmm.R
===================================================================
--- pkg/gmm/R/Methods.gmm.R 2015-10-28 19:28:30 UTC (rev 85)
+++ pkg/gmm/R/Methods.gmm.R 2015-10-29 18:48:28 UTC (rev 86)
@@ -13,7 +13,7 @@
summary.gmm <- function(object, ...)
- {
+ {
z <- object
se <- sqrt(diag(z$vcov))
par <- z$coefficients
@@ -21,57 +21,57 @@
ans <- list(met=z$met,kernel=z$kernel,algo=z$algo,call=z$call)
names(ans$met) <- "GMM method"
names(ans$kernel) <- "kernel for cov matrix"
-
+
ans$coefficients <- cbind(par,se, tval, 2 * pnorm(abs(tval), lower.tail = FALSE))
dimnames(ans$coefficients) <- list(names(z$coefficients),
- c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
+ c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
ans$stest <- specTest(z)
ans$algoInfo <- z$algoInfo
if(z$met=="cue")
- ans$cue <- object$cue
+ ans$cue <- object$cue
if (!is.null(object$initTheta))
- {
+ {
ans$initTheta <- object$initTheta
names(ans$initTheta) <- names(z$coefficients)
- }
+ }
ans$specMod <- object$specMod
ans$bw <- attr(object$w0,"Spec")$bw
ans$weights <- attr(object$w0,"Spec")$weights
if(object$infVcov == "iid")
- ans$kernel <- NULL
+ ans$kernel <- NULL
class(ans) <- "summary.gmm"
ans
- }
+ }
summary.tsls <- function(object, vcov = NULL, ...)
- {
+ {
if (!is.null(vcov))
- object$vcov=vcov
+ object$vcov=vcov
else
- object$vcov=vcov(object)
+ object$vcov=vcov(object)
ans <- summary.gmm(object)
ans$met <- paste(ans$met, "(Meat type = ", attr(object$vcov, "vcovType"), ")",sep="")
k <- object$dat$k
if (!is.null(object$fsRes))
- {
+ {
fstat <- vector()
fsRes <- object$fsRes
if (class(fsRes) == "listof")
{
- nendo <- length(fsRes)
- } else {
- nendo <- 1
- }
+ nendo <- length(fsRes)
+ } else {
+ nendo <- 1
+ }
if (nendo == 1)
{
- fstat[1] <- fsRes$fstatistic[1]
- df1 <- fsRes$fstatistic[2]
- df2 <- fsRes$fstatistic[3]
- } else {
- fstat[1] <- fsRes[[1]]$fstatistic[1]
- df1 <- fsRes[[1]]$fstatistic[2]
- df2 <- fsRes[[1]]$fstatistic[3]
- }
+ fstat[1] <- fsRes$fstatistic[1]
+ df1 <- fsRes$fstatistic[2]
+ df2 <- fsRes$fstatistic[3]
+ } else {
+ fstat[1] <- fsRes[[1]]$fstatistic[1]
+ df1 <- fsRes[[1]]$fstatistic[2]
+ df2 <- fsRes[[1]]$fstatistic[3]
+ }
if (nendo > 1){
for (i in 2:nendo)
fstat[i] <- fsRes[[i]]$fstatistic[1]
@@ -79,49 +79,49 @@
pvfstat <- 1-pf(fstat,df1, df2)
names(fstat) <- attr(fsRes,"Endo")
ans$fstatistic <- list(fstat = fstat, pvfstat = pvfstat, df1 = df1, df2 = df2)
- }
+ }
ans$specMod <- object$specMod
class(ans) <- "summary.tsls"
return(ans)
- }
+ }
print.summary.tsls <- function(x, digits = 5, ...)
- {
+ {
print.summary.gmm(x,digits)
if (!is.null(x$fstatistic))
- {
+ {
cat("\n First stage F-statistics: \n")
if(names(x$fstatistic$fstat)[1]=="(Intercept)")
- start=2
+ start=2
else
- start=1
+ start=1
for (i in start:length(x$fstatistic$fstat))
- cat(names(x$fstatistic$fstat)[i],
+ cat(names(x$fstatistic$fstat)[i],
": F(",x$fstatistic$df1,", ",x$fstatistic$df2,") = ",x$fstatistic$fstat[i],
" (P-Vavue = ",x$fstatistic$pvfstat[i],")\n")
- } else {
+ } else {
cat("\n No first stage F-statistics (just identified model)\n")
- }
- }
+ }
+ }
print.summary.gmm <- function(x, digits = 5, ...)
- {
+ {
cat("\nCall:\n")
cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
cat("\nMethod: ", x$met,"\n")
if (x$met=="cue")
- cat(" (",x$cue$message,")\n\n")
+ cat(" (",x$cue$message,")\n\n")
else
- cat("\n")
+ cat("\n")
if( !is.null(x$kernel))
- {
+ {
cat("Kernel: ", x$kernel)
if (!is.null(x$bw))
- cat("(with bw = ", round(x$bw,5),")\n\n")
+ cat("(with bw = ", round(x$bw,5),")\n\n")
else
- cat("\n\n")
- }
+ cat("\n\n")
+ }
cat("Coefficients:\n")
print.default(format(x$coefficients, digits=digits),
print.gap = 2, quote = FALSE)
@@ -131,28 +131,28 @@
print.gap = 2, quote = FALSE)
cat("\n")
if(!is.null(x$initTheta))
- {
+ {
cat("Initial values of the coefficients\n")
print(x$initTheta)
cat("\n")
- }
+ }
cat(x$specMod)
if(!is.null(x$algoInfo))
- {
+ {
cat("#############\n")
cat("Information related to the numerical optimization\n")
- }
+ }
if(!is.null(x$algoInfo$convergence))
- cat("Convergence code = ", x$algoInfo$convergence,"\n")
+ cat("Convergence code = ", x$algoInfo$convergence,"\n")
if(!is.null(x$algoInfo$counts))
- {
+ {
cat("Function eval. = ",x$algoInfo$counts[1],"\n")
cat("Gradian eval. = ",x$algoInfo$counts[2],"\n")
- }
+ }
if(!is.null(x$algoInfo$message))
- cat("Message: ",x$algoInfo$message,"\n")
+ cat("Message: ",x$algoInfo$message,"\n")
invisible(x)
- }
+ }
formula.gmm <- function(x, ...)
@@ -164,87 +164,90 @@
}
confint.gmm <- function(object, parm, level=0.95, ...)
- {
- z <- object
- se <- sqrt(diag(z$vcov))
- par <- z$coefficients
-
- zs <- qnorm((1-level)/2,lower.tail=FALSE)
- ch <- zs*se
- ans <- cbind(par-ch,par+ch)
- dimnames(ans) <- list(names(par),c((1-level)/2,0.5+level/2))
- if(!missing(parm))
- ans <- ans[parm,]
- ans
- }
-
+ {
+ ntest <- "Wald type confidence interval"
+ z <- object
+ se <- sqrt(diag(z$vcov))
+ par <- z$coefficients
+
+ zs <- qnorm((1-level)/2,lower.tail=FALSE)
+ ch <- zs*se
+ ans <- cbind(par-ch,par+ch)
+ dimnames(ans) <- list(names(par),c((1-level)/2,0.5+level/2))
+ if(!missing(parm))
+ ans <- ans[parm,]
+ ans <- list(test=ans, ntest=ntest)
+ class(ans) <- "confint"
+ ans
+ }
+
residuals.gmm <- function(object,...)
- {
+ {
if(is.null(object$model))
- stop("The residuals method is valid only for g=formula")
+ stop("The residuals method is valid only for g=formula")
object$residuals
- }
+ }
fitted.gmm <- function(object,...)
- {
+ {
if(is.null(object$model))
- stop("The residuals method is valid only for g=formula")
+ stop("The residuals method is valid only for g=formula")
object$fitted.value
- }
+ }
print.gmm <- function(x, digits=5, ...)
- {
+ {
cat("Method\n", x$met,"\n\n")
cat("Objective function value: ",x$objective,"\n\n")
print.default(format(coef(x), digits=digits),
print.gap = 2, quote = FALSE)
cat("\n")
if(!is.null(x$algoInfo$convergence))
- cat("Convergence code = ", x$algoInfo$convergence,"\n")
+ cat("Convergence code = ", x$algoInfo$convergence,"\n")
cat(x$specMod)
invisible(x)
- }
+ }
coef.gmm <- function(object,...) object$coefficients
vcov.gmm <- function(object,...) object$vcov
estfun.gmmFct <- function(x, y = NULL, theta = NULL, ...)
- {
+ {
if (is(x, "function"))
- {
+ {
gmat <- x(theta, y)
return(gmat)
- }
+ }
else
- return(x)
- }
+ return(x)
+ }
estfun.tsls <- function(x, ...)
- {
+ {
model.matrix(x)*c(residuals(x))
- }
+ }
model.matrix.tsls <- function(object, ...)
{
-dat <- object$dat
-ny <- dat$ny
-nh <- dat$nh
-k <- dat$k
-x <- dat$x
-n <- nrow(x)
-hm <- as.matrix(x[,(ny+k+1):(ny+k+nh)])
-xm <- as.matrix(x[,(ny+1):(ny+k)])
-xhat <- lm(xm~hm-1)$fitted
-assign <- 1:ncol(xhat)
-if (attr(object$terms,"intercept")==1)
+ dat <- object$dat
+ ny <- dat$ny
+ nh <- dat$nh
+ k <- dat$k
+ x <- dat$x
+ n <- nrow(x)
+ hm <- as.matrix(x[,(ny+k+1):(ny+k+nh)])
+ xm <- as.matrix(x[,(ny+1):(ny+k)])
+ xhat <- lm(xm~hm-1)$fitted
+ assign <- 1:ncol(xhat)
+ if (attr(object$terms,"intercept")==1)
assign <- assign-1
-attr(xhat,"assign") <- assign
-xhat
+ attr(xhat,"assign") <- assign
+ xhat
}
vcov.tsls <- function(object, type=c("Classical","HC0","HC1","HAC"), hacProp = list(), ...)
- {
+ {
type <- match.arg(type)
if (type == "Classical")
- {
+ {
sig <- sum(c(residuals(object))^2)/(nrow(object$dat$x)-object$dat$k)
ny <- object$dat$ny
nh <- object$dat$nh
@@ -254,16 +257,16 @@
Omega <- crossprod(hm)*sig/nrow(object$dat$x)
vcovType <- "Classical"
V <- solve(crossprod(object$G,solve(Omega,object$G)))/nrow(object$dat$x)
- }
+ }
else if (strtrim(type,2) == "HC")
- {
+ {
meat <- meatHC(object, type)
bread <- bread(object)
vcovType <- paste("HCCM: ", type, sep="")
V <- crossprod(bread, meat%*%bread)/nrow(object$dat$x)
- }
+ }
else
- {
+ {
object$centeredVcov <- TRUE
gt <- model.matrix(object)*c(residuals(object))
gt <- lm(gt~1)
@@ -273,27 +276,27 @@
vcovType <- paste("HAC: ", KType, sep="")
bread <- bread(object)
V <- crossprod(bread, meat%*%bread)/nrow(object$dat$x)
- }
+ }
attr(V, "vcovType") <- vcovType
return(V)
- }
+ }
estfun.gmm <- function(x, ...)
- {
- foc <- x$gt %*% x$w %*% x$G
- return(foc)
- }
+ {
+ foc <- x$gt %*% x$w %*% x$G
+ return(foc)
+ }
bread.gmm <- function(x, ...)
- {
- GWG <- crossprod(x$G, x$w %*% x$G)
- b <- try(solve(GWG), silent = TRUE)
- if (class(b) == "try-error")
- stop("The bread matrix is singular")
- return(b)
- }
+ {
+ GWG <- crossprod(x$G, x$w %*% x$G)
+ b <- try(solve(GWG), silent = TRUE)
+ if (class(b) == "try-error")
+ stop("The bread matrix is singular")
+ return(b)
+ }
bread.tsls <- function(x, ...)
- {
+ {
dat <- x$dat
ny <- dat$ny
nh <- dat$nh
@@ -304,7 +307,7 @@
xm <- as.matrix(x[,(ny+1):(ny+k)])
xhat <- lm(xm~hm-1)$fitted
solve(crossprod(xhat)/n)
- }
+ }
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2015-10-28 19:28:30 UTC (rev 85)
+++ pkg/gmm/R/getModel.R 2015-10-29 18:48:28 UTC (rev 86)
@@ -141,7 +141,7 @@
if (!is.null(dim(object$eqConst)))
stop("eqConst must be a vector which indicates which parameters to fix")
if (length(object$eqConst)>=length(object$tet0))
- stop("Too many constraints; use evalGel() if all coefficients are fixed")
+ stop("Too many constraints; use evalGel() if all coefficients are fixed")
if (is.character(object$eqConst))
{
if (is.null(names(object$tet0)))
@@ -149,7 +149,7 @@
if (any(!(object$eqConst %in% names(object$tet0))))
stop("Wrong coefficient names in eqConst")
object$eqConst <- sort(match(object$eqConst,names(object$tet0)))
- }
+ }
restTet <- object$tet0[object$eqConst]
obj$tet0 <- object$tet0[-object$eqConst]
object$eqConst <- cbind(object$eqConst,restTet)
@@ -173,123 +173,126 @@
getModel.baseGel <- function(object, ...)
{
- object$allArg <- c(object, list(...))
- if(is(object$g, "formula"))
- {
- dat <- getDat(object$g, object$x, data = object$data)
- k <- dat$k
- if (is.null(object$tet0))
- {
- if (!is.null(object$eqConst))
- stop("You have to provide tet0 with equality constrains")
- if (object$optfct == "optimize")
- stop("For optimize, you must provide the 2x1 vector tet0")
- res0 <- gmm(object$g, object$x, data=object$data)
- object$tet0 <- res0$coefficients
- if (object$smooth)
- gt <- res0$gt
- } else {
- if (object$optfct == "optimize")
- {
- if (k != 1)
- stop("optimize() is for univariate optimization")
- if (length(object$tet0) != 2)
- stop("For optimize(), tet0 must be a 2x1 vector")
- } else {
- if (k != length(object$tet0))
- stop("The number of starting values does not correspond to the number of regressors")
- }
- if (object$smooth)
- gt <- gmm(object$g, object$x, data=object$data)$gt
- }
- clname <- paste(class(object), ".modFormula", sep = "")
- object$gradv <- .DmomentFct
- object$gradvf <- FALSE
- object$x <- dat
- object$gform<-object$g
- 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])
- 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 {
- namesCoef <- namex
- object$namesgt <- nameh
- }
- if (is.null(names(object$tet0)))
- object$namesCoef <- namesCoef
- else
- object$namesCoef <- names(object$tet0)
- attr(object$x,"ModelType") <- "linear"
- attr(object$x, "k") <- k
- attr(object$x, "q") <- object$x$ny*object$x$nh
- attr(object$x, "n") <- NROW(object$x$x)
- } else {
- if (is.null(object$tet0))
- stop("You must provide the starting values tet0 for nonlinear moments")
- if(any(object$optfct == c("optim", "nlminb")))
- k <- length(object$tet0)
- else
- k <- 1
- attr(object$x,"ModelType") <- "nonlinear"
- attr(object$x, "momentfct") <- object$g
- attr(object$x, "k") <- k
- attr(object$x, "q") <- NCOL(object$g(object$tet0, object$x))
- attr(object$x, "n") <- NROW(object$x)
- if(is.null(names(object$tet0)))
- object$namesCoef <- paste("Theta[" ,1:attr(object$x, "k"), "]", sep = "")
- else
- object$namesCoef <- names(object$tet0)
- if (!is.function(object$gradv) | object$smooth)
- {
- object$gradvf <- FALSE
- } else {
- attr(object$x, "gradv") <- object$gradv
- object$gradvf <- TRUE
- }
- object$gradv <- .DmomentFct
- if (object$smooth)
- gt <- gmm(object$g, object$x, object$tet0, wmatrix = "ident", ...)$gt
- clname <- paste(class(object), ".mod", sep = "")
- }
- if (object$smooth)
- {
- if (is.function(object$gradv))
- warning("The provided gradv is not used when smooth=TRUE",
- call. = FALSE)
- if(object$kernel == "Truncated")
+ object$allArg <- c(object, list(...))
+ object$allArg$weights <- NULL
+ object$allArg$call <- NULL
+ if(is(object$g, "formula"))
{
- object$wkernel <- "Bartlett"
- object$k1 <- 2
- object$k2 <- 2
+ dat <- getDat(object$g, object$x, data = object$data)
+ k <- dat$k
+ if (is.null(object$tet0))
+ {
+ if (!is.null(object$eqConst))
+ stop("You have to provide tet0 with equality constrains")
+ if (object$optfct == "optimize")
+ stop("For optimize, you must provide the 2x1 vector tet0")
+ res0 <- gmm(object$g, object$x, data=object$data)
+ object$tet0 <- res0$coefficients
+ if (object$smooth)
+ gt <- res0$gt
+ } else {
+ if (object$optfct == "optimize")
+ {
+ if (k != 1)
+ stop("optimize() is for univariate optimization")
+ if (length(object$tet0) != 2)
+ stop("For optimize(), tet0 must be a 2x1 vector")
+ } else {
+ if (k != length(object$tet0))
+ stop("The number of starting values does not correspond to the number of regressors")
+ }
+ if (object$smooth)
+ gt <- gmm(object$g, object$x, data=object$data)$gt
+ }
+ clname <- paste(class(object), ".modFormula", sep = "")
+ object$gradv <- .DmomentFct
+ object$gradvf <- FALSE
+ object$x <- dat
+ object$gform<-object$g
+ 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])
+ 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 {
+ namesCoef <- namex
+ object$namesgt <- nameh
+ }
+ if (is.null(names(object$tet0)))
+ object$namesCoef <- namesCoef
+ else
+ object$namesCoef <- names(object$tet0)
+ attr(object$x,"ModelType") <- "linear"
+ attr(object$x, "k") <- k
+ attr(object$x, "q") <- object$x$ny*object$x$nh
+
+ attr(object$x, "n") <- NROW(object$x$x)
+ } else {
+ if (is.null(object$tet0))
+ stop("You must provide the starting values tet0 for nonlinear moments")
+ if(any(object$optfct == c("optim", "nlminb")))
+ k <- length(object$tet0)
+ else
+ k <- 1
+ attr(object$x,"ModelType") <- "nonlinear"
+ attr(object$x, "momentfct") <- object$g
+ attr(object$x, "k") <- k
+ attr(object$x, "q") <- NCOL(object$g(object$tet0, object$x))
+ attr(object$x, "n") <- NROW(object$x)
+ if(is.null(names(object$tet0)))
+ object$namesCoef <- paste("Theta[" ,1:attr(object$x, "k"), "]", sep = "")
+ else
+ object$namesCoef <- names(object$tet0)
+ if (!is.function(object$gradv) | object$smooth)
+ {
+ object$gradvf <- FALSE
+ } else {
+ attr(object$x, "gradv") <- object$gradv
+ object$gradvf <- TRUE
+ }
+ object$gradv <- .DmomentFct
+ if (object$smooth)
+ gt <- gmm(object$g, object$x, object$tet0, wmatrix = "ident", ...)$gt
+ clname <- paste(class(object), ".mod", sep = "")
}
- if(object$kernel == "Bartlett")
+ if (object$smooth)
{
- object$wkernel <- "Parzen"
+ if (is.function(object$gradv))
+ warning("The provided gradv is not used when smooth=TRUE",
+ call. = FALSE)
+ if(object$kernel == "Truncated")
+ {
+ object$wkernel <- "Bartlett"
+ object$k1 <- 2
+ object$k2 <- 2
+ }
+ if(object$kernel == "Bartlett")
+ {
+ object$wkernel <- "Parzen"
+ object$k1 <- 1
+ object$k2 <- 2/3
+ }
+ gt <- scale(gt, scale=FALSE)
+ class(gt) <- "gmmFct"
+ if (is.function(object$bw))
+ object$bwVal <- object$bw(gt, kernel = object$wkernel, prewhite = object$prewhite,
+ ar.method = object$ar.method, approx = object$approx)
+ else
+ object$bwVal <- object$bw
+ object$w <- smoothG(gt, bw = object$bwVal, kernel = object$kernel, tol = object$tol_weights)$kern_weights
+ attr(object$x,"smooth") <- list(bw=object$bwVal, w=object$w, kernel=object$kernel)
+ } else {
object$k1 <- 1
- object$k2 <- 2/3
+ object$k2 <- 1
+ object$w <- kernel(1)
+ object$bwVal <- 1
}
- gt <- scale(gt, scale=FALSE)
- class(gt) <- "gmmFct"
- if (is.function(object$bw))
- object$bwVal <- object$bw(gt, kernel = object$wkernel, prewhite = object$prewhite,
- ar.method = object$ar.method, approx = object$approx)
- else
- object$bwVal <- object$bw
- object$w <- smoothG(gt, bw = object$bwVal, kernel = object$kernel, tol = object$tol_weights)$kern_weights
- attr(object$x,"smooth") <- list(bw=object$bwVal, w=object$w, kernel=object$kernel)
- } else {
- object$k1 <- 1
- object$k2 <- 1
- object$w <- kernel(1)
- object$bwVal <- 1
+ object$g <- .momentFct
+ object$CGEL <- object$alpha
+ object$typeDesc <- object$type
+ class(object) <- clname
+ return(object)
}
- object$g <- .momentFct
- object$CGEL <- object$alpha
- object$typeDesc <- object$type
- class(object) <- clname
- return(object)
-}
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2015-10-28 19:28:30 UTC (rev 85)
+++ pkg/gmm/R/momentEstim.R 2015-10-29 18:48:28 UTC (rev 86)
@@ -599,8 +599,10 @@
{
P <- object
g <- P$g
- q <- attr(object$x, "q")
- n <- attr(object$x, "n")
+ q <- attr(P$x, "q")
+ n <- attr(P$x, "n")
+ k <- attr(P$x, "k")
+ df <- q-k*P$x$ny
l0Env <- new.env()
assign("l0",rep(0,q),envir=l0Env)
if (!P$constraint)
@@ -638,6 +640,7 @@
z$CGEL <- P$CGEL
z$typeDesc <- P$typeDesc
z$specMod <- P$specMod
+ z$df <- df
names(z$coefficients) <- object$namesCoef
if (!is.null(object$namesgt))
@@ -659,7 +662,8 @@
ncoef[eqConst[,1]] <- rownames(eqConst)
names(coef) <- ncoef
z$coefficients <- coef
- attr(P$x, "k") <- attr(P$x, "k") + nrow(eqConst)
+ attr(P$x, "k") <- attr(P$x, "k") + nrow(eqConst)
+ z$df <- z$df - nrow(eqConst)
attr(P$x,"eqConst") <- NULL
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 86
More information about the Gmm-commits
mailing list