[Gmm-commits] r87 - in pkg/gmm: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 2 19:58:56 CET 2015
Author: chaussep
Date: 2015-11-02 19:58:56 +0100 (Mon, 02 Nov 2015)
New Revision: 87
Modified:
pkg/gmm/NEWS
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/momentEstim.R
Log:
Fixed a bug with nonlinear GEL with non-null eqConst
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2015-10-29 18:48:28 UTC (rev 86)
+++ pkg/gmm/NEWS 2015-11-02 18:58:56 UTC (rev 87)
@@ -19,6 +19,7 @@
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.
+o On the process to create a sysGmm for system og linear equations. The function sysGmm() will do the job. Not yet working.
Changes in version 1.5-2
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2015-10-29 18:48:28 UTC (rev 86)
+++ pkg/gmm/R/getModel.R 2015-11-02 18:58:56 UTC (rev 87)
@@ -16,6 +16,88 @@
UseMethod("getModel")
}
+getModel.sysGmm <- function(object, ...)
+ {
+ object$allArg <- c(object, list(...))
+ if (!is.list(object$g))
+ stop("g must be list of formulas")
+ if (length(object$g) == 1)
+ stop("For single equation GMM, use the function gmm()")
+ if (!all(sapply(1:length(object$g), function(i) is(object$g[[i]], "formula"))))
+ stop("g must be a list of formulas")
+ if (!is.list(object$h))
+ {
+ if(!is(object$h, "formula"))
+ stop("h is either a list of formulas or a formula")
+ else
+ object$h <- list(object$h)
+ } else {
+ if (!all(sapply(1:length(object$h), function(i) is(object$h[[i]], "formula"))))
+ stop("h is either a list of formulas or a formula")
+ }
+ 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))
+ typeDesc <- "System Gmm with common instruments"
+ } 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))
+ typeDesc <- "System Gmm"
+ }
+ 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]])
+ }
+ }
+ 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"
+ 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]))
+ 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)
+ 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)
+ object$TypeGmm <- class(object)
+ attr(object$x, "weight") <- list(w=object$weightsMatrix,
+ centeredVcov=object$centeredVcov)
+ attr(object$x, "weight")$WSpec <- list()
+ attr(object$x, "weight")$WSpec$sandwich <- list(kernel = object$kernel, bw = object$bw,
+ prewhite = object$prewhite,
+ ar.method = object$ar.method,
+ approx = object$approx, tol = object$tol)
+ attr(object$x, "weight")$vcov <- object$vcov
+ object$g <- .momentFct
+ class(object) <- clname
+ return(object)
+ }
+
+
getModel.constGmm <- function(object, ...)
{
class(object) <- "baseGmm"
@@ -63,7 +145,8 @@
getModel.baseGmm <- function(object, ...)
{
- object$allArg <- c(object, list(...))
+
+ object$allArg <- c(object, list(...))
if(is(object$g, "formula"))
{
object$gradv <- .DmomentFct
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2015-10-29 18:48:28 UTC (rev 86)
+++ pkg/gmm/R/gmm.R 2015-11-02 18:58:56 UTC (rev 87)
@@ -12,72 +12,107 @@
# 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"),
- 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,
- eqConstFullVcov = FALSE, ...)
-{
+ 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,
+ eqConstFullVcov = FALSE, ...)
+ {
+ type <- match.arg(type)
+ kernel <- match.arg(kernel)
+ vcov <- match.arg(vcov)
+ wmatrix <- match.arg(wmatrix)
+ optfct <- match.arg(optfct)
-type <- match.arg(type)
-kernel <- match.arg(kernel)
-vcov <- match.arg(vcov)
-wmatrix <- match.arg(wmatrix)
-optfct <- match.arg(optfct)
+ if (!is.null(eqConst))
+ TypeGmm <- "constGmm"
+
+ if(vcov=="TrueFixed" & is.null(weightsMatrix))
+ stop("TrueFixed vcov only for fixed weighting matrix")
+ if(!is.null(weightsMatrix))
+ wmatrix <- "optimal"
-if (!is.null(eqConst))
- TypeGmm <- "constGmm"
+ if(missing(data))
+ data<-NULL
+ all_args<-list(data = data, g = g, x = x, t0 = t0, gradv = gradv, type = type, 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)
+ class(all_args)<-TypeGmm
+ Model_info<-getModel(all_args, ...)
+ z <- momentEstim(Model_info, ...)
-if(vcov=="TrueFixed" & is.null(weightsMatrix))
- stop("TrueFixed vcov only for fixed weighting matrix")
-if(!is.null(weightsMatrix))
- wmatrix <- "optimal"
+ z <- FinRes(z, Model_info)
+ z
+ }
-if(missing(data))
- data<-NULL
-all_args<-list(data = data, g = g, x = x, t0 = t0, gradv = gradv, type = type, 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)
-class(all_args)<-TypeGmm
-Model_info<-getModel(all_args, ...)
-z <- momentEstim(Model_info, ...)
+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, ...)
+ {
+ type <- match.arg(type)
+ kernel <- match.arg(kernel)
+ vcov <- match.arg(vcov)
+ wmatrix <- match.arg(wmatrix)
+ optfct <- match.arg(optfct)
-z <- FinRes(z, Model_info)
-z
-}
+ 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,
+ 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)
+ class(all_args)<-TypeGmm
+ Model_info<-getModel(all_args, ...)
+ z <- momentEstim(Model_info, ...)
+ z <- FinRes(z, Model_info)
+ 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,
- model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, weightsMatrix = NULL, data)
-{
-TypeGmm = "baseGmm"
-type <- "eval"
-kernel <- match.arg(kernel)
-vcov <- match.arg(vcov)
-wmatrix <- match.arg(wmatrix)
-if (is.null(tetw) & is.null(weightsMatrix))
- stop("If the weighting matrix is not provided, you need to provide the vector of parameters tetw")
-
-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, x = x, t0 = t0, tetw = tetw, gradv = gradv, type = type, 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 = 100,
- optfct = NULL, model = model, X = X, Y = Y, call = match.call(), traceIter = NULL,
- eqConst = NULL, eqConstFullVcov = FALSE)
-class(all_args)<-TypeGmm
-Model_info<-getModel(all_args)
-class(Model_info) <- "baseGmm.eval"
-z <- momentEstim(Model_info)
-z <- FinRes(z, Model_info)
-z
-}
+ 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)
+ {
+ TypeGmm = "baseGmm"
+ type <- "eval"
+ kernel <- match.arg(kernel)
+ vcov <- match.arg(vcov)
+ wmatrix <- match.arg(wmatrix)
+ if (is.null(tetw) & is.null(weightsMatrix))
+ stop("If the weighting matrix is not provided, you need to provide the vector of parameters tetw")
+
+ 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, x = x, t0 = t0, tetw = tetw, gradv = gradv, type = type, 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 = 100,
+ optfct = NULL, model = model, X = X, Y = Y, call = match.call(), traceIter = NULL,
+ eqConst = NULL, eqConstFullVcov = FALSE)
+ class(all_args)<-TypeGmm
+ Model_info<-getModel(all_args)
+ class(Model_info) <- "baseGmm.eval"
+ z <- momentEstim(Model_info)
+ z <- FinRes(z, Model_info)
+ z
+ }
tsls <- function(g,x,data)
{
@@ -116,226 +151,226 @@
getDat <- function (formula, h, data)
{
- cl <- match.call()
- mf <- match.call(expand.dots = FALSE)
- m <- match(c("formula", "data"), names(mf), 0L)
- mf <- mf[c(1L, m)]
- mf$drop.unused.levels <- TRUE
- mf$na.action <- "na.pass"
- mf[[1L]] <- as.name("model.frame")
- mf <- eval(mf, parent.frame())
- mt <- attr(mf, "terms")
- 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 <- h
- mfh$h <- NULL
- mfh$drop.unused.levels <- TRUE
- mfh$na.action <- "na.pass"
- mfh[[1L]] <- as.name("model.frame")
- mfh <- eval(mfh, parent.frame())
- mth <- attr(mfh, "terms")
- h <- as.matrix(model.matrix(mth, mfh, NULL))
- }
- else
- {
- 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
- attr(h,'dimnames')[[2]][1] <- "h.(Intercept)"
- if (attr(mt,"intercept")==0)
- {
- h <- as.matrix(h[,2:ncol(h)])
- }
- }
- 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 (is.null(colnames(y)))
- {
- if (ny>1)
- colnames(y) <- paste("y",1:ncol(y),sep="")
- if (ny == 1)
- colnames(y) <- "y"
- }
- rownames(xt) <- rownames(y)
- rownames(h) <- rownames(y)
- x <- cbind(y,xt,h)
- if(any(is.na(x)))
- {
- 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")
- }
- 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))
+ cl <- match.call()
+ mf <- match.call(expand.dots = FALSE)
+ m <- match(c("formula", "data"), names(mf), 0L)
+ mf <- mf[c(1L, m)]
+ mf$drop.unused.levels <- TRUE
+ mf$na.action <- "na.pass"
+ mf[[1L]] <- as.name("model.frame")
+ mf <- eval(mf, parent.frame())
+ mt <- attr(mf, "terms")
+ 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 <- h
+ mfh$h <- NULL
+ mfh$drop.unused.levels <- TRUE
+ mfh$na.action <- "na.pass"
+ mfh[[1L]] <- as.name("model.frame")
+ mfh <- eval(mfh, parent.frame())
+ mth <- attr(mfh, "terms")
+ h <- as.matrix(model.matrix(mth, mfh, NULL))
+ }
+ else
+ {
+ 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
+ attr(h,'dimnames')[[2]][1] <- "h.(Intercept)"
+ if (attr(mt,"intercept")==0)
+ {
+ h <- as.matrix(h[,2:ncol(h)])
+ }
+ }
+ 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 (is.null(colnames(y)))
+ {
+ if (ny>1)
+ colnames(y) <- paste("y",1:ncol(y),sep="")
+ if (ny == 1)
+ colnames(y) <- "y"
+ }
+ rownames(xt) <- rownames(y)
+ rownames(h) <- rownames(y)
+ x <- cbind(y,xt,h)
+ if(any(is.na(x)))
+ {
+ 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")
+ }
+ 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))
}
.tetlin <- function(dat, w, type=NULL)
- {
- x <- dat$x
- g <- .momentFct
- gradv <- .DmomentFct
- ny <- dat$ny
- nh <- dat$nh
- k <- dat$k
- n <- nrow(x)
- ym <- as.matrix(x[,1:ny])
- xm <- as.matrix(x[,(ny+1):(ny+k)])
- hm <- as.matrix(x[,(ny+k+1):(ny+k+nh)])
- if (!is.null(attr(dat, "eqConst")))
- {
- resTet <- attr(dat,"eqConst")$eqConst
- y2 <- xm[, resTet[,1],drop=FALSE]%*%resTet[,2]
- ym <- ym-c(y2)
- xm <- xm[,-resTet[,1],drop=FALSE]
- k <- ncol(xm)
- }
- includeExo <- which(colnames(xm)%in%colnames(hm))
- inv <- attr(w, "inv")
- 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)
- fsls <- xm
- fsls[, -includeExo] <- restsls$fitted
- } else {
- fsls <- xm
- restsls <- NULL
- }
- } else {
- restsls <- lm(xm~hm-1)
- fsls <- restsls$fitted
- endoName <- colnames(xm)
- }
- par <- lm.fit(as.matrix(fsls), ym)$coefficients
- if (ny == 1)
- {
- e2sls <- ym-xm%*%par
- v2sls <- lm.fit(as.matrix(hm), e2sls)$fitted
- value <- sum(v2sls^2)/sum(e2sls^2)
- }
- else
- {
- par <- c(t(par))
- g2sls <- g(par, dat)
- w <- crossprod(g2sls)/n
- gb <- matrix(colMeans(g2sls), ncol = 1)
- value <- crossprod(gb, solve(w, gb))
- }
- }
- }
- 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
- {
- 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
- {
- 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)
- }
- 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(type))
- {
- if (type == "2sls")
- res$firstStageReg <- restsls
- if (!is.null(restsls))
- {
- res$fsRes <- summary(restsls)
- attr(res$fsRes, "Endo") <- endoName
- }
- }
- return(res)
- }
+ {
+ x <- dat$x
+ g <- .momentFct
+ gradv <- .DmomentFct
+ ny <- dat$ny
+ nh <- dat$nh
+ k <- dat$k
+ n <- nrow(x)
+ ym <- as.matrix(x[,1:ny])
+ xm <- as.matrix(x[,(ny+1):(ny+k)])
+ hm <- as.matrix(x[,(ny+k+1):(ny+k+nh)])
+ if (!is.null(attr(dat, "eqConst")))
+ {
+ resTet <- attr(dat,"eqConst")$eqConst
+ y2 <- xm[, resTet[,1],drop=FALSE]%*%resTet[,2]
+ ym <- ym-c(y2)
+ xm <- xm[,-resTet[,1],drop=FALSE]
+ k <- ncol(xm)
+ }
+ includeExo <- which(colnames(xm)%in%colnames(hm))
+ inv <- attr(w, "inv")
+ 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)
+ fsls <- xm
+ fsls[, -includeExo] <- restsls$fitted
+ } else {
+ fsls <- xm
+ restsls <- NULL
+ }
+ } else {
+ restsls <- lm(xm~hm-1)
+ fsls <- restsls$fitted
+ endoName <- colnames(xm)
+ }
+ par <- lm.fit(as.matrix(fsls), ym)$coefficients
+ if (ny == 1)
+ {
+ e2sls <- ym-xm%*%par
+ v2sls <- lm.fit(as.matrix(hm), e2sls)$fitted
+ value <- sum(v2sls^2)/sum(e2sls^2)
+ }
+ else
+ {
+ par <- c(t(par))
+ g2sls <- g(par, dat)
+ w <- crossprod(g2sls)/n
+ gb <- matrix(colMeans(g2sls), ncol = 1)
+ value <- crossprod(gb, solve(w, gb))
+ }
+ }
+ }
+ 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
+ {
+ 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
+ {
+ 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)
+ }
+ 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(type))
+ {
+ if (type == "2sls")
+ res$firstStageReg <- restsls
+ if (!is.null(restsls))
+ {
+ res$fsRes <- summary(restsls)
+ attr(res$fsRes, "Endo") <- endoName
+ }
+ }
+ return(res)
+ }
.obj1 <- function(thet, x, w)
{
- gt <- .momentFct(thet, x)
- gbar <- as.vector(colMeans(gt))
- INV <- attr(w, "inv")
- if (INV)
- obj <- crossprod(gbar, solve(w, gbar))
- else
- obj <- crossprod(gbar,w)%*%gbar
- return(obj)
- }
+ gt <- .momentFct(thet, x)
+ gbar <- as.vector(colMeans(gt))
+ INV <- attr(w, "inv")
+ if (INV)
+ obj <- crossprod(gbar, solve(w, gbar))
+ else
+ obj <- crossprod(gbar,w)%*%gbar
+ return(obj)
+ }
.Gf <- function(thet, x, pt = NULL)
- {
- myenv <- new.env()
- assign('x', x, envir = myenv)
- assign('thet', thet, envir = myenv)
- barg <- function(x, thet)
{
- gt <- .momentFct(thet, x)
- if (is.null(pt))
- gbar <- as.vector(colMeans(gt))
- else
- gbar <- as.vector(colSums(c(pt)*gt))
+ myenv <- new.env()
+ assign('x', x, envir = myenv)
+ assign('thet', thet, envir = myenv)
+ barg <- function(x, thet)
+ {
+ gt <- .momentFct(thet, x)
+ if (is.null(pt))
+ gbar <- as.vector(colMeans(gt))
+ else
+ gbar <- as.vector(colSums(c(pt)*gt))
- return(gbar)
+ return(gbar)
+ }
+ G <- attr(numericDeriv(quote(barg(x, thet)), "thet", myenv), "gradient")
+ return(G)
}
- G <- attr(numericDeriv(quote(barg(x, thet)), "thet", myenv), "gradient")
- return(G)
- }
.objCue <- function(thet, x, type=c("HAC", "iid", "ident", "fct", "fixed"))
{
- type <- match.arg(type)
- gt <- .momentFct(thet, x)
- gbar <- as.vector(colMeans(gt))
- w <- .weightFct(thet, x, type)
- obj <- crossprod(gbar,solve(w,gbar))
- return(obj)
-}
+ type <- match.arg(type)
+ gt <- .momentFct(thet, x)
+ gbar <- as.vector(colMeans(gt))
+ w <- .weightFct(thet, x, type)
+ obj <- crossprod(gbar,solve(w,gbar))
+ return(obj)
+ }
.momentFct <- function(tet, dat)
@@ -395,15 +430,17 @@
dgb <- dgb[,-attr(dat,"eqConst")$eqConst[,1], drop=FALSE]
} else {
if (is.null(attr(dat,"gradv")))
- dgb <- .Gf(tet2, dat, pt)
- else
- dgb <- attr(dat,"gradv")(tet2, dat)
- if (!is.null(attr(dat, "eqConst")))
- dgb <- dgb[,-attr(dat,"eqConst")$eqConst[,1], drop=FALSE]
+ {
+ dgb <- .Gf(tet, dat, pt)
+ } else {
+ dgb <- attr(dat,"gradv")(tet2, dat)
+ if (!is.null(attr(dat, "eqConst")))
+ dgb <- dgb[,-attr(dat,"eqConst")$eqConst[,1], drop=FALSE]
+ }
}
return(as.matrix(dgb))
}
-
+
.weightFct <- function(tet, dat, type=c("HAC", "iid", "ident", "fct", "fixed"))
{
type <- match.arg(type)
@@ -412,16 +449,16 @@
w <- attr(dat, "weight")$w
attr(w, "inv") <- FALSE
}
- else if (type == "ident")
- {
- w <- diag(attr(dat, "q"))
- attr(w, "inv") <- FALSE
- } else {
- gt <- .momentFct(tet,dat)
- if(attr(dat, "weight")$centeredVcov)
- gt <- residuals(lm(gt~1))
- n <- NROW(gt)
- }
+ else if (type == "ident")
+ {
+ w <- diag(attr(dat, "q"))
+ attr(w, "inv") <- FALSE
+ } else {
+ gt <- .momentFct(tet,dat)
+ if(attr(dat, "weight")$centeredVcov)
+ gt <- residuals(lm(gt~1))
+ n <- NROW(gt)
+ }
if (type == "HAC")
{
obj <- attr(dat, "weight")
@@ -460,5 +497,5 @@
e <- y-yhat
return(list(residuals=e, yhat=yhat))
}
-
-
+
+
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2015-10-29 18:48:28 UTC (rev 86)
+++ pkg/gmm/R/momentEstim.R 2015-11-02 18:58:56 UTC (rev 87)
@@ -783,8 +783,7 @@
G <- P$gradv(z$coefficients, x, z$pt)
z$G <- G
khat <- crossprod(c(z$pt)*z$gt, z$gt)/(P$k2)*P$bwVal
- G <- G/P$k1
-
+ G <- G/P$k1
kg <- solve(khat, G)
z$vcov_par <- try(solve(crossprod(G, kg))/n, silent=TRUE)
if (class(z$vcov_par) == "try-error")
More information about the Gmm-commits
mailing list