[Gmm-commits] r88 - in pkg/gmm: . R data man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 11 22:44:07 CET 2015
Author: chaussep
Date: 2015-11-11 22:44:07 +0100 (Wed, 11 Nov 2015)
New Revision: 88
Added:
pkg/gmm/data/wage.rda
pkg/gmm/man/sysGmm.Rd
pkg/gmm/man/wage.Rd
Modified:
pkg/gmm/NAMESPACE
pkg/gmm/NEWS
pkg/gmm/R/FinRes.R
pkg/gmm/R/Methods.gmm.R
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/momentEstim.R
pkg/gmm/R/specTest.R
pkg/gmm/man/getModel.Rd
pkg/gmm/man/gmm.Rd
pkg/gmm/man/momentEstim.Rd
pkg/gmm/man/print.Rd
pkg/gmm/man/summary.Rd
Log:
Added functions to estimate system of equations (in development
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/NAMESPACE 2015-11-11 21:44:07 UTC (rev 88)
@@ -12,13 +12,18 @@
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, print.confint)
+ momentEstim.baseGel.eval, evalGel, confint.gmm, print.confint, sysGmm, getModel.sysGmm, momentEstim.sysGmm.twoStep.formula,
+ momentEstim.tsls.twoStep.formula, getModel.tsls, summary.sysGmm, print.sysGmm, print.summary.sysGmm,
+ sur, threeSLS, randEffect, five)
S3method(summary, gmm)
+S3method(summary, sysGmm)
S3method(summary, tsls)
S3method(summary, gel)
S3method(print, gmm)
+S3method(print, sysGmm)
S3method(print, summary.gmm)
+S3method(print, summary.sysGmm)
S3method(coef, gmm)
S3method(vcov, gmm)
S3method(vcov, tsls)
@@ -45,11 +50,15 @@
S3method(FinRes, baseGmm.res)
S3method(getModel, baseGmm)
S3method(getModel, baseGel)
+S3method(getModel, sysGmm)
S3method(getModel, constGmm)
S3method(getModel, constGel)
+S3method(getModel, tsls)
S3method(momentEstim, baseGmm.twoStep)
S3method(momentEstim, baseGmm.eval)
S3method(momentEstim, baseGmm.twoStep.formula)
+S3method(momentEstim, tsls.twoStep.formula)
+S3method(momentEstim, sysGmm.twoStep.formula)
S3method(momentEstim, baseGmm.iterative.formula)
S3method(momentEstim, baseGmm.iterative)
S3method(momentEstim, baseGmm.cue.formula)
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/NEWS 2015-11-11 21:44:07 UTC (rev 88)
@@ -1,5 +1,6 @@
Changes in version 1.6-0
+o tsls is now a real 2-stage least square. Before it was a 2-step optimal GMM with HCCM weighting matrix.
o Fixed a typo in FinRes.R file. It was preventing to compute the proper vcov matrix for a very special case (fixed weights)
o Added Helinger distance and Exponentially tilted Helinger estimator to gel()
o Fixed the LR test of ETEL in gel()
@@ -20,6 +21,9 @@
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.
o On the process to create a sysGmm for system og linear equations. The function sysGmm() will do the job. Not yet working.
+o getDat is modified so that the name of the dependent variable is the one included in the formula.
+o A set of functions have been added to estimate system of linear equations (this is a beta version not tested yet).
+o A set of labor data as been added to test the sysGmm function
Changes in version 1.5-2
Modified: pkg/gmm/R/FinRes.R
===================================================================
--- pkg/gmm/R/FinRes.R 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/R/FinRes.R 2015-11-11 21:44:07 UTC (rev 88)
@@ -23,7 +23,6 @@
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
@@ -112,4 +111,76 @@
return(z)
}
+FinRes.sysGmm.res <- function(z, object, ...)
+ {
+ P <- object
+ x <- z$dat
+ z$G <- z$gradv(z$coefficients, x)
+ n <- z$n
+ G <- z$G
+ v <- .weightFct_Sys(z$coefficient, x, P$vcov)
+ nk <- z$k
+ 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,nk,nk)
+ 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,nk,nk)
+ warning("The covariance matrix of the coefficients is singular")
+ }
+ } else {
+ if (is.null(P$weightsMatrix))
+ w <- .weightFct_Sys(z$coefficients, x, "ident")
+ 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, nk, nk)
+ warning("The covariance matrix of the coefficients is singular")
+ } else {
+ z$vcov <- T1%*%v%*%t(T1)/n
+ }
+ }
+ if (attr(x, "sysInfo")$commonCoef)
+ dimnames(z$vcov) <- list(P$namesCoef[[1]], P$namesCoef[[1]])
+ else
+ dimnames(z$vcov) <- list(do.call(c, P$namesCoef), do.call(c, P$namesCoef))
+ z$call <- P$call
+ if(is.null(P$weightsMatrix))
+ {
+ if(P$wmatrix == "ident")
+ {
+ z$w <- .weightFct_Sys(z$coefficient, x, "ident")
+ } 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 <- paste("System of Equations: ", P$type, sep="")
+ z$kernel <- P$kernel
+ z$coefficients <- c(z$coefficients)
+ class(z) <- c("sysGmm", "gmm")
+ return(z)
+ }
+
Modified: pkg/gmm/R/Methods.gmm.R
===================================================================
--- pkg/gmm/R/Methods.gmm.R 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/R/Methods.gmm.R 2015-11-11 21:44:07 UTC (rev 88)
@@ -43,6 +43,50 @@
ans
}
+summary.sysGmm <- function(object, ...)
+ {
+ z <- object
+ se <- sqrt(diag(z$vcov))
+ k <- attr(z$dat, "k")
+ if (attr(z$dat, "sysInfo")$commonCoef)
+ {
+ seList <- rep(list(se), length(z$dat))
+ } else {
+ seList <- list()
+ for (i in 1:length(z$dat))
+ {
+ seList[[i]] <- se[1:k[[i]]]
+ se <- se[-(1:k[[i]])]
+ }
+ }
+ par <- z$coefficients
+ tval <- lapply(1:length(z$dat), function(i) par[[i]]/seList[[i]])
+ 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 <- lapply(1:length(z$dat), function(i) cbind(par[[i]],seList[[i]], tval[[i]], 2 * pnorm(abs(tval[[i]]), lower.tail = FALSE)))
+ ans$stest <- specTest(z)
+ ans$algoInfo <- z$algoInfo
+ ans$initTheta <- object$initTheta
+
+ for (i in 1:length(z$dat))
+ {
+ dimnames(ans$coefficients[[i]]) <- list(names(z$coefficients[[i]]),
+ c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
+ names(ans$initTheta[[i]]) <- names(z$coefficients[[i]])
+ }
+ ans$specMod <- object$specMod
+ ans$bw <- attr(object$w0,"Spec")$bw
+ ans$weights <- attr(object$w0,"Spec")$weights
+ ans$Sysnames <- names(z$dat)
+ ans$met <- attr(object$dat, "sysInfo")$typeDesc
+ if(object$infVcov != "HAC")
+ ans$kernel <- NULL
+ class(ans) <- "summary.sysGmm"
+ ans
+ }
+
summary.tsls <- function(object, vcov = NULL, ...)
{
if (!is.null(vcov))
@@ -155,6 +199,63 @@
}
+print.summary.sysGmm <- function(x, digits = 5, ...)
+ {
+ cat("\nCall:\n")
+ cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
+ cat("Method\n", x$met,"\n\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")
+ else
+ cat("\n\n")
+ }
+ cat("Coefficients:\n")
+ m <- length(x$coefficients)
+ for (i in 1:m)
+ {
+ cat(x$Sysnames[i], "\n")
+ cat("#########\n")
+ print.default(format(x$coefficients[[i]], digits=digits),
+ print.gap = 2, quote = FALSE)
+ cat("\n")
+ }
+ cat(x$stest$ntest,"\n")
+ print.default(format(x$stest$test, digits=digits),
+ print.gap = 2, quote = FALSE)
+ cat("\n")
+ if(!is.null(x$initTheta))
+ {
+ cat("Initial values of the coefficients\n")
+ for (i in 1:m)
+ {
+ cat(x$Sysnames[i], "\n")
+ print(x$initTheta[[i]])
+ }
+ 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")
+ 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")
+ invisible(x)
+ }
+
+
formula.gmm <- function(x, ...)
{
if(is.null(x$terms))
@@ -208,6 +309,23 @@
invisible(x)
}
+print.sysGmm <- function(x, digits=5, ...)
+ {
+ cat("Method\n", attr(x$dat, "sysInfo")$typeDesc,"\n\n")
+ cat("Objective function value: ",x$objective,"\n\n")
+ for (i in 1:length(x$coefficients))
+ {
+ cat(names(x$dat)[[i]], ": \n")
+ print.default(format(coef(x)[[i]], digits=digits),
+ print.gap = 2, quote = FALSE)
+ }
+ cat("\n")
+ if(!is.null(x$algoInfo$convergence))
+ cat("Convergence code = ", x$algoInfo$convergence,"\n")
+ cat(x$specMod)
+ invisible(x)
+ }
+
coef.gmm <- function(object,...) object$coefficients
vcov.gmm <- function(object,...) object$vcov
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/R/getModel.R 2015-11-11 21:44:07 UTC (rev 88)
@@ -16,9 +16,16 @@
UseMethod("getModel")
}
+getModel.tsls <- function(object, ...)
+ {
+ obj <- getModel.baseGmm(object, ...)
+ return(obj)
+ }
+
getModel.sysGmm <- function(object, ...)
{
object$allArg <- c(object, list(...))
+ object$formula <- list(g=object$g, h=object$h)
if (!is.list(object$g))
stop("g must be list of formulas")
if (length(object$g) == 1)
@@ -37,52 +44,52 @@
}
if (length(object$h) == 1)
{
- dat <- lapply(1:length(object$g), function(i) try(getDat(object$g[[i]], object$h[[1]], data = object$data), silent=TRUE))
+ object$h <- rep(object$h, length(object$g))
typeDesc <- "System Gmm with common instruments"
+ typeInst <- "Common"
} else {
if (length(object$h) != length(object$g))
- stop("The number of formulas in h should be the same as the number of formulas in g, \nunless the instruments are the same in each equation, \nin which case the number of equation in h should be 1")
- dat <- lapply(1:length(object$g), function(i) try(getDat(object$g[[i]], object$h[[i]], data = object$data), silent=TRUE))
+ stop("The number of formulas in h should be the same as the number of formulas in g, \nunless the instruments are the same in each equation, \nin which case the number of equations in h should be 1")
typeDesc <- "System Gmm"
+ typeInst <- "nonCommon"
}
+ if (object$commonCoef)
+ typeDesc <- paste(typeDesc, " (Common Coefficients)")
+ dat <- lapply(1:length(object$g), function(i) try(getDat(object$g[[i]], object$h[[i]], data = object$data), silent=TRUE))
+ chk <- sapply(1:length(dat), function(i) class(dat[[i]]) == "try-error")
+ chk <- which(chk)
+ mess <- vector()
+ for (i in chk)
+ {
+ mess <- paste(mess, "\nSystem of equations:", i, "\n###############\n", sep="")
+ mess <- paste(mess, attr(dat[[i]], "condition")[[1]])
+ }
+ if (length(chk)>0)
+ stop(mess)
if (is.null(names(object$g)))
names(dat) <- paste("System_", 1:length(dat), sep="")
else
names(dat) <- names(object$g)
-#### Need to set the gradiant function #### object$gradv <- .DmomentFct
- if (length(object$h) == 1)
- dat <- lapply(1:length(object$g), function(i) try(getDat(object$g[[i]], object$h[[1]], data = object$data), silent=TRUE))
- for (i in 1:length(dat))
- {
- if (class(dat[[i]]) == "try-error")
- {
- cat("\nSystem of equations:", i, "\n###############\n", sep="")
- stop(dat[[i]])
- }
- }
+ object$gradv <- .DmomentFct_Sys
+ object$formula <- list(g=object$g, h=object$h)
if (!all(sapply(1:length(dat), function(i) dat[[i]]$ny == 1)))
stop("The number of dependent variable must be one in every equation")
-#### What to do with fixed weights???
-# if(is.null(object$weightsMatrix))
-# {
-# clname <- paste(class(object), ".", object$type, ".formula", sep = "")
-# } else {
-# clname <- "fixedW.formula"
-# object$type <- "One step GMM with fixed W"
-# }
- clname <- "sysGmm.formula"
+ clname <- "sysGmm.twoStep.formula"
object$x <- dat
namex <- lapply(1:length(dat), function(i) colnames(dat[[i]]$x[,2:(1+dat[[i]]$k), drop=FALSE]))
- nameh <- lapply(1:length(dat), function(i) colnames(dat[[i]]$x[,(2+dat$k):(1+dat[[i]]$k+dat[[i]]$nh), drop=FALSE]))
+ nameh <- lapply(1:length(dat), function(i) colnames(dat[[i]]$x[,(2+dat[[i]]$k):(1+dat[[i]]$k+dat[[i]]$nh), drop=FALSE]))
namey <- lapply(1:length(dat), function(i) colnames(dat[[i]]$x[,1, drop=FALSE]))
- object$namesCoef <- do.call(c, namex)
- object$namesgt <- do.call(c, nameh)
- object$namesy <- do.call(c, namey)
+ object$namesCoef <- namex
+ namesgt <- lapply(1:length(dat), function(i) paste(namey[[i]], "_", nameh[[i]], sep=""))
+ object$namesgt <- namesgt
+ object$namesy <- namey
attr(object$x,"ModelType") <- "linear"
- attr(object$x, "k") <- length(object$namesCoef)
- attr(object$x, "q") <- length(object$namesgt)
- attr(object$x, "n") <- NROW(object$x[[1]]$x)
+ #for (i in 1:length(object$x))
+ # attr(object$x[[i]], c("linear")) <- attr(object$x, "modelType")
+ attr(object$x, "k") <- lapply(1:length(dat), function(i) length(object$namesCoef[[i]]))
+ attr(object$x, "q") <- lapply(1:length(dat), function(i) length(object$namesgt[[i]]))
+ attr(object$x, "n") <- lapply(1:length(dat), function(i) nrow(object$x[[i]]$x))
object$TypeGmm <- class(object)
attr(object$x, "weight") <- list(w=object$weightsMatrix,
centeredVcov=object$centeredVcov)
@@ -92,7 +99,17 @@
ar.method = object$ar.method,
approx = object$approx, tol = object$tol)
attr(object$x, "weight")$vcov <- object$vcov
- object$g <- .momentFct
+ k <- lapply(1:length(dat), function(i) dat[[i]]$k)
+ q <- lapply(1:length(dat), function(i) dat[[i]]$nh)
+ df <- lapply(1:length(dat), function(i) dat[[i]]$nh-dat[[i]]$k)
+ if (object$commonCoef)
+ {
+ k <- do.call(c,k)
+ if (!all(k==k[1]))
+ stop("In a common coefficient model, the number of regressors is the same in each equation")
+ }
+ attr(object$x, "sysInfo") <- list(k=k, df=df, q=q, typeInst=typeInst, typeDesc=typeDesc, commonCoef=object$commonCoef)
+ object$g <- .momentFct_Sys
class(object) <- clname
return(object)
}
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/R/gmm.R 2015-11-11 21:44:07 UTC (rev 88)
@@ -11,7 +11,7 @@
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
-gmm <- function(g,x,t0=NULL,gradv=NULL, type=c("twoStep","cue","iterative"), wmatrix = c("optimal","ident"), vcov=c("HAC","iid","TrueFixed"),
+gmm <- function(g,x,t0=NULL,gradv=NULL, type=c("twoStep","cue","iterative"), wmatrix = c("optimal","ident"), vcov=c("HAC","MDS","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, 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,
@@ -46,39 +46,81 @@
z
}
-sysGmm <- function(g, h, type=c("twoStep","cue","iterative"), wmatrix = c("optimal","ident"), vcov=c("HAC","mdf","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, itermax=100,optfct=c("optim","optimize","nlminb", "constrOptim"),
- model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, weightsMatrix = NULL, traceIter = FALSE, data, eqConst = NULL,
- eqConstFullVcov = FALSE, ...)
+five <- function(g, h, commonCoef = FALSE, data = NULL)
{
- type <- match.arg(type)
+ res <- sysGmm(g, h, wmatrix = "optimal", vcov = "CondHom", commonCoef=commonCoef, data=data)
+ attr(res$dat, "sysInfo")$typeDesc <- "Full-Information IV (FIVE)"
+ res$call <- match.call()
+ res
+ }
+
+threeSLS <- function(g, h, commonCoef = FALSE, data = NULL)
+ {
+ if (!is(h, "formula"))
+ if (length(h) != 1)
+ stop("In 3SLS, h is a single equation since the same instruments are used in each equation")
+ res <- sysGmm(g, h, vcov = "CondHom", commonCoef=commonCoef, data=data)
+ attr(res$dat, "sysInfo")$typeDesc <- "3-stage Least Squares"
+ res$call <- match.call()
+ res
+ }
+sur <- function(g, commonCoef = FALSE, data = NULL)
+ {
+ if (!is.list(g))
+ stop("g must be list of formulas")
+ if (length(g) == 1)
+ stop("For single equation GMM, use the function gmm()")
+ if (!all(sapply(1:length(g), function(i) is(g[[i]], "formula"))))
+ stop("g must be a list of formulas")
+ tm <- lapply(1:length(g), function(i) terms(g[[i]]))
+ reg <- lapply(1:length(g), function(i) attr(tm[[i]], "term.labels"))
+ reg <- unique(do.call(c,reg))
+ h <- paste(reg, collapse = "+")
+ if (all(sapply(1:length(g), function(i) attr(tm[[i]], "intercept") == 0)))
+ h <- paste(h, "-1")
+ h <- as.formula(paste("~", h))
+ res <- sysGmm(g, h, vcov = "CondHom", commonCoef=commonCoef, data=data)
+ attr(res$dat, "sysInfo")$typeDesc <- "Seemingly Unrelated Regression (SUR)"
+ res$call <- match.call()
+ res
+ }
+
+randEffect <- function(g, data = NULL)
+ {
+ res <- sur(g, commonCoef = TRUE, data = data)
+ attr(res$dat, "sysInfo")$typeDesc <- "Random Effect Estimator"
+ res$call <- match.call()
+ res
+ }
+
+
+sysGmm <- function(g, h, wmatrix = c("optimal","ident"),
+ vcov=c("MDS", "HAC", "CondHom", "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, commonCoef = FALSE)
+ {
kernel <- match.arg(kernel)
vcov <- match.arg(vcov)
wmatrix <- match.arg(wmatrix)
- optfct <- match.arg(optfct)
+ TypeGmm = "sysGmm"
- if (!is.null(eqConst))
- TypeGmm <- "constSysGmm"
- else
- TypeGmm = "sysGmm"
-
if(vcov=="TrueFixed" & is.null(weightsMatrix))
stop("TrueFixed vcov only for fixed weighting matrix")
if(!is.null(weightsMatrix))
wmatrix <- "optimal"
if(missing(data))
data<-NULL
- all_args<-list(data = data, g = g, h = h, type = type, wmatrix = wmatrix, vcov = vcov, kernel = kernel,
+ all_args<-list(data = data, g = g, h = h, wmatrix = wmatrix, vcov = vcov, kernel = kernel,
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)
+ weightsMatrix = weightsMatrix, centeredVcov = centeredVcov, tol = tol,
+ model = model, X = X, Y = Y, call = match.call(), commonCoef=commonCoef)
class(all_args)<-TypeGmm
- Model_info<-getModel(all_args, ...)
- z <- momentEstim(Model_info, ...)
+ Model_info<-getModel(all_args)
+ z <- momentEstim(Model_info)
z <- FinRes(z, Model_info)
- z
+ return(z)
}
@@ -118,7 +160,7 @@
{
if(class(g) != "formula")
stop("2SLS is for linear models expressed as formula only")
- ans <- gmm(g,x,data=data,vcov="iid")
+ ans <- gmm(g,x,data=data,vcov="iid", TypeGmm="tsls")
ans$met <- "Two Stage Least Squares"
ans$call <- match.call()
class(ans) <- c("tsls","gmm")
@@ -161,6 +203,9 @@
mf <- eval(mf, parent.frame())
mt <- attr(mf, "terms")
y <- as.matrix(model.response(mf, "numeric"))
+ namey <- as.character(formula)[2]
+ if (ncol(y)>1)
+ namey <- paste(namey, ".", 1:ncol(y), sep="")
xt <- as.matrix(model.matrix(mt, mf, NULL))
n <- NROW(y)
if (inherits(h,'formula'))
@@ -199,12 +244,7 @@
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)))
- {
- if (ny>1)
- colnames(y) <- paste("y",1:ncol(y),sep="")
- if (ny == 1)
- colnames(y) <- "y"
- }
+ colnames(y) <- namey
rownames(xt) <- rownames(y)
rownames(h) <- rownames(y)
x <- cbind(y,xt,h)
@@ -219,7 +259,6 @@
return(list(x=x,nh=nh,ny=ny,k=k,mf=mf,mt=mt,cl=cl))
}
-
.tetlin <- function(dat, w, type=NULL)
{
x <- dat$x
@@ -372,6 +411,21 @@
return(obj)
}
+.momentFct_Sys <- function(tet, dat)
+ {
+ q <- length(dat)
+ f <- function(i, dat)
+ {
+ d <- dat[[i]]
+ attr(d, "eqConst") <- attr(dat, "eqConst")
+ attr(d, "ModelType") <- attr(dat, "ModelType")
+ attr(d,"momentfct") <- attr(dat,"momentfct")
+ attr(d, "smooth") <- attr(dat, "smooth")
+ .momentFct(tet[[i]], d)
+ }
+ mom <- lapply(1:q, function(i) f(i, dat))
+ do.call(cbind, mom)
+ }
.momentFct <- function(tet, dat)
{
@@ -408,6 +462,25 @@
return(as.matrix(gt))
}
+.DmomentFct_Sys <- function(tet, dat, pt=NULL)
+ {
+ q <- length(dat)
+ f <- function(i, dat)
+ {
+ d <- dat[[i]]
+ attr(d, "eqConst") <- attr(dat, "eqConst")
+ attr(d, "ModelType") <- attr(dat, "ModelType")
+ attr(d,"momentfct") <- attr(dat,"momentfct")
+ attr(d, "smooth") <- attr(dat, "smooth")
+ .DmomentFct(tet[[i]], d, pt)
+ }
+ dmom <- lapply(1:q, function(i) f(i, dat))
+ if (attr(dat, "sysInfo")$commonCoef)
+ do.call(rbind,dmom)
+ else
+ .diagMatrix(dmom)
+ }
+
.DmomentFct <- function(tet, dat, pt=NULL)
{
if (!is.null(attr(dat, "eqConst")))
@@ -441,7 +514,7 @@
return(as.matrix(dgb))
}
-.weightFct <- function(tet, dat, type=c("HAC", "iid", "ident", "fct", "fixed"))
+.weightFct <- function(tet, dat, type=c("HAC", "MDS", "iid", "ident", "fct", "fixed"))
{
type <- match.arg(type)
if (type == "fixed")
@@ -466,11 +539,24 @@
w <- .myKernHAC(gt, obj)
attr(w, "inv") <- TRUE
}
- if (type == "iid")
+ if (type == "MDS")
{
w <- crossprod(gt)/n
attr(w, "inv") <- TRUE
}
+ if (type == "iid")
+ {
+ if ((attr(dat, "ModelType") == "linear") & (dat$ny == 1))
+ {
+ e <- .residuals(tet, dat)$residuals
+ sig <- mean(scale(e,scale=FALSE)^2)
+ z <- dat$x[,(1+dat$ny+dat$k):ncol(dat$x)]
+ w <- sig*crossprod(z)/length(e)
+ } else {
+ w <- crossprod(gt)/n
+ }
+ attr(w, "inv") <- TRUE
+ }
if (type == "fct")
{
w <- attr(dat, "weight")$fct(gt, attr(dat, "weight")$fctArg)
@@ -479,6 +565,61 @@
return(w)
}
+
+.weightFct_Sys<- function(tet, dat, type=c("MDS", "HAC", "CondHom", "ident", "fct", "fixed"))
+ {
+ type <- match.arg(type)
+ if (type == "fixed")
+ {
+ w <- attr(dat, "weight")$w
+ attr(w, "inv") <- FALSE
+ } else if (type == "ident") {
+ w <- diag(attr(dat, "q"))
+ attr(w, "inv") <- FALSE
+ } else {
+ if (type == "HAC")
+ {
+ gt <- .momentFct_Sys(tet,dat)
+ if(attr(dat, "weight")$centeredVcov)
+ gt <- residuals(lm(gt~1))
+ n <- NROW(gt)
+ obj <- attr(dat, "weight")
+ obj$centeredVcov <- FALSE
+ w <- .myKernHAC(gt, obj)
+ attr(w, "inv") <- TRUE
+ }
+ if (type == "MDS")
+ {
+ gt <- .momentFct_Sys(tet,dat)
+ n <- NROW(gt)
+ if(attr(dat, "weight")$centeredVcov)
+ gt <- scale(gt, scale=FALSE)
+ w <- crossprod(gt)/n
+ attr(w, "inv") <- TRUE
+ }
+ if (type == "CondHom")
+ {
+ e <- lapply(1:length(dat), function(i) .residuals(tet[[i]], dat[[i]])$residuals)
+ e <- do.call(cbind, e)
+ Sig <- crossprod(scale(e, scale=FALSE))/nrow(e)
+ Z <- lapply(1:length(dat), function(i) dat[[i]]$x[,(2+dat[[i]]$k):ncol(dat[[i]]$x)])
+ Z <- do.call(cbind, Z)
+ w <- crossprod(Z)/nrow(e)
+ for (i in 1:length(dat))
+ for (j in 1:length(dat))
+ {
+ s1 <- 1+(i-1)*dat[[i]]$nh
+ e1 <- i*dat[[i]]$nh
+ s2 <- 1+(j-1)*dat[[j]]$nh
+ e2 <- j*dat[[j]]$nh
+ w[s1:e1, s2:e2] <- w[s1:e1, s2:e2]*Sig[i,j]
+ }
+ attr(w, "inv") <- TRUE
+ }
+ }
+ return(w)
+ }
+
.residuals <- function(tet, dat)
{
if (!is.null(attr(dat, "eqConst")))
@@ -497,5 +638,24 @@
e <- y-yhat
return(list(residuals=e, yhat=yhat))
}
+
+.diagMatrix <- function(xlist)
+ {
+ # Create block diagonal matrix from matrices with the same number of rows.
+ m <- length(xlist)
+ n <- NROW(xlist[[1]])
+ l <- sapply(1:m, function(i) dim(as.matrix(xlist[[i]])))
+ dimX <- rowSums(l)
+ X <- matrix(0, dimX[1], dimX[2])
+ for (i in 1:m)
+ {
+ s1 <- 1 + (i-1)*n
+ e1 <- n*i
+ s2 <- 1 + sum(l[2,][-(i:m)])
+ e2 <- sum(l[2,][1:i])
+ X[s1:e1, s2:e2] <- xlist[[i]]
+ }
+ X
+ }
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2015-11-02 18:58:56 UTC (rev 87)
+++ pkg/gmm/R/momentEstim.R 2015-11-11 21:44:07 UTC (rev 88)
@@ -17,6 +17,102 @@
UseMethod("momentEstim")
}
+momentEstim.sysGmm.twoStep.formula <- function(object, ...)
+ {
+ dat <- object$x
+ if (attr(dat, "sysInfo")$commonCoef)
+ {
+ ## Getting first step estimate ##
+ y <- lapply(1:length(dat), function(i) dat[[i]]$x[,1])
+ y <- do.call(c, y)
+ x <- lapply(1:length(dat), function(i) dat[[i]]$x[,2:(dat[[i]]$k+1)])
+ x <- do.call(rbind, x)
+ z <- lapply(1:length(dat), function(i)
+ dat[[i]]$x[,(2+dat[[i]]$k):ncol(dat[[i]]$x)])
+ z <- .diagMatrix(z)
+ names(y) <- dimnames(x) <- dimnames(z) <- NULL
+ data <- list(y=y, x=x, z=z)
+ dat2 <- getDat(y~x-1, ~z-1, data=data)
+ attr(dat2, "ModelType") <- "linear"
+ res0 <- .tetlin(dat2, 1, "2sls")
+ tet0 <- res0$par
+ fsRes <- res0$Res
+ #tet0 <- tsls(y~x-1, ~z-1, data=data)$coefficients
+ tet0 <- rep(list(tet0), length(dat))
+ ## Getting the weights
+ w <- .weightFct_Sys(tet=tet0, dat=dat, type=object$vcov)
+ res <- .tetlin(dat2, w)
+ res$par <- c(res$par)
+ par <- list()
+ for (i in 1:length(dat))
+ {
+ par[[i]] <- res$par
+ names(par[[i]]) <- object$namesCoef[[i]]
+ }
+ names(par) <- names(dat)
+ df <- ncol(z) - ncol(x)
+ k <- ncol(x)
+ q <- ncol(z)
+ n <- nrow(dat[[1]]$x)
+ df.residuals <- n - k
+ } else {
+ ## Getting first step estimate ##
+ res0 <- lapply(1:length(dat), function(i) tsls(object$formula$g[[i]], object$formula$h[[i]], data=object$data))
+ tet0 <- lapply(1:length(dat), function(i) res0[[i]]$coefficients)
+ fsRes <- lapply(1:length(dat), function(i) res0[[i]]$fsRes)
+ # Getting the weights
+ w <- .weightFct_Sys(tet=tet0, dat=dat, type=object$vcov)
+ y <- lapply(1:length(dat), function(i) dat[[i]]$x[,1])
+ y <- do.call(c, y)
+ x <- lapply(1:length(dat), function(i) dat[[i]]$x[,2:(dat[[i]]$k+1)])
+ x <- .diagMatrix(x)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 88
More information about the Gmm-commits
mailing list