[Gmm-commits] r52 - in pkg/gmm: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 16 21:26:34 CEST 2012
Author: chaussep
Date: 2012-04-16 21:26:33 +0200 (Mon, 16 Apr 2012)
New Revision: 52
Added:
pkg/gmm/R/gmmTests.R
pkg/gmm/man/KTest.Rd
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NAMESPACE
pkg/gmm/NEWS
pkg/gmm/R/getModel.R
pkg/gmm/R/momentEstim.R
Log:
Added test of Kleibergen (2005)
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2012-04-13 17:37:25 UTC (rev 51)
+++ pkg/gmm/DESCRIPTION 2012-04-16 19:26:33 UTC (rev 52)
@@ -1,6 +1,6 @@
Package: gmm
Version: 1.4-0
-Date: 2012-04-12
+Date: 2012-04-16
Title: Generalized Method of Moments and Generalized Empirical
Likelihood
Author: Pierre Chausse <pchausse at uwaterloo.ca>
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2012-04-13 17:37:25 UTC (rev 51)
+++ pkg/gmm/NAMESPACE 2012-04-16 19:26:33 UTC (rev 52)
@@ -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)
+ momentEstim.baseGel.modFormula,tsls,summary.tsls, print.summary.tsls, KTest, print.gmmTests)
S3method(summary, gmm)
S3method(summary, tsls)
@@ -25,6 +25,7 @@
S3method(print, gel)
S3method(print, summary.gel)
S3method(print, summary.tsls)
+S3method(print, gmmTests)
S3method(coef, gel)
S3method(vcov, gel)
S3method(confint, gel)
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2012-04-13 17:37:25 UTC (rev 51)
+++ pkg/gmm/NEWS 2012-04-16 19:26:33 UTC (rev 52)
@@ -28,6 +28,7 @@
o Cleaned the codes. The data are in object$dat and we can get the moment matrix by calling gt <- object$g(object$coef,object$dat) for linear and
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.
Changes in version 1.3-8
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2012-04-13 17:37:25 UTC (rev 51)
+++ pkg/gmm/R/getModel.R 2012-04-16 19:26:33 UTC (rev 52)
@@ -59,7 +59,7 @@
obj$g <- function(tet, dat)
{
resTet <- attr(dat,"eqConst")$eqConst
- tet2 <- vector(length=length(tet)+length(resTet))
+ tet2 <- vector(length=length(tet)+nrow(resTet))
tet2[resTet[,1]] <- resTet[,2]
tet2[-resTet[,1]] <- tet
attr(dat,"eqConst")$unConstg(tet2, dat)
@@ -67,7 +67,7 @@
obj$gradv <- function(tet, dat)
{
resTet <- attr(dat,"eqConst")$eqConst
- tet2 <- vector(length=length(tet)+length(resTet))
+ tet2 <- vector(length=length(tet)+nrow(resTet))
tet2[resTet[,1]] <- resTet[,2]
tet2[-resTet[,1]] <- tet
attr(dat,"eqConst")$unConstgradv(tet2, dat)[,-resTet[,1]]
@@ -185,7 +185,8 @@
object$iid<-iid
object$TypeGmm <- class(object)
object$gradv <- gradv
-
+ object$WSpec <- list(vcov = object$vcov, sandwich = list(kernel = object$kernel, bw = object$bw, prewhite = object$prewhite,
+ ar.method = object$ar.method, approx = object$approx, tol = object$tol))
class(object) <- clname
return(object)
}
Added: pkg/gmm/R/gmmTests.R
===================================================================
--- pkg/gmm/R/gmmTests.R (rev 0)
+++ pkg/gmm/R/gmmTests.R 2012-04-16 19:26:33 UTC (rev 52)
@@ -0,0 +1,166 @@
+# 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/
+
+
+# This function compute what is needed for the K statistics of KleiBergen (2005)
+#####################################################################################
+
+.BigCov <- function(obj,theta0)
+ {
+ insertRC <- function(A,w,v)
+ {
+ NewA <- matrix(ncol=ncol(A)+length(w),nrow=nrow(A)+length(w))
+ NewA[-w,-w] <- A
+ NewA[w,] <- v
+ NewA[,w] <- v
+ NewA
+ }
+ dg <- function(obj)
+ {
+ dat <- obj$dat
+ x <- dat$x
+ ny <- dat$ny
+ nh <- dat$nh
+ k <- dat$k
+ 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}
+
+ if (!is.list(obj$dat))
+ {
+ Myenv <- new.env()
+ assign("obj", obj, envir=Myenv)
+ assign("theta", theta0, envir=Myenv)
+ res <- numericDeriv(quote(g(theta,obj)),"theta",Myenv)
+ qT <- attr(res,"gradient")
+ } else {
+ qT <- dg(obj)}
+
+ qTmat <- apply(qT,3,colSums)
+ qT <- matrix(qT,nrow=dim(qT)[1])
+ gt <- obj$g(theta0,obj$dat)
+ fT <- colSums(gt)
+ n <- nrow(gt)
+ q <- ncol(gt)
+ All <- cbind(gt,qT)
+ All <- sweep(All,2,colMeans(All),FUN="-")
+ f <- function(x)
+ all(abs(x)<1e-7)
+ w <- which(apply(All,2,f))
+ if (length(w) != 0)
+ All <- All[,-w]
+ if (obj$WSpec$vcov == "iid")
+ {
+ V <- crossprod(All)/nrow(All)
+ } else {
+ class(All) <- "gmmFct"
+ argSand <- obj$WSpec$sandwich
+ argSand$x <- All
+ argSand$sandwich <- FALSE
+ V <- do.call(kernHAC,argSand)
+ }
+ if (length(w) != 0)
+ V <- insertRC(V,w,0)
+ Vff <- V[1:q,1:q]
+ Vthetaf <- V[(q+1):nrow(V),1:q]
+ list(Vff=Vff,Vthetaf=Vthetaf,qT=qTmat,fT=fT,n=n,q=q)
+ }
+
+KTest <- function(obj, theta0=NULL, alphaK = 0.04, alphaJ = 0.01)
+ {
+ if (class(obj) != "gmm")
+ stop("KTest is only for gmm type objects")
+
+ if (!is.null(attr(obj$dat,"eqConst")))
+ {
+ if (!is.null(theta0))
+ warning("setting a value for theta0 has no effect when the gmm is already constrained")
+ resTet <- attr(obj$dat,"eqConst")$eqConst
+ tet <- obj$coefficients
+ theta0 <- vector(length=length(tet)+nrow(resTet))
+ theta0[resTet[,1]] <- resTet[,2]
+ theta0[-resTet[,1]] <- tet
+ testName <- paste(rownames(resTet), " = ", resTet[,2], collapse="\n")
+ if (is.list(obj$dat))
+ {
+ x <- model.matrix(obj$dat$mt,obj$dat$mf)
+ y <- model.response(obj$dat$mf)
+ obj$dat$x <- cbind(y,x,obj$dat$x[,(obj$dat$ny+obj$dat$k+1):ncol(obj$dat$x)])
+ obj$dat$k <- ncol(x)
+ } else {
+ obj$g <- attr(obj$dat,"eqConst")$unConstg
+ }
+ dfK <- nrow(resTet)
+ } 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")
+ }
+ V <- .BigCov(obj, theta0)
+ Vff <- V$Vff
+ Vtf <- V$Vthetaf
+ qT <- V$qT
+ fT <- V$fT
+ dfJ <- V$q-length(theta0)
+ # the following is vec(D)
+ D <- c(qT)-Vtf%*%solve(Vff,fT)
+ D <- matrix(D,ncol=length(theta0))
+ 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
+ pvJ <- 1-pchisq(J,dfJ)
+ type <- c("K statistics","J statistics")
+ test <- c(K,J)
+ test <- cbind(test,c(pv,pvJ),c(dfK,dfJ))
+ dist <- paste("Chi_sq with ", c(dfK,dfJ), " degrees of freedom", sep="")
+ if(dfJ>0)
+ ans <- list(test=test,dist=dist,type=type,testName=testName)
+ else
+ ans <- list(test=matrix(test[1,],nrow=1),dist=dist[1],type=type[1],testName=testName)
+ if (pvJ<alphaJ) {
+ message <- "reject"
+ } else {
+ if (pv < alphaK)
+ message <- "reject"
+ else
+ message <- "do not reject"
+ }
+ ans$KJ <- (message == "do not reject")
+ ans$message <- paste("KJ-test result: We ", message, " H0 (alphaJ = ", alphaJ, ", alphaK = ", alphaK, ")", sep="")
+ class(ans) <- "gmmTests"
+ ans
+ }
+
+print.gmmTests <- function(x, digits = 5, ...)
+ {
+ lab <- paste("%0.",digits,"f",sep="")
+ cat("\nTest robust to weak identification\n")
+ cat("**********************************\n\n")
+ cat("The Null Hypothesis\n")
+ cat(x$testName,"\n\n")
+ for (i in 1:length(x$type))
+ {
+ cat(x$type[i],"\n")
+ cat("Test: ", sprintf(lab,x$test[i,1]), "(",x$dist[i],")\n")
+ cat("P-value: ", sprintf(lab,x$test[i,2]),"\n\n")
+ }
+ cat(x$message,"\n")
+ }
+
+
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2012-04-13 17:37:25 UTC (rev 51)
+++ pkg/gmm/R/momentEstim.R 2012-04-16 19:26:33 UTC (rev 52)
@@ -153,6 +153,7 @@
z$gradv <- P$gradv
z$iid <- P$iid
z$g <- P$g
+ z$WSpec <- P$WSpec
names(z$coefficients) <- P$namesCoef
if (is.null(colnames(z$gt)))
@@ -221,7 +222,8 @@
z$iid <- P$iid
z$g <- P$g
z$G <- P$gradv(dat)
-
+ z$WSpec <- P$WSpec
+
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -306,6 +308,7 @@
z$iid <- P$iid
z$g <- P$g
z$G <- P$gradv(dat)
+ z$WSpec <- P$WSpec
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -473,7 +476,8 @@
z$gradv <- P$gradv
z$iid <- P$iid
z$g <- P$g
-
+ z$WSpec <- P$WSpec
+
names(z$coefficients) <- P$namesCoef
if (is.null(colnames(z$gt)))
colnames(z$gt) <- paste("gt",1:ncol(z$gt),sep="")
@@ -573,7 +577,8 @@
z$G <- P$gradv(dat)
z$specMod <- P$specMod
z$cue <- list(weights=P$fixedKernW,message=P$weightMessage)
-
+ z$WSpec <- P$WSpec
+
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -653,6 +658,7 @@
names(z$coefficients) <- P$namesCoef
if (is.null(colnames(z$gt)))
colnames(z$gt) <- paste("gt",1:ncol(z$gt),sep="")
+ z$WSpec <- P$WSpec
z$specMod <- P$specMod
class(z) <- paste(P$TypeGmm, ".res", sep = "")
@@ -886,6 +892,7 @@
z$iid <- P$iid
z$g <- P$g
z$G <- P$gradv(dat)
+ z$WSpec <- P$WSpec
names(z$coefficients) <- P$namesCoef
colnames(z$gt) <- P$namesgt
@@ -978,6 +985,8 @@
z$gradv <- P$gradv
z$iid <- P$iid
z$g <- P$g
+ z$WSpec <- P$WSpec
+
names(z$coefficients) <- P$namesCoef
if (is.null(colnames(z$gt)))
colnames(z$gt) <- paste("gt",1:ncol(z$gt),sep="")
Added: pkg/gmm/man/KTest.Rd
===================================================================
--- pkg/gmm/man/KTest.Rd (rev 0)
+++ pkg/gmm/man/KTest.Rd 2012-04-16 19:26:33 UTC (rev 52)
@@ -0,0 +1,61 @@
+\name{KTest}
+\alias{KTest}
+\alias{print.gmmTests}
+
+\title{Compute the K statistics of Kleibergen}
+\description{The test is proposed by Kleibergen (2005). It is robust to weak identification.
+}
+\usage{
+KTest(obj, theta0 = NULL, alphaK = 0.04, alphaJ = 0.01)
+\method{print}{gmmTests}(x, digits = 5, ...)
+}
+\arguments{
+ \item{obj}{Object of class "gmm" returned by \code{\link{gmm}}}
+ \item{theta0}{The null hypothesis being tested. See details.}
+ \item{alphaK, alphaJ}{The size of the J and K tests when combining the two. The overall size is alphaK+alphaJ.}
+ \item{x}{An object of class \code{gmmTests} returned by \code{KTest}}
+ \item{digits}{The number of digits to be printed}
+ \item{...}{Other arguments when \code{print} is applied to another class object}
+}
+
+\details{
+The function produces the J-test and K-statistics which are robust to weak identification. The test is either \eqn{H0:\theta=theta_0}, in which case theta0 must be provided, or \eqn{\beta=\beta_0}, where \eqn{\theta=(\alpha', \beta')'}, and \eqn{\alpha} is assumed to be identified. In the latter case, theta0 is NULL and obj is a restricted estimation in which \eqn{\beta} is fixed to \eqn{\beta_0}. See \code{\link{gmm}} and the option "eqConst" for more details.
+}
+
+\value{
+Tests and p-values
+}
+
+\references{
+
+ Keibergen, F. (2005),
+ Testing Parameters in GMM without assuming that they are identified.
+ \emph{Econometrica}, \bold{73},
+ 1103-1123,
+
+}
+
+\examples{
+
+library(mvtnorm)
+sig <- matrix(c(1,.5,.5,1),2,2)
+n <- 400
+e <- rmvnorm(n,sigma=sig)
+x4 <- rnorm(n)
+w <- exp(-x4^2) + e[,1]
+y <- 0.1*w + e[,2]
+h <- cbind(x4, x4^2, x4^3, x4^6)
+g3 <- y~w
+res <- gmm(g3,h)
+
+# Testing the whole vector:
+
+KTest(res,theta0=c(0,.1))
+
+# Testing a subset of the vector (See \code{\link{gmm}})
+
+res2 <- gmm(g3, h, eqConst=matrix(c(2,.1),1,2))
+res2
+KTest(res2)
+
+}
More information about the Gmm-commits
mailing list