[Gmm-commits] r53 - in pkg/gmm: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 30 22:21:16 CEST 2012
Author: chaussep
Date: 2012-04-30 22:21:16 +0200 (Mon, 30 Apr 2012)
New Revision: 53
Added:
pkg/gmm/man/KConfid.Rd
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NAMESPACE
pkg/gmm/NEWS
pkg/gmm/R/FinRes.R
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/gmmTests.R
pkg/gmm/R/momentEstim.R
pkg/gmm/man/gmm.Rd
Log:
Keep adding functionalities see NEWS
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/DESCRIPTION 2012-04-30 20:21:16 UTC (rev 53)
@@ -12,7 +12,7 @@
methods that belong to the Generalized Empirical Likelihood
(GEL) family of estimators, as presented by Smith(1997),
Kitamura(1997), Newey-Smith(2004) and Anatolyev(2005).
-Depends: R (>= 2.10.0), sandwich
+Depends: R (>= 2.10.0), sandwich, multicore
Suggests: mvtnorm, car, fBasics, MASS, timeDate, timeSeries
Imports: stats
License: GPL (>= 2)
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/NAMESPACE 2012-04-30 20:21:16 UTC (rev 53)
@@ -8,7 +8,7 @@
specTest.gmm, specTest.gel, print.specTest, momentEstim.baseGmm.twoStep, momentEstim.baseGmm.twoStep.formula,
momentEstim.baseGmm.iterative.formula, momentEstim.baseGmm.iterative, momentEstim.baseGmm.cue.formula,
momentEstim.baseGmm.cue, getModel.baseGmm, getModel.baseGel, getModel.constGmm, FinRes.baseGmm.res, momentEstim.baseGel.mod,
- momentEstim.baseGel.modFormula,tsls,summary.tsls, print.summary.tsls, KTest, print.gmmTests)
+ momentEstim.baseGel.modFormula,tsls,summary.tsls, print.summary.tsls, KTest, print.gmmTests, KConfid, gmmWithConst)
S3method(summary, gmm)
S3method(summary, tsls)
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/NEWS 2012-04-30 20:21:16 UTC (rev 53)
@@ -29,6 +29,11 @@
non-linear models, where object is of class gmm.
o We can now impose equality constraints of the type theta_i=c_i with the option eqConst=. See help(gmm) which includes examples.
o The K statistics for weakly identified models of Kleiberben (2005) has been added. See ?KTest.
+o (gmm only) The bandwidth for the HAC estimation is set by the first step estimate. Its value, along with the weights,
+ is saved in attr(obj$w0,"Spec"). If the user is not satisfied with this automatic bw selection, bw can be set to any fixed number when
+ calling gmm(). See ?gmm for more details and examples.
+o The function gmmWithConst() reestimate an unrestricted model by imposing linear constraints, using the same specification. It also use the same bandwidth for
+ the estimation os the HAC matrix
Changes in version 1.3-8
Modified: pkg/gmm/R/FinRes.R
===================================================================
--- pkg/gmm/R/FinRes.R 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/R/FinRes.R 2012-04-30 20:21:16 UTC (rev 53)
@@ -34,15 +34,9 @@
}
else if(P$vcov == "HAC")
{
- if(P$centeredVcov)
- gmat <- lm(z$gt~1)
- else
- {
- gmat <- z$gt
- class(gmat) <- "gmmFct"
- }
- v <- kernHAC(gmat, kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
- ar.method = P$ar.method, approx = P$approx, tol = P$tol, sandwich = FALSE)
+ if (!is.null(attr(z$w0,"Spec")))
+ object$WSpec$sandwich$bw <- attr(z$w0,"Spec")$bw
+ v <- .myKernHAC(z$gt, object)
z$v <- v
}
@@ -100,6 +94,7 @@
z$weightsMatrix <- P$weightsMatrix
z$infVcov <- P$vcov
z$infWmatrix <- P$wmatrix
+ z$allArg <- P$allArg
z$met <- P$type
z$kernel <- P$kernel
z$coefficients <- c(z$coefficients)
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/R/getModel.R 2012-04-30 20:21:16 UTC (rev 53)
@@ -20,7 +20,6 @@
{
class(object) <- "baseGmm"
obj <- getModel(object)
-
if (!is.null(object$t0))
{
if (!is.null(dim(object$eqConst)))
@@ -47,13 +46,38 @@
{
if (obj$x$ny>1)
stop("Constrained GMM not implemented yet for system of equations")
+ obj$g2 <- function(tet, dat)
+ {
+ x <- attr(dat,"eqConst")$Xunc
+ y <- attr(dat,"eqConst")$Yunc
+ ny <- 1
+ nh <- dat$nh
+ tet <- matrix(tet, ncol = ncol(x))
+ e <- y - x %*% t(tet)
+ gt <- e * dat$x[, ny+dat$k+1]
+ if(nh > 1)
+ for (i in 2:nh) gt <- cbind(gt, e*x[, (ny+dat$k+i)])
+ return(gt)
+ }
+ obj$gradv2 <- function(dat)
+ {
+ x <- attr(dat,"eqConst")$Xunc
+ y <- attr(dat,"eqConst")$Yunc
+ nh <- dat$nh
+ k <- ncol(x)
+ ny <- 1
+ dat$x <- cbind(y,x,dat$x[, dat$ny+dat$k+(1:nh)])
+ dgb <- -(t(x[,(ny+k+1):(ny+k+nh)]) %*% x[,(ny+1):(ny+k)]) %x% diag(rep(1,ny))/nrow(x)
+ return(dgb)
+ }
+ attr(obj$x,"eqConst") <- list(unConstg = obj$g2, unConstgradv = obj$gradv2, eqConst = object$eqConst,
+ Yunc = obj$x$x[,1], Xunc = as.matrix(obj$x$x[,1+(1:obj$x$k)]))
x <- as.matrix(obj$x$x[,1+(object$eqConst[,1])])%*%object$eqConst[,2]
obj$x$x <- obj$x$x[,-(1+(object$eqConst[,1]))]
obj$x$x[,1] <- obj$x$x[,1]-x
obj$x$k <- obj$x$k-nrow(object$eqConst)
if (obj$x$k<=0)
stop("Nothing to estimate")
- attr(obj$x,"eqConst") <- list(eqConst = object$eqConst)
} else {
attr(obj$x,"eqConst") <- list(unConstg = obj$g, unConstgradv = obj$gradv, eqConst = object$eqConst)
obj$g <- function(tet, dat)
@@ -70,20 +94,24 @@
tet2 <- vector(length=length(tet)+nrow(resTet))
tet2[resTet[,1]] <- resTet[,2]
tet2[-resTet[,1]] <- tet
- attr(dat,"eqConst")$unConstgradv(tet2, dat)[,-resTet[,1]]
+ if (!is.null(as.list(args(attr(dat,"eqConst")$unConstgradv))$g))
+ attr(dat,"eqConst")$unConstgradv(tet2, dat, g=attr(dat,"eqConst")$unConstg)[,-resTet[,1]]
+ else
+ attr(dat,"eqConst")$unConstgradv(tet2, dat)[,-resTet[,1]]
}
}
obj$namesCoef <- obj$namesCoef[-object$eqConst[,1]]
obj$type <- paste(obj$type,"(with equality constraints)",sep=" ")
- clname <- strsplit(class(obj),"constGmm")
mess <- paste(rownames(object$eqConst), " = " , object$eqConst[,2], "\n",collapse="")
mess <- paste("#### Equality constraints ####\n",mess,"##############################\n\n",sep="")
- obj$specMod <- mess
+ obj$specMod <- mess
+
return(obj)
}
getModel.baseGmm <- function(object, ...)
{
+ object$allArg <- c(object, list(...))
if(is(object$g, "formula"))
{
object$gradvf <- FALSE
@@ -98,9 +126,11 @@
{
clname <- "baseGmm.twoStep.formula"
object$type <- "Linear model with iid errors: Regular IV or 2SLS"
- }
- else
+ }
+ else
+ {
clname <- paste(class(object), ".", object$type, ".formula", sep = "")
+ }
}
else
{
@@ -145,9 +175,11 @@
}
object$g <- g
object$x <- dat
+ attr(object$x,"ModelType") <- "linear"
}
else
{
+ attr(object$x,"ModelType") <- "nonlinear"
k <- length(object$t0)
if(is.null(names(object$t0)))
object$namesCoef <- paste("Theta[" ,1:k, "]", sep = "")
@@ -155,7 +187,9 @@
object$namesCoef <- names(object$t0)
if(is.null(object$weightsMatrix))
+ {
clname <- paste(class(object), "." ,object$type, sep = "")
+ }
else
{
clname <- "fixedW"
@@ -182,6 +216,7 @@
return(v)
}
+
object$iid<-iid
object$TypeGmm <- class(object)
object$gradv <- gradv
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/R/gmm.R 2012-04-30 20:21:16 UTC (rev 53)
@@ -54,6 +54,29 @@
return(ans)
}
+
+.myKernHAC <- function(gmat, obj)
+ {
+ if(obj$centeredVcov)
+ gmat <- lm(gmat~1)
+ else
+ class(gmat) <- "gmmFct"
+ AllArg <- obj$WSpec$sandwich
+ AllArg$x <- gmat
+ if (is.function(AllArg$bw))
+ {
+ bw <- AllArg$bw(gmat, order.by = AllArg$order.by, kernel = AllArg$kernel,
+ prewhite = AllArg$prewhite, ar.method = AllArg$ar.method)
+ AllArg$bw <- bw
+ }
+ weights <- do.call(weightsAndrews,AllArg)
+ AllArg$sandwich <- FALSE
+ AllArg$weights <- weights
+ w <- do.call(vcovHAC, AllArg)
+ attr(w,"Spec") <- list(weights = weights, bw = AllArg$bw, kernel = AllArg$kernel)
+ w
+ }
+
getDat <- function (formula, h, data)
{
cl <- match.call()
@@ -241,28 +264,12 @@
.objCue <- function(thet, x, P)
{
gt <- P$g(thet,x)
- if (is.null(P$fixedKernW))
- w <- weightsAndrews(lm(gt~1), prewhite=P$prewhite,
- bw = P$bw, kernel = P$kernel, approx = P$approx,
- ar.method = P$ar.method, tol = P$tol)
- else
- w <- P$fixedKernW
-
gbar <- as.vector(colMeans(gt))
if (P$vcov == "iid")
- w2 <- P$iid(thet, x, P$g, P$centeredVcov)
+ w <- P$iid(thet, x, P$g, P$centeredVcov)
if (P$vcov == "HAC")
- {
- if(P$centeredVcov)
- gmat <- lm(P$g(thet,x)~1)
- else
- {
- gmat <- P$g(thet,x)
- class(gmat) <- "gmmFct"
- }
- w2 <- vcovHAC(gmat, weights=w, sandwich = FALSE)
- }
- obj <- crossprod(gbar,solve(w2,gbar))
+ w <- .myKernHAC(gt, P)
+ obj <- crossprod(gbar,solve(w,gbar))
return(obj)
}
Modified: pkg/gmm/R/gmmTests.R
===================================================================
--- pkg/gmm/R/gmmTests.R 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/R/gmmTests.R 2012-04-30 20:21:16 UTC (rev 53)
@@ -28,21 +28,28 @@
dg <- function(obj)
{
dat <- obj$dat
- x <- dat$x
- ny <- dat$ny
- nh <- dat$nh
- k <- dat$k
+ if (!is.null(attr(dat,"eqConst")))
+ x <- attr(dat,"eqConst")$Xunc
+ else
+ x <- dat$x[,(dat$ny+1):(dat$ny+dat$k)]
+ k <- ncol(x)
+ h <- dat$x[,(dat$ny+dat$k+1):ncol(dat$x)]
qt <- array(dim=c(dim(obj$gt),k))
for (i in 1:k)
- qt[,,i] <- -x[, (ny + i)]*x[, (ny + k + 1):(ny+k+nh)]
+ qt[,,i] <- -x[,i]*h
qt}
- if (!is.list(obj$dat))
+ if (attr(obj$dat,"ModelType") == "nonlinear")
{
Myenv <- new.env()
assign("obj", obj, envir=Myenv)
assign("theta", theta0, envir=Myenv)
- res <- numericDeriv(quote(g(theta,obj)),"theta",Myenv)
+ gFunct <- if (!is.null(attr(obj$dat,"eqConst")))
+ attr(obj$dat,"eqConst")$unConstg
+ else
+ obj$g
+ assign("g",gFunct,envir=Myenv)
+ res <- numericDeriv(quote(g(theta,obj$dat)),"theta",Myenv)
qT <- attr(res,"gradient")
} else {
qT <- dg(obj)}
@@ -60,6 +67,8 @@
w <- which(apply(All,2,f))
if (length(w) != 0)
All <- All[,-w]
+ if (dim(All)[2] >= dim(All)[1])
+ stop("Too many moment conditions. Cannot estimate V")
if (obj$WSpec$vcov == "iid")
{
V <- crossprod(All)/nrow(All)
@@ -102,13 +111,15 @@
obj$g <- attr(obj$dat,"eqConst")$unConstg
}
dfK <- nrow(resTet)
- } else {
+ which <- resTet[,1]
+ } else {
if (is.null(theta0))
stop("You must either estimate a restricted model first or set theta0 under H0")
if (length(theta0) != length(obj$coef))
stop("theta0 is only for tests on the whole vector theta when obj is an unrestricted GMM")
dfK <- length(theta0)
testName <- paste(names(obj$coef), " = ", theta0, collapse="\n")
+ which <- 1:length(theta0)
}
V <- .BigCov(obj, theta0)
Vff <- V$Vff
@@ -119,8 +130,8 @@
# the following is vec(D)
D <- c(qT)-Vtf%*%solve(Vff,fT)
D <- matrix(D,ncol=length(theta0))
+ meat <- t(D)%*%solve(Vff,D)
bread <- t(fT)%*%solve(Vff,D)
- meat <- t(D)%*%solve(Vff,D)
K <- bread%*%solve(meat,t(bread))/V$n
pv <- 1-pchisq(K,dfK)
J <- t(fT)%*%solve(Vff,fT)/V$n-K
@@ -141,7 +152,8 @@
else
message <- "do not reject"
}
- ans$KJ <- (message == "do not reject")
+ ans$KJ <- (message == "do not reject")
+ ans$Matrix <- list(D=D,bread=bread,meat=meat,qT=qT,fT=fT)
ans$message <- paste("KJ-test result: We ", message, " H0 (alphaJ = ", alphaJ, ", alphaK = ", alphaK, ")", sep="")
class(ans) <- "gmmTests"
ans
@@ -163,4 +175,173 @@
cat(x$message,"\n")
}
+gmmWithConst <- function(obj, which, value)
+ {
+ argCall <- obj$allArg
+ argCall$call = NULL
+ if (!is.null(attr(obj$w0,"Spec")))
+ if (is.function(argCall$bw))
+ argCall$bw <- attr(obj$w0,"Spec")$bw
+ if (length(which)>=length(obj$coefficients))
+ stop("Too many constraints")
+ if (is.character(which))
+ {
+ if (any(!(which %in% names(obj$coefficients))))
+ stop("Wrong coefficient names in eqConst")
+ if (attr(obj$dat,"ModelType") == "linear")
+ which <- match(which,names(obj$coefficients))
+ }
+ if (!is.null(argCall$t0))
+ {
+ argCall$t0 <- obj$coefficients
+ argCall$t0[which] <- value
+ }
+
+ if (attr(obj$dat,"ModelType") == "nonlinear")
+ {
+ eqConst <- which
+ } else {
+ eqConst <- cbind(which,value)
+ }
+ argCall$eqConst <- eqConst
+ res <- do.call(gmm,argCall)
+ res$call <- match.call()
+ return(res)
+ }
+
+KConfid <- function(obj, which, type = c("K", "KJ"), alpha = 0.05, alphaJ = 0.01, n = 4)
+ {
+ type <- match.arg(type)
+ if ( (obj$df == 0) & (type == "KJ"))
+ stop("Only K type is available for just identified models")
+
+ if ( (alphaJ >= alpha) | (alphaJ <= 0) )
+ stop("We must have 0 < alphaJ < alpha")
+ if ( (alpha <= 0) | (alpha >=1) )
+ stop("We must have 0 < alpha < 1")
+ alphaK <- alpha-alphaJ
+
+ if (!is.null(attr(obj$dat,"eqConst")))
+ stop("Confidence intervals are constructed from unrestricted models")
+ if (length(which) > 2)
+ stop("This function computes confidence intervals for 1 or 2 coefficients only")
+ if (length(which)>length(obj$coefficients))
+ stop("length(which) must not exceed length(coefficients)")
+ if (is.character(which))
+ {
+ if (any(!(which %in% names(obj$coef))))
+ stop("names in which do not match names in coefficients")
+ else
+ which <- sort(which(which %in% names(obj$coef)))
+ } else {
+ if (any(which>length(obj$coefficients)) | any(which<=0))
+ stop("indices in which are not valid")
+ }
+ getUniInt <- function(tetx, obj, which, value=NULL, type , alpha, alphaJ, alphaK)
+ {
+ # value is the value of the other coefficient when we have 2 restrictions
+ # tetx is coef(obj)[which[1]] when we have two restrictions
+
+ if (is.null(value))
+ value <- tetx
+ else
+ value <- c(tetx,value)
+ res <- gmmWithConst(obj, which = which, value = value)
+ test <- KTest(res)
+ if (type == "K")
+ {
+ ifelse(alpha > test$test[1,2], 1, -1)
+ } else {
+ if (test$test[2,2]<alphaJ)
+ {
+ return(1)
+ } else {
+ ifelse(test$test[1,2]<alphaK, 1, -1)
+ }
+ }
+ }
+ gForMulti <- function(theta, obj, which, type , alpha, alphaJ, alphaK)
+ getUniInt(theta[1], obj, which, theta[2], type , alpha, alphaJ, alphaK)
+
+ if (length(which) == 1)
+ {
+ step <- 5*sqrt(diag(vcov(obj)))[which]
+ X0 <- obj$coef[which]
+ getBoth <- function(d)
+ {
+ res <- try(uniroot(getUniInt,c(X0,X0+d*step), obj = obj, which = which, type = type,
+ alpha = alpha, alphaJ = alphaJ, alphaK = alphaK, tol=1e-4), silent=TRUE)
+ if (class(res) == "try-error")
+ return(NA)
+ else
+ return(res$root)
+ }
+ res <- mclapply(c(-1,1), getBoth)
+ c(res[[1]],res[[2]])
+ } else {
+ step <- 5*sqrt(diag(vcov(obj)))[which]
+ sol <- .getCircle(obj$coef[which][1],obj$coef[which][2],gForMulti,n , step, obj=obj,
+ which=which, type=type, alpha=alpha, alphaJ=alphaJ, alphaK=alphaK)
+ colnames(sol) <- names(obj$coef)[which]
+ sol
+ }
+ }
+
+.getCircle <- function(x0,y0,g,n,b, trace=FALSE, ...)
+ {
+ tol=1e-4
+ if (any(b<=0))
+ stop("b must be strictly positive")
+
+
+ lambda <- seq(1,0,length=n)[-c(1,n)]
+ f <- function(x, y0, x0 = NULL, xi = NULL, yi = NULL)
+ {
+ if (is.null(xi))
+ y <- y0
+ else
+ y <- y0 + (yi-y0)/(xi-x0)*(x-x0)
+ g(c(x,y), ...)
+ }
+ f2 <- function(y, y0)
+ g(c(y0,y), ...)
+ f3 <- function(x)
+ ifelse(class(x)=="try-error",NA,x)
+
+ # get the four points of the cross
+ selectf <- list(f,f,f2,f2)
+ selectd <- c(1,-1,1,-1)
+ xy12 <- mclapply(1:4, function(i) try(uniroot(selectf[[i]],c(x0,x0+selectd[i]*b[1]),y0=y0,tol=tol)$root,silent=TRUE))
+
+ x1 <- xy12[[1]]
+ x2 <- xy12[[2]]
+ y1 <- xy12[[3]]
+ y2 <- xy12[[4]]
+
+ getAll <- function(lambda, dir=1, x)
+ {
+ yi <- y1*lambda + y0*(1-lambda)
+ xi <- x0*lambda + x*(1-lambda)
+ xt <- try(uniroot(f,c(x0,x0+dir*b[1]),y0=y0, x0=x0, xi=xi, yi=yi,tol=tol)$root,silent=TRUE)
+ xt <- f3(xt)
+ yt <- y0 + (yi-y0)/(xi-x0)*(xt-x0)
+ return(c(xt,yt))
+ }
+
+ res1 <- mclapply(lambda,getAll,dir=1,x=x1)
+ res2 <- mclapply(lambda,getAll,dir=-1,x=x2)
+ res1 <- t(simplify2array(res1))
+ res2 <- t(simplify2array(res2))
+ solU <- rbind(c(x2,y0),res2[length(lambda):1,],c(f3(x0), f3(y1)),res1,c(x1,y0))
+
+ res1 <- mclapply(lambda,getAll,dir=-1,x=x1)
+ res2 <- mclapply(lambda,getAll,dir=1,x=x2)
+ res1 <- t(simplify2array(res1))
+ res2 <- t(simplify2array(res2))
+ solD <- rbind(res2[length(lambda):1,],c(f3(x0), f3(y2)),res1)
+
+ sol <- rbind(solU,solD)
+ return(sol)
+ }
+
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/R/momentEstim.R 2012-04-30 20:21:16 UTC (rev 53)
@@ -85,16 +85,8 @@
if (P$vcov == "HAC")
{
- if(P$centeredVcov)
- gmat <- lm(P$g(res$par, P$x)~1)
- else
- {
- gmat <- P$g(res$par, P$x)
- class(gmat) <- "gmmFct"
- }
- w <- kernHAC(gmat, kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
- ar.method = P$ar.method, approx = P$approx, tol = P$tol, sandwich = FALSE)
-
+ gmat <- P$g(res$par, P$x)
+ w <- .myKernHAC(gmat, P)
}
if (P$optfct == "optim")
@@ -143,10 +135,10 @@
z$algoInfo <- list(convergence = res2$convergence, counts = res2$evaluations, message = res2$message)
}
- if(P$gradvf)
- z$G <- P$gradv(z$coefficients, P$x)
- else
- z$G <- P$gradv(z$coefficients, P$x, g = P$g)
+ if (length(as.list(args(P$gradv))) == 3)
+ z$G <- P$gradv(z$coefficients, P$x)
+ else
+ z$G <- P$gradv(z$coefficients, P$x, g = P$g)
z$dat <- P$x
z$gt <- P$g(z$coefficients, P$x)
@@ -160,6 +152,7 @@
colnames(z$gt) <- paste("gt",1:ncol(z$gt),sep="")
class(z) <- paste(P$TypeGmm,".res",sep="")
z$specMod <- P$specMod
+ z$w0 <- w
return(z)
}
@@ -175,10 +168,10 @@
n <- nrow(x)
q <- dat$ny*dat$nh
df <- q-k*dat$ny
-
+ w <- diag(q)
+
if (q == k2 | P$wmatrix == "ident")
{
- w <- diag(q)
res <- .tetlin(dat, w, P$gradv, P$g)
z = list(coefficients = res$par, objective = res$value, dat = dat, k = k, k2 = k2, n = n, q = q, df = df)
}
@@ -186,23 +179,16 @@
{
if (P$vcov == "iid")
{
- res2 <- .tetlin(dat, diag(q), P$gradv, P$g, type="2sls")
+ res2 <- .tetlin(dat, w, P$gradv, P$g, type="2sls")
initTheta <- NULL
}
if (P$vcov == "HAC")
{
- res1 <- .tetlin(dat, diag(q), P$gradv, P$g, type="2sls")
+ res1 <- .tetlin(dat, w, P$gradv, P$g, type="2sls")
initTheta <- res1$par
- if(P$centeredVcov)
- gmat <- lm(g(res1$par, dat)~1)
- else
- {
- gmat <- g(res1$par, dat)
- class(gmat) <- "gmmFct"
- }
- w <- kernHAC(gmat, kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
- ar.method = P$ar.method, approx = P$approx, tol = P$tol, sandwich = FALSE)
- res2 <- .tetlin(dat, w, P$gradv, g)
+ gmat <- g(res1$par, dat)
+ w <- .myKernHAC(gmat, P)
+ res2 <- .tetlin(dat, w, P$gradv, g)
}
z = list(coefficients = res2$par, objective = res2$value, dat=dat, k=k, k2=k2, n=n, q=q, df=df, initTheta = initTheta)
}
@@ -223,6 +209,7 @@
z$g <- P$g
z$G <- P$gradv(dat)
z$WSpec <- P$WSpec
+ z$w0 <- w
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -239,7 +226,6 @@
{
P <- object
g <- P$g
-
dat <- P$x
x <- dat$x
k <- dat$k
@@ -247,16 +233,15 @@
n <- nrow(x)
q <- dat$ny*dat$nh
df <- q-k*dat$ny
+ w <- diag(q)
if (q == k2 | P$wmatrix == "ident")
{
- w <- diag(q)
res <- .tetlin(dat, w, P$gradv, g)
z = list(coefficients = res$par, objective = res$value, dat = dat, k = k, k2 = k2, n = n, q = q, df = df)
}
else
{
- w <- diag(q)
res <- .tetlin(dat, w, P$gradv, g, type="2sls")
initTheta <- res$par
ch <- 100000
@@ -266,18 +251,10 @@
tet <- res$par
if (P$vcov == "HAC")
{
- if (P$centeredVcov)
- gmat <- lm(g(tet, dat)~1)
- else
- {
- gmat <- g(tet, dat)
- class(gmat) <- "gmmFct"
- }
- if (j==1)
- fixedKernW <- weightsAndrews(gmat, prewhite=P$prewhite,
- bw = P$bw, kernel = P$kernel, approx = P$approx,
- ar.method = P$ar.method, tol = P$tol)
- w <- vcovHAC(gmat, weights = fixedKernW, sandwich = FALSE)
+ gmat <- g(tet, dat)
+ if (!is.null(attr(w,"Spec")))
+ P$WSpec$sandwich$bw <- attr(w,"Spec")$bw
+ w <- .myKernHAC(gmat, P)
}
res <- .tetlin(dat, w, P$gradv, g)
ch <- crossprod(abs(tet- res$par)/tet)^.5
@@ -309,6 +286,7 @@
z$g <- P$g
z$G <- P$gradv(dat)
z$WSpec <- P$WSpec
+ z$w0 <- w
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -394,19 +372,10 @@
w <- P$iid(tet, P$x, P$g, P$centeredVcov)
if (P$vcov == "HAC")
{
- if (P$centeredVcov)
- gmat <- lm(P$g(tet, P$x)~1)
- else
- {
- gmat <- P$g(tet, P$x)
- class(gmat) <- "gmmFct"
- }
- if (j==1)
- fixedKernW <- weightsAndrews(gmat, prewhite=P$prewhite,
- bw = P$bw, kernel = P$kernel, approx = P$approx,
- ar.method = P$ar.method, tol = P$tol)
- w <- vcovHAC(gmat, weights = fixedKernW, sandwich = FALSE)
-
+ gmat <- P$g(tet, P$x)
+ if (!is.null(attr(w,"Spec")))
+ P$WSpec$sandwich$bw <- attr(w,"Spec")$bw
+ w <- .myKernHAC(gmat, P)
}
if (P$optfct == "optim")
@@ -465,18 +434,18 @@
}
+ if (length(as.list(args(P$gradv))) == 3)
+ z$G <- P$gradv(z$coefficients, P$x)
+ else
+ z$G <- P$gradv(z$coefficients, P$x, g = P$g)
- if(P$gradvf)
- z$G <- P$gradv(z$coefficients, P$x)
- else
- z$G <- P$gradv(z$coefficients, P$x, g = P$g)
-
z$dat <- P$x
z$gt <- P$g(z$coefficients, P$x)
z$gradv <- P$gradv
z$iid <- P$iid
z$g <- P$g
z$WSpec <- P$WSpec
+ z$w0 <- w
names(z$coefficients) <- P$namesCoef
if (is.null(colnames(z$gt)))
@@ -489,7 +458,6 @@
momentEstim.baseGmm.cue.formula <- function(object, ...)
{
- fixedKernWeights <- TRUE # to be changed or included as an option in gmm() in future version
P <- object
g <- P$g
@@ -500,10 +468,11 @@
n <- nrow(x)
q <- dat$ny*dat$nh
df <- q-k*dat$ny
-
+ w <- diag(q)
+ P$weightMessage <- "Weights for kernel estimate of the covariance are fixed and based on the first step estimate of Theta"
+
if (q == k2 | P$wmatrix == "ident")
{
- w <- diag(q)
res <- .tetlin(dat, w, P$gradv, g)
z = list(coefficients = res$par, objective = res$value, dat = dat, k = k, k2 = k2, n = n, q = q, df = df)
P$weightMessage <- "No CUE needed because the model is just identified"
@@ -512,32 +481,18 @@
{
if (is.null(P$t0))
{
- P$t0 <- .tetlin(dat,diag(q), P$gradv, g, type="2sls")$par
+ P$t0 <- .tetlin(dat, w, P$gradv, g, type="2sls")$par
initTheta <- P$t0
- if (fixedKernWeights)
- P$weightMessage <- "Weights for kernel estimate of the covariance are fixed and based on the first step estimate of Theta"
- else
- P$weightMessage <- "Weights for kernel estimate are variable inside the objective"
}
else
{
initTheta <- P$t0
- if (fixedKernWeights)
- P$weightMessage <- "Weights for kernel estimate of the covariance are fixed and based on the initial value of Theta provided by the user"
- else
- P$weightMessage <- "Weights for kernel estimate are variable inside the objective"
}
- if (fixedKernWeights)
- {
- gt0 <- g(P$t0,dat)
- gt0 <- lm(gt0~1)
- P$fixedKernW <- weightsAndrews(gt0, prewhite=P$prewhite,
- bw = P$bw, kernel = P$kernel, approx = P$approx,
- ar.method = P$ar.method, tol = P$tol)
- }
+ gt0 <- g(P$t0,dat)
+ w <- .myKernHAC(gt0, P)
+ P$WSpec$sandwich$bw <- attr(w,"Spec")$bw
-
if (P$optfct == "optim")
res2 <- optim(P$t0,.objCue, x = dat, P = P, ...)
if (P$optfct == "nlminb")
@@ -578,7 +533,7 @@
z$specMod <- P$specMod
z$cue <- list(weights=P$fixedKernW,message=P$weightMessage)
z$WSpec <- P$WSpec
-
+ z$w0 <- w
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -598,6 +553,7 @@
initTheta <- res$coef
n <- nrow(res$gt)
q <- ncol(res$gt)
+ w <- diag(q)
if (P$optfct == "optimize")
k = 1
@@ -615,10 +571,8 @@
else
{
gt0 <- P$g(res$coef,x) #Should we compute the weigths with the initial value provided?
- gt0 <- lm(gt0~1)
- P$fixedKernW <- weightsAndrews(gt0, prewhite=P$prewhite,
- bw = P$bw, kernel = P$kernel, approx = P$approx,
- ar.method = P$ar.method, tol = P$tol)
+ w <- .myKernHAC(gt0, P)
+ P$WSpec$sandwich$bw <- attr(w,"Spec")$bw
P$weightMessage <- "Weights for kernel estimate of the covariance are fixed and based on the first step estimate of Theta"
if (P$optfct == "optim")
@@ -644,10 +598,10 @@
}
- if(P$gradvf)
- z$G <- P$gradv(z$coefficients, P$x)
+ if (length(as.list(args(P$gradv))) == 3)
+ z$G <- P$gradv(z$coefficients, P$x)
else
- z$G <- P$gradv(z$coefficients, P$x, g = P$g)
+ z$G <- P$gradv(z$coefficients, P$x, g = P$g)
z$dat <- P$x
z$gradv <- P$gradv
@@ -975,10 +929,10 @@
else if(P$optfct == "nlminb")
z$algoInfo <- list(convergence = res2$convergence, counts = res2$evaluations, message = res2$message)
- if(P$gradvf)
- z$G <- P$gradv(z$coefficients, P$x)
+ if (length(as.list(args(P$gradv))) == 3)
+ z$G <- P$gradv(z$coefficients, P$x)
else
- z$G <- P$gradv(z$coefficients, P$x, g = P$g)
+ z$G <- P$gradv(z$coefficients, P$x, g = P$g)
z$dat <- P$x
z$gt <- P$g(z$coefficients, P$x)
Added: pkg/gmm/man/KConfid.Rd
===================================================================
--- pkg/gmm/man/KConfid.Rd (rev 0)
+++ pkg/gmm/man/KConfid.Rd 2012-04-30 20:21:16 UTC (rev 53)
@@ -0,0 +1,51 @@
+\name{KConfid}
+\alias{KConfid}
+
+\title{Confidence interval using the K statistics of Kleibergen}
+\description{The confidence is either an interval or an ellipse.
+}
+\usage{
+KConfid(obj, which, type = c("K", "KJ"), alpha = 0.05, alphaJ = 0.01, n = 4)
+}
+\arguments{
+ \item{obj}{Object of class "gmm" returned by \code{\link{gmm}} (not restricted)}
+ \item{type}{Should we base the confidence interval on the K or K-J statistics.}
+ \item{which}{A 2x1 vector or a scalar. The interval is computed for \code{coef(obj)[which]}. }
+ \item{alpha, alphaJ}{The overall size and the size for the J-test when type is "KS".}
+ \item{n}{The number of points to compute the confidence region is 4(n-1). It must be greater than 2.}
+}
+
+
+\value{
+Interval for \code{lenght(which)=1} or a series of points if \code{lenght(which)>1}.
+}
+
+\references{
+
+ Kleibergen, F. (2005),
+ Testing Parameters in GMM without assuming that they are identified.
+ \emph{Econometrica}, \bold{73},
+ 1103-1123,
+
+}
+
+\examples{
+
+data(Finance)
+r <- Finance[1:300, 1]
+rf <- Finance[1:300, "rf"]
+z <- as.matrix(r-rf)
+zm <- Finance[1:300, "rm"]-rf
+f1 <- zm
+f2 <- Finance[1:300, "hml"] - rf
+f3 <- Finance[1:300, "smb"] - rf
+res <- gmm(z ~ f1 + f2 + f3, ~ f1 + f2 + f3)
+
+KConfid(res,2)
+sol <- KConfid(res,c(2,3))
+plot(sol, main="Confidence Region")
+polygon(sol,col="grey")
+points(res$coef[2],res$coef[3],pch=21,bg=1)
+text(res$coef[2],res$coef[3],expression(hat(theta)),pos=3)
+
+}
Modified: pkg/gmm/man/gmm.Rd
===================================================================
--- pkg/gmm/man/gmm.Rd 2012-04-16 19:26:33 UTC (rev 52)
+++ pkg/gmm/man/gmm.Rd 2012-04-30 20:21:16 UTC (rev 53)
@@ -1,7 +1,8 @@
\name{gmm}
\alias{gmm}
-
+\alias{gmmWithConst}
+
\title{Generalized method of moment estimation}
\description{
@@ -15,6 +16,7 @@
prewhite = FALSE, ar.method = "ols", approx = "AR(1)", tol = 1e-7,
itermax = 100, optfct = c("optim","optimize","nlminb"), model=TRUE, X=FALSE, Y=FALSE,
TypeGmm = "baseGmm", centeredVcov = TRUE, weightsMatrix = NULL, traceIter = FALSE, data, eqConst = NULL, ...)
+gmmWithConst(obj, which, value)
}
\arguments{
\item{g}{A function of the form \eqn{g(\theta,x)} and which returns a \eqn{n \times q} matrix with typical element \eqn{g_i(\theta,x_t)} for \eqn{i=1,...q} and \eqn{t=1,...,n}. This matrix is then used to build the q sample moment conditions. It can also be a formula if the model is linear (see details below).}
@@ -63,6 +65,10 @@
\item{eqConst}{Either a named vector (if "g" is a function), a simple vector for the nonlinear case indicating which of the \eqn{\theta_0} is restricted, or a qx2 vector defining equality constraints of the form \eqn{\theta_i=c_i}. See below for an example.}
+\item{which, value}{The equality constraint is of the form which=value. "which" can be a vector of type characters with the names of the coefficients being constrained, or a vector of type numeric with the position of the coefficient in the whole vector.}
+
+\item{obj}{Object of class "gmm"}
+
\item{...}{More options to give to \code{\link{optim}}.}
}
@@ -96,6 +102,8 @@
\deqn{ 1~~~~~~~~~~ 0 }
\deqn{ -2\mu+2x ~~~~~ 2\sigma }
\deqn{-3\mu^{2}-3\sigma^{2} ~~~~ -6\mu\sigma}
+
+\code{gmmWithConst()} re-estimates an unrestricted model by adding an equality constraint.
}
\value{
@@ -291,6 +299,13 @@
res3 <- gmm(z ~ f1 + f2 + f3, ~ f1 + f2 + f3, t0=c(0,1,.5,.5), type="cue", eqConst = c(1,2))
res3
+### Examples with equality constraint using gmmWithConst
+###########################################################
+res2 <- gmm(z ~ f1 + f2 + f3, ~ f1 + f2 + f3)
+gmmWithConst(res2,c("f2","f3"),c(.5,.5))
+gmmWithConst(res2,c(2,3),c(.5,.5))
+
+
}
More information about the Gmm-commits
mailing list