[Gmm-commits] r89 - in pkg/gmm: . R data man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 18 21:25:53 CET 2015
Author: chaussep
Date: 2015-11-18 21:25:53 +0100 (Wed, 18 Nov 2015)
New Revision: 89
Added:
pkg/gmm/R/Methods.sysGmm.R
pkg/gmm/R/sysGmm.R
pkg/gmm/data/Growth.rda
pkg/gmm/man/growth.Rd
Modified:
pkg/gmm/NEWS
pkg/gmm/R/Methods.gmm.R
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/momentEstim.R
pkg/gmm/man/getDat.Rd
pkg/gmm/man/sysGmm.Rd
Log:
Working on sysGmm functions
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2015-11-11 21:44:07 UTC (rev 88)
+++ pkg/gmm/NEWS 2015-11-18 20:25:53 UTC (rev 89)
@@ -24,6 +24,7 @@
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
+o A panel of Macro data has been added.
Changes in version 1.5-2
Modified: pkg/gmm/R/Methods.gmm.R
===================================================================
--- pkg/gmm/R/Methods.gmm.R 2015-11-11 21:44:07 UTC (rev 88)
+++ pkg/gmm/R/Methods.gmm.R 2015-11-18 20:25:53 UTC (rev 89)
@@ -43,50 +43,7 @@
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))
@@ -198,64 +155,6 @@
invisible(x)
}
-
-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))
@@ -309,23 +208,6 @@
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
Added: pkg/gmm/R/Methods.sysGmm.R
===================================================================
--- pkg/gmm/R/Methods.sysGmm.R (rev 0)
+++ pkg/gmm/R/Methods.sysGmm.R 2015-11-18 20:25:53 UTC (rev 89)
@@ -0,0 +1,138 @@
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# A copy of the GNU General Public License is available at
+# http://www.r-project.org/Licenses/
+
+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
+ }
+
+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)
+ printCoefmat(x$coefficients[[i]], digits=digits, ...)
+ 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)
+ }
+
+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)
+ }
+
+
+
+
+
+
+
+
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2015-11-11 21:44:07 UTC (rev 88)
+++ pkg/gmm/R/getModel.R 2015-11-18 20:25:53 UTC (rev 89)
@@ -24,7 +24,12 @@
getModel.sysGmm <- function(object, ...)
{
- object$allArg <- c(object, list(...))
+ if (object$commonCoef & !is.null(object$crossEquConst))
+ {
+ object$commonCoef <- FALSE
+ warning("When crossEquConst is not NULL, commonCoef=TRUE is ignore and set to FALSE")
+ }
+ object$allArg <- c(object, ...)
object$formula <- list(g=object$g, h=object$h)
if (!is.list(object$g))
stop("g must be list of formulas")
@@ -55,7 +60,8 @@
}
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))
+ dat <- lapply(1:length(object$g), function(i) try(getDat(object$g[[i]], object$h[[i]], data = object$data,
+ error=!object$commonCoef), silent=TRUE))
chk <- sapply(1:length(dat), function(i) class(dat[[i]]) == "try-error")
chk <- which(chk)
mess <- vector()
@@ -102,13 +108,30 @@
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)
+ k2 <- do.call(c,k)
+ if (object$commonCoef | !is.null(object$crossEquConst))
{
- 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")
+ if (!all(k2==k2[1]))
+ stop("In a common coefficient model the number of regressors is the same in each equation")
+ if (object$commonCoef)
+ totK <- k2[1]
+ else
+ totK <- length(dat)*(k2[1]-length(object$crossEquConst)) + length(object$crossEquConst)
+ if (sum(do.call(c,q))<totK)
+ stop("The number of moment conditions is less than the number of coefficients")
+ if (!is.null(object$crossEquConst))
+ {
+ object$crossEquConst <- sort(object$crossEquConst)
+ if (length(object$crossEquConst) == k2[1])
+ if (all(object$crossEquConst==(1:k2[1])))
+ {
+ object$crossEquConst <- NULL
+ object$commonCoef <- TRUE
+ }
+ }
}
- attr(object$x, "sysInfo") <- list(k=k, df=df, q=q, typeInst=typeInst, typeDesc=typeDesc, commonCoef=object$commonCoef)
+ attr(object$x, "sysInfo") <- list(k=k, df=df, q=q, typeInst=typeInst, typeDesc=typeDesc, commonCoef=object$commonCoef,
+ crossEquConst=object$crossEquConst)
object$g <- .momentFct_Sys
class(object) <- clname
return(object)
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2015-11-11 21:44:07 UTC (rev 88)
+++ pkg/gmm/R/gmm.R 2015-11-18 20:25:53 UTC (rev 89)
@@ -46,84 +46,6 @@
z
}
-five <- function(g, h, commonCoef = FALSE, data = NULL)
- {
- 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)
- 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, wmatrix = wmatrix, vcov = vcov, kernel = kernel,
- crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method, approx = approx,
- 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)
- z <- FinRes(z, Model_info)
- return(z)
- }
-
-
evalGmm <- function(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,
@@ -191,7 +113,7 @@
w
}
-getDat <- function (formula, h, data)
+getDat <- function (formula, h, data, error=TRUE)
{
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
@@ -211,6 +133,7 @@
if (inherits(h,'formula'))
{
tmp <- as.character(formula)
+ termsh <- terms(h)
h <- paste(tmp[2], "~", as.character(h)[2], sep="")
h <- as.formula(h)
mfh <- match.call(expand.dots = FALSE)
@@ -228,23 +151,61 @@
else
{
h <- as.matrix(h)
- h <- cbind(1,h)
+ chkInt <- sapply(1:ncol(h), function(i) all(h[,i]/mean(h[,i]) == 1))
+ if (sum(chkInt) > 1)
+ stop("Too many intercept in the matrix h")
+ if (any(chkInt))
+ {
+ h <- h[,!chkInt, drop=FALSE]
+ h <- cbind(1,h)
+ intercept_h <- TRUE
+ } else {
+ if (attr(mt,"intercept")==1)
+ {
+ h <- cbind(1, h)
+ intercept_h <- TRUE
+ } else {
+ intercept_h <- FALSE
+ }
+ }
if(is.null(colnames(h)))
- colnames(h) <- c("h.(Intercept)",paste("h",1:(ncol(h)-1),sep=""))
- else
- attr(h,'dimnames')[[2]][1] <- "h.(Intercept)"
- if (attr(mt,"intercept")==0)
{
- h <- as.matrix(h[,2:ncol(h)])
+ if (intercept_h)
+ colnames(h) <- c("(Intercept)",paste("h",1:(ncol(h)-1),sep=""))
+ else
+ colnames(h) <- paste("h",1:ncol(h),sep="")
+ } else {
+ if (intercept_h)
+ colnames(h)[1] <- "(Intercept)"
+ coln_h <- colnames(h)
+ coln_h <- gsub(" ", "", coln_h)
+ chk <- which(coln_h == "")
+ if (length(chk) >0)
+ coln_h[chk] <- paste("h", 1:length(chk), sep="")
+ if (any(duplicated(coln_h)))
+ stop("colnames of the matrix h must be unique")
+ colnames(h) <- coln_h
}
+ if (!intercept_h)
+ {
+ hFormula <- paste(colnames(h), collapse="+")
+ hFormula <- as.formula(paste("~", hFormula, "-1", sep=""))
+ } else {
+ hFormula <- paste(colnames(h)[-1], collapse="+")
+ hFormula <- as.formula(paste("~", hFormula, sep=""))
+ }
+ termsh <- terms(hFormula)
}
ny <- ncol(y)
k <- ncol(xt)
nh <- ncol(h)
if (nh<k)
- stop("The number of moment conditions must be at least equal to the number of coefficients to estimate")
+ {
+ if (error)
+ stop("The number of moment conditions must be at least equal to the number of coefficients to estimate")
+ }
if (is.null(colnames(y)))
- colnames(y) <- namey
+ colnames(y) <- namey
rownames(xt) <- rownames(y)
rownames(h) <- rownames(y)
x <- cbind(y,xt,h)
@@ -253,10 +214,13 @@
warning("There are missing values. Associated observations have been removed")
x <- na.omit(x)
if (nrow(x)<=k)
- stop("The number of observations must be greater than the number of coefficients")
+ {
+ if (error)
+ stop("The number of observations must be greater than the number of coefficients")
+ }
}
colnames(x)<-c(colnames(y),colnames(xt),colnames(h))
- return(list(x=x,nh=nh,ny=ny,k=k,mf=mf,mt=mt,cl=cl))
+ return(list(x=x,nh=nh,ny=ny,k=k,mf=mf,mt=mt,cl=cl,termsh=termsh,termsx=mt))
}
.tetlin <- function(dat, w, type=NULL)
@@ -284,14 +248,17 @@
if (!is.null(type))
{
if(type=="2sls")
- {
+ {
if (length(includeExo) > 0)
{
endo <- xm[, -includeExo, drop = FALSE]
endoName <- colnames(endo)
if (ncol(endo) != 0)
{
- restsls <- lm(endo~hm-1)
+ if (attr(dat$termsh, "intercept") == 1)
+ restsls <- lm(endo~hm[,-1])
+ else
+ restsls <- lm(endo~hm-1)
fsls <- xm
fsls[, -includeExo] <- restsls$fitted
} else {
@@ -299,7 +266,10 @@
restsls <- NULL
}
} else {
- restsls <- lm(xm~hm-1)
+ if (attr(dat$termsh, "intercept") == 1)
+ restsls <- lm(xm~hm[,-1])
+ else
+ restsls <- lm(xm~hm-1)
fsls <- restsls$fitted
endoName <- colnames(xm)
}
@@ -363,8 +333,9 @@
res$firstStageReg <- restsls
if (!is.null(restsls))
{
- res$fsRes <- summary(restsls)
- attr(res$fsRes, "Endo") <- endoName
+ chk <- .chkPerfectFit(restsls)
+ res$fsRes <- suppressWarnings(summary(restsls))[!chk]
+ attr(res$fsRes, "Endo") <- endoName[!chk]
}
}
return(res)
@@ -411,22 +382,6 @@
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)
{
if (!is.null(attr(dat, "eqConst")))
@@ -462,25 +417,6 @@
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")))
@@ -565,61 +501,6 @@
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")))
@@ -638,24 +519,4 @@
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-11 21:44:07 UTC (rev 88)
+++ pkg/gmm/R/momentEstim.R 2015-11-18 20:25:53 UTC (rev 89)
@@ -19,77 +19,39 @@
momentEstim.sysGmm.twoStep.formula <- function(object, ...)
{
- dat <- object$x
+ dat <- object$x
+ y <- lapply(1:length(dat), function(i) dat[[i]]$x[,1])
+ y <- do.call(c, y)
+ z <- lapply(1:length(dat), function(i)
+ dat[[i]]$x[,(2+dat[[i]]$k):ncol(dat[[i]]$x)])
+ z <- .diagMatrix(z)
+ x <- lapply(1:length(dat), function(i) dat[[i]]$x[,2:(dat[[i]]$k+1)])
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 if (!is.null(attr(dat, "sysInfo")$crossEquConst)) {
+ k <- attr(dat, "k")[[1]]
+ x <- .diagMatrix(x, (1:k)[-attr(dat, "sysInfo")$crossEquConst])
} 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)
- 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"
- res <- .tetlin(dat2, w)
- par <- list()
- k <- attr(dat, "sysInfo")$k
- for (i in 1:length(dat))
- {
- par[[i]] <- res$par[1:k[[i]]]
- names(par[[i]]) <- object$namesCoef[[i]]
- res$par <- res$par[-(1:k[[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
}
- z = list(coefficients = par, objective = res$value, dat=dat, k=k, q=q, df=df, df.residual=df.residual, n=n)
+ names(y) <- rownames(x) <- rownames(z) <- 1:length(y)
+ 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 <- .getThetaList(res0$par, dat)
+ fsRes <- res0$fsRes
+ w <- .weightFct_Sys(tet=tet0, dat=dat, type=object$vcov)
+ res <- .tetlin(dat2, w)
+ par <- .getThetaList(res$par, dat)
+ names(par) <- names(dat)
+ df <- ncol(z) - ncol(x)
+ k <- ncol(x)
+ q <- ncol(z)
+ n <- nrow(x)
+ df.residuals <- n - k
+ z = list(coefficients = par, objective = res$value, dat=dat, k=k, q=q, df=df, df.residual=df.residual, n=n)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 89
More information about the Gmm-commits
mailing list