[Gmm-commits] r11 - in pkg/gmm: . R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 17 06:48:11 CET 2009
Author: chaussep
Date: 2009-12-17 06:48:10 +0100 (Thu, 17 Dec 2009)
New Revision: 11
Added:
pkg/gmm/R/FinRes.R
pkg/gmm/R/getModel.R
pkg/gmm/R/momentEstim.R
pkg/gmm/man/FinRes.Rd
pkg/gmm/man/getModel.Rd
pkg/gmm/man/momentEstim.Rd
Modified:
pkg/gmm/DESCRIPTION
pkg/gmm/NAMESPACE
pkg/gmm/R/PlotMethods.R
pkg/gmm/R/gel.R
pkg/gmm/R/gmm.R
pkg/gmm/inst/doc/empir.bib
pkg/gmm/inst/doc/gmm_with_R.pdf
pkg/gmm/inst/doc/gmm_with_R.rnw
pkg/gmm/man/gmm.Rd
pkg/gmm/man/print.Rd
pkg/gmm/man/summary.Rd
Log:
Complete modification of the structure.
Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION 2009-12-03 06:30:24 UTC (rev 10)
+++ pkg/gmm/DESCRIPTION 2009-12-17 05:48:10 UTC (rev 11)
@@ -1,6 +1,6 @@
Package: gmm
-Version: 1.1-2
-Date: 2009-12-01
+Version: 1.1-3
+Date: 2009-12-15
Title: Generalized Method of Moments and Generalized Empirical Likelihood
Author: Pierre Chausse <pierre.chausse at uqam.ca>
Maintainer: Pierre Chausse <pierre.chausse at uqam.ca>
Modified: pkg/gmm/NAMESPACE
===================================================================
--- pkg/gmm/NAMESPACE 2009-12-03 06:30:24 UTC (rev 10)
+++ pkg/gmm/NAMESPACE 2009-12-17 05:48:10 UTC (rev 11)
@@ -2,7 +2,10 @@
export(HAC,gmm,weightsAndrews2,bwAndrews2,summary.gmm,rho,smoothG,getDat,bwNeweyWest2,summary.gel,getLamb,gel,
print.gmm,coef.gmm,vcov.gmm,print.summary.gmm, confint.gel, print.gel, print.summary.gel, vcov.gel, coef.gel, fitted.gmm,
- residuals.gmm, fitted.gel, residuals.gel, plot.gmm, plot.gel,formula.gmm, formula.gel, charStable)
+ residuals.gmm, fitted.gel, residuals.gel, plot.gmm, plot.gel,formula.gmm, formula.gel, charStable, specTest,
+ 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, FinRes.baseGmm.res)
S3method(summary, gmm)
S3method(summary, gel)
@@ -24,8 +27,16 @@
S3method(residuals, gel)
S3method(plot, gel)
S3method(formula, gel)
+S3method(specTest, gmm)
+S3method(specTest, gel)
+S3method(print, specTest)
+S3method(FinRes, baseGmm.res)
+S3method(getModel, baseGmm)
+S3method(momentEstim, baseGmm.twoStep)
+S3method(momentEstim, baseGmm.twoStep.formula)
+S3method(momentEstim, baseGmm.iterative.formula)
+S3method(momentEstim, baseGmm.iterative)
+S3method(momentEstim, baseGmm.cue.formula)
+S3method(momentEstim, baseGmm.cue)
-
-
-
Added: pkg/gmm/R/FinRes.R
===================================================================
--- pkg/gmm/R/FinRes.R (rev 0)
+++ pkg/gmm/R/FinRes.R 2009-12-17 05:48:10 UTC (rev 11)
@@ -0,0 +1,66 @@
+# 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/
+
+
+FinRes <- function(z, object, ...)
+ {
+# object is computed by the getModel method #
+ UseMethod("FinRes")
+ }
+
+FinRes.baseGmm.res <- function(z, object, ...)
+ {
+ P <- object
+ if(!is.null(object$gform))
+ {
+ dat <- z$dat
+ x <- dat$x
+ }
+ else
+ x <- z$x
+
+ n <- z$n
+ gradv <- z$gradv
+ iid <- z$iid
+
+ if(P$gradvf)
+ G <- gradv(z$coefficients, x)
+ else
+ G <- gradv(z$coefficients, x, g = object$g)
+
+ if (P$vcov == "iid")
+ v <- iid(z$coefficients, x, z$g)/n
+ else
+ {
+ v <- HAC(z$gt, kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
+ ar.method = P$ar.method, approx = P$approx, tol = P$tol)/n
+ }
+ if (P$wmatrix == "optimal")
+ {
+ z$vcov <- solve(crossprod(G, solve(v, G)))
+ }
+ else
+ {
+ GGG <- solve(crossprod(G), t(G))
+ z$vcov <- GGG %*% v %*% t(GGG)
+ }
+ dimnames(z$vcov) <- list(names(z$coefficients), names(z$coefficients))
+ z$call <- P$call
+ z$met <- P$type
+ z$kernel <- P$kernel
+ z$coefficients <- c(z$coefficients)
+ class(z) <- "gmm"
+ return(z)
+ }
+
+
Modified: pkg/gmm/R/PlotMethods.R
===================================================================
--- pkg/gmm/R/PlotMethods.R 2009-12-03 06:30:24 UTC (rev 10)
+++ pkg/gmm/R/PlotMethods.R 2009-12-17 05:48:10 UTC (rev 11)
@@ -53,7 +53,7 @@
ylim <- extendrange(r= ylim, f = 0.08)
plot(yh, r, xlab = "Fitted", ylab = "Residuals", main = main[1L],
ylim = ylim, type = "n", ...)
- panel(yh, r, ...)
+ panel(yh, r)
abline(h = 0, lty = 3, col = "gray")
}
if (show[2L]) { ## Normal
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2009-12-03 06:30:24 UTC (rev 10)
+++ pkg/gmm/R/gel.R 2009-12-17 05:48:10 UTC (rev 11)
@@ -344,22 +344,8 @@
names(z$badrho) <- "Number_of_bad_rho"
}
- Gf <- function(thet)
- {
- myenv <- new.env()
- assign('x',x,envir=myenv)
- assign('thet',thet,envir=myenv)
- barg <- function(x,thet)
- {
- gt <- g(thet,x)
- gbar <- as.vector(colMeans(gt))
- return(gbar)
- }
- G <- attr(numericDeriv(quote(barg(x,thet)),"thet",myenv),"gradient")
- return(G)
- }
if (!is.function(gradv))
- G <- Gf(z$coefficients)
+ G <- .Gf(z$coefficients, x, g)
else
G <- gradv(z$coefficients,x)
if (vcov == "iid")
Added: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R (rev 0)
+++ pkg/gmm/R/getModel.R 2009-12-17 05:48:10 UTC (rev 11)
@@ -0,0 +1,76 @@
+# 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/
+
+getModel <- function(object, ...)
+ {
+ UseMethod("getModel")
+ }
+
+getModel.baseGmm <- function(object, ...)
+ {
+ if(is(object$g, "formula"))
+ {
+ object$gradvf <- FALSE
+ dat <- getDat(object$g, object$x)
+ clname <- paste(class(object), ".", object$type, ".formula", sep = "")
+ object$gform<-object$g
+ g <- function(tet, x, ny = dat$ny, nh = dat$nh, k = dat$k)
+ {
+ tet <- matrix(tet, ncol = k)
+ e <- x[,1:ny] - x[,(ny+1):(ny+k)] %*% t(tet)
+ gt <- e * x[, ny+k+1]
+ if(nh > 1)
+ for (i in 2:nh) gt <- cbind(gt, e*x[, (ny+k+i)])
+ return(gt)
+ }
+ gradv <- function(tet, x, ny = dat$ny, nh = dat$nh, k = dat$k, g = NULL)
+ {
+ a <- g
+ tet <- NULL
+ dgb <- -(t(x[,(ny+k+1):(ny+k+nh)]) %*% x[,(ny+1):(ny+k)]) %x% diag(rep(1,ny))/nrow(x)
+ return(dgb)
+ }
+ object$g <- g
+ }
+ else
+ {
+ clname<-paste(class(object), "." ,object$type, sep = "")
+ if (!is.function(object$gradv))
+ {
+ gradv <- .Gf
+ object$gradvf <- FALSE
+ }
+ else
+ {
+ gradv <- object$gradv
+ object$gradvf <- TRUE
+ }
+ }
+
+ iid <- function(thet, x, g)
+ {
+ gt <- g(thet,x)
+ n <- ifelse(is.null(nrow(x)), length(x), nrow(x))
+ v <- crossprod(gt,gt)/n
+ return(v)
+ }
+
+ object$iid<-iid
+ object$TypeGmm <- class(object)
+ object$gradv <- gradv
+
+ class(object) <- clname
+ return(object)
+ }
+
+
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2009-12-03 06:30:24 UTC (rev 10)
+++ pkg/gmm/R/gmm.R 2009-12-17 05:48:10 UTC (rev 11)
@@ -13,369 +13,22 @@
gmm <- function(g,x,t0=NULL,gradv=NULL, type=c("twoStep","cue","iterative"), wmatrix = c("optimal","ident"), vcov=c("HAC","iid"),
kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen", "Tukey-Hanning"),crit=10e-7,bw = bwAndrews2,
- prewhite = FALSE, ar.method = "ols", approx="AR(1)",tol = 1e-7, itermax=100,optfct=c("optim","optimize","nlminb"), model=TRUE, X=FALSE, Y=FALSE, ...)
+ 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", ...)
{
+
type <- match.arg(type)
kernel <- match.arg(kernel)
vcov <- match.arg(vcov)
wmatrix <- match.arg(wmatrix)
optfct <- match.arg(optfct)
-typeg=0
-
- if (is(g,"formula"))
- {
- typeg=1
- dat <- getDat(g,x)
- x <- dat$x
- g <- function(tet,x,ny=dat$ny,nh=dat$nh,k=dat$k)
- {
- tet <- matrix(tet,ncol=k)
- e <- x[,1:ny] - x[,(ny+1):(ny+k)]%*%t(tet)
- gt <- e*x[,ny+k+1]
- if (nh > 1)
- {
- for (i in 2:nh)
- {
- gt <- cbind(gt,e*x[,(ny+k+i)])
- }
- }
- return(gt)
- }
- gradv <- function(x,ny=dat$ny,nh=dat$nh,k=dat$k)
- {
- dgb <- -(t(x[,(ny+k+1):(ny+k+nh)])%*%x[,(ny+1):(ny+k)])%x%diag(rep(1,ny))/nrow(x)
- return(dgb)
- }
- tetlin <- function(x,w,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)])
- whx <- solve(w,(crossprod(hm,xm)%x%diag(ny)))
- wvecyh <- solve(w,matrix(crossprod(ym,hm),ncol=1))
- dg <- gradv(x)
- xx <- crossprod(dg,whx)
- par <- solve(xx,crossprod(dg,wvecyh))
- gb <- matrix(colSums(g(par,x))/n,ncol=1)
- value <- crossprod(gb,solve(w,gb))
- res <- list(par=par,value=value)
- return(res)
- }
- }
-if (optfct == "optimize")
- {
- n = nrow(g(t0[1],x))
- q = ncol(g(t0[1],x))
- k = 1
- k2 <- k
- df <- q-k
- }
-else
- {
- if (typeg)
- {
- k <- dat$k
- k2 <- k*dat$ny
- n <- nrow(x)
- q <- dat$ny*dat$nh
- df <- q-k*dat$ny
- }
- else
- {
- n = nrow(g(t0,x))
- q = ncol(g(t0,x))
- k = length(t0)
- k2 <- k
- df <- q-k
- }
- }
-obj1 <- function(thet)
- {
- gt <- g(thet,x)
- gbar <- as.vector(colMeans(gt))
- obj <- crossprod(gbar,solve(w,gbar))
- return(obj)
- }
-iid <- function(thet)
- {
- gt <- g(thet,x)
- v <- crossprod(gt,gt)/n
- return(v)
- }
-Gf <- function(thet)
- {
- myenv <- new.env()
- assign('x',x,envir=myenv)
- assign('thet',thet,envir=myenv)
- barg <- function(x,thet)
- {
- gt <- g(thet,x)
- gbar <- as.vector(colMeans(gt))
- return(gbar)
- }
- G <- attr(numericDeriv(quote(barg(x,thet)),"thet",myenv),"gradient")
- return(G)
- }
-
-if (q == k2 | wmatrix == "ident")
- {
- if (typeg)
- {
- w <- diag(q)
- res <- tetlin(x,w)
- z = list(coefficients=res$par,objective=res$value)
- }
- else
- {
- w=diag(rep(1,q))
- if (optfct == "optim")
- res <- optim(t0,obj1, ...)
- if (optfct == "nlminb")
- {
- res <- nlminb(t0,obj1, ...)
- res$value <- res$objective
- }
- if (optfct == "optimize")
- {
- res <- optimize(obj1,t0, ...)
- res$par <- res$minimum
- res$value <- res$objective
- }
- z = list(coefficients=res$par,objective=res$value)
- }
- }
-else
- {
- if (type=="twoStep")
- {
- w=diag(rep(1,q))
- if (typeg)
- {
- res1 <- tetlin(x,w)
- }
- else
- {
- if (optfct == "optim")
- res1 <- optim(t0,obj1, ...)
- if (optfct == "nlminb")
- {
- res1 <- nlminb(t0,obj1, ...)
- res1$value <- res1$objective
- }
- if (optfct == "optimize")
- {
- res1 <- optimize(obj1,t0, ...)
- res1$par <- res1$minimum
- res1$value <- res1$objective
- }
- }
- if (vcov == "iid")
- w <- iid(res1$par)
- if (vcov == "HAC")
- w <- HAC(g(res1$par,x), kernel=kernel, bw=bw,prewhite=prewhite,ar.method=ar.method,approx=approx,tol=tol)
- if (typeg)
- {
- res2 <- tetlin(x,w)
- }
- else
- {
- if (optfct == "optim")
- res2 <- optim(res1$par,obj1, ...)
- if (optfct == "nlminb")
- {
- res2 <- nlminb(res1$par,obj1, ...)
- res2$value <- res2$objective
- }
- if (optfct == "optimize")
- {
- res2 <- optimize(obj1,t0, ...)
- res2$par <- res2$minimum
- res2$value <- res2$objective
- }
- }
- z = list(coefficients=res2$par,objective=res2$value)
- }
- if (type=="cue")
- {
- objCue <- function(thet)
- {
- gt <- g(thet,x)
- gbar <- as.vector(colMeans(gt))
- if (vcov == "iid")
- w2 <- iid(thet)
- if (vcov == "HAC")
- w2 <- HAC(g(thet,x), kernel=kernel, bw=bw,prewhite=prewhite,ar.method=ar.method,approx=approx,tol=tol)
- obj <- crossprod(gbar,solve(w2,gbar))
- return(obj)
- }
- if (typeg)
- {
- if (is.null(t0))
- t0 <- tetlin(x,diag(rep(1,q)))$par
- if (optfct == "optim")
- res2 <- optim(t0,objCue, ...)
- if (optfct == "nlminb")
- {
- res2 <- nlminb(t0,objCue, ...)
- res2$value <- res2$objective
- }
- if (optfct == "optimize")
- {
- res2 <- optimize(objCue,t0, ...)
- res2$par <- res2$minimum
- res2$value <- res2$objective
- }
- }
- else
- {
- if (optfct == "optim")
- res2 <- optim(t0,objCue, ...)
- if (optfct == "nlminb")
- {
- res2 <- nlminb(t0,objCue, ...)
- res2$value <- res2$objective
- }
- if (optfct == "optimize")
- {
- res2 <- optimize(objCue,t0, ...)
- res2$par <- res2$minimum
- res2$value <- res2$objective
- }
- }
- z = list(coefficients=res2$par,objective=res2$value)
- }
- if (type=="iterative")
- {
- w=diag(rep(1,q))
- if (typeg)
- res <- tetlin(x,w)
- else
- {
- if (optfct == "optim")
- res <- optim(t0,obj1, ...)
- if (optfct == "nlminb")
- {
- res <- nlminb(t0,obj1, ...)
- res$value <- res$objective
- }
- if (optfct == "optimize")
- {
- res <- optimize(obj1,t0, ...)
- res$par <- res$minimum
- res$value <- res$objective
- }
- }
- ch <- 100000
- j <- 1
- while(ch>crit)
- {
- tet <- res$par
- if (vcov == "iid")
- w <- iid(tet)
- if (vcov == "HAC")
- w <- HAC(g(tet,x), kernel=kernel, bw=bw,prewhite=prewhite,ar.method=ar.method,approx=approx,tol=tol)
- if (typeg)
- res <- tetlin(x,w)
- else
- {
- if (optfct == "optim")
- res <- optim(tet,obj1, ...)
- if (optfct == "nlminb")
- {
- res <- nlminb(tet,obj1, ...)
- res$value <- res$objective
- }
- if (optfct == "optimize")
- {
- res <- optimize(obj1,t0, ...)
- res$par <- res$minimum
- res$value <- res$objective
- }
- }
- ch <- crossprod(abs(tet-res$par)/tet,abs(tet-res$par)/tet)
- if (j>itermax)
- {
- cat("No convergence after ",itermax," iterations")
- ch <- crit
- }
- j <- j+1
- }
-
- z = list(coefficients=res$par,objective=res$value)
- }
- }
-
- if (!is.function(gradv))
- G <- Gf(z$coefficients)
- else
- if (typeg)
- G <- gradv(x)
- else
- G <- gradv(z$coefficients,x)
-
- if (vcov == "iid")
- v <- iid(z$coefficients)/n
- else
- v <- HAC(g(z$coefficients,x), kernel=kernel, bw=bw,prewhite=prewhite,ar.method=ar.method,approx=approx,tol=tol)/n
-
- if (wmatrix == "optimal")
- {
- z$vcov <- solve(crossprod(G,solve(v,G)))
- }
- else
- {
- GGG <- solve(crossprod(G),t(G))
- z$vcov <- GGG%*%v%*%t(GGG)
- }
-
-z$gt <- g(z$coefficients,x)
-if (typeg==0)
- names(z$coefficients) <- paste("Theta[",1:k,"]",sep="")
-if (typeg == 1)
- {
- namex <- colnames(dat$x[,(dat$ny+1):(dat$ny+dat$k)])
- nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
- if (dat$ny > 1)
- {
- namey <- colnames(dat$x[,1:dat$ny])
- names(z$coefficients) <- paste(rep(namey,dat$k),"_",rep(namex,rep(dat$ny,dat$k)),sep="")
- colnames(z$gt) <- paste(rep(namey,dat$nh),"_",rep(nameh,rep(dat$ny,dat$nh)),sep="")
- }
- if (dat$ny == 1)
- {
- names(z$coefficients) <- namex
- colnames(z$gt) <- nameh
- }
- }
-
-dimnames(z$vcov) <- list(names(z$coefficients),names(z$coefficients))
-z$df <- df
-z$k <- k
-z$n <- n
-if(typeg==1)
- {
- b <- z$coefficients
- y <- as.matrix(model.response(dat$mf, "numeric"))
- ny <- dat$ny
- b <- t(matrix(b,nrow=ny))
- x <- as.matrix(model.matrix(dat$mt, dat$mf, NULL))
- yhat <- x%*%b
- z$fitted.values<-yhat
- z$residuals<-y-yhat
- z$terms<- dat$mt
- if(model) z$model<-dat$mf
- if(X) z$x<-x
- if(Y) z$y<-y
- }
-if(typeg==0)
- if(X) z$x <- x
-
-z$call<-match.call()
-z$met <- type
-z$kernel <- kernel
-z$coefficients <- c(z$coefficients)
-class(z) <- "gmm"
+all_args<-list(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,
+ tol = tol, itermax = itermax, optfct = optfct, model = model, X = X, Y = Y, call = match.call())
+class(all_args)<-TypeGmm
+Model_info<-getModel(all_args)
+z <- momentEstim(Model_info, ...)
+z <- FinRes(z, Model_info)
z
}
@@ -420,3 +73,57 @@
}
+.tetlin <- function(x, w, ny, nh, k, gradv, g)
+ {
+ 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)])
+ whx <- solve(w, (crossprod(hm, xm) %x% diag(ny)))
+ wvecyh <- solve(w, matrix(crossprod(ym, hm), ncol = 1))
+ dg <- gradv(NULL,x, ny, nh, k)
+ xx <- crossprod(dg, whx)
+ par <- solve(xx, crossprod(dg, wvecyh))
+ gb <- matrix(colSums(g(par, x, ny, nh, k))/n, ncol = 1)
+ value <- crossprod(gb, solve(w, gb))
+ res <- list(par = par, value = value)
+ return(res)
+ }
+
+
+.obj1 <- function(thet, x, w, gf)
+ {
+ gt <- gf(thet, x)
+ gbar <- as.vector(colMeans(gt))
+ obj <- crossprod(gbar, solve(w, gbar))
+ return(obj)
+ }
+
+.Gf <- function(thet, x, g)
+ {
+ myenv <- new.env()
+ assign('x', x, envir = myenv)
+ assign('thet', thet, envir = myenv)
+ barg <- function(x, thet)
+ {
+ gt <- g(thet, x)
+ gbar <- as.vector(colMeans(gt))
+ return(gbar)
+ }
+ G <- attr(numericDeriv(quote(barg(x, thet)), "thet", myenv), "gradient")
+ return(G)
+ }
+
+.objCue <- function(thet, x, P)
+ {
+ gt <- P$g(thet,x)
+ gbar <- as.vector(colMeans(gt))
+ if (P$vcov == "iid")
+ w2 <- P$iid(thet, x, P$g)
+ if (P$vcov == "HAC")
+ w2 <- HAC(P$g(thet,x), kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
+ ar.method = P$ar.method, approx = P$approx, tol = P$tol)
+ obj <- crossprod(gbar,solve(w2,gbar))
+ return(obj)
+}
+
Added: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R (rev 0)
+++ pkg/gmm/R/momentEstim.R 2009-12-17 05:48:10 UTC (rev 11)
@@ -0,0 +1,464 @@
+# 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/
+
+momentEstim <- function(object, ...)
+ {
+ UseMethod("momentEstim")
+ }
+
+momentEstim.baseGmm.twoStep <- function(object, ...)
+ {
+ P <- object
+ x <- P$x
+ if (P$optfct == "optimize")
+ {
+ n = nrow(P$g(P$t0[1], x))
+ q = ncol(P$g(P$t0[1], x))
+ k = 1
+ }
+ else
+ {
+ n = nrow(P$g(P$t0, x))
+ q = ncol(P$g(P$t0, x))
+ k = length(P$t0)
+ }
+ k2 <- k
+ df <- q - k
+ w=diag(q)
+ if (P$optfct == "optim")
+ {
+ res <- optim(P$t0, .obj1, x = P$x, w = w, gf = P$g, ...)
+ }
+ if (P$optfct == "nlminb")
+ {
+ res <- nlminb(P$t0, .obj1, x = P$x, w = w, gf = P$g, ...)
+ res$value <- res$objective
+ }
+ if (P$optfct == "optimize")
+ {
+ res <- optimize(.obj1, P$t0, x = P$x, w = w, gf = P$g, ...)
+ res$par <- res$minimum
+ res$value <- res$objective
+ }
+ if (q == k2 | P$wmatrix == "ident")
+ z = list(coefficients = res$par, objective = res$value, k=k, k2=k2, n=n, q=q, df=df)
+ else
+ {
+ if (P$vcov == "iid")
+ w <- P$iid(res$par, P$x, P$g)
+
+ if (P$vcov == "HAC")
+ w <- HAC(P$g(res$par, P$x), kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
+ ar.method = P$ar.method, approx = P$approx, tol = P$tol)
+
+ if (P$optfct == "optim")
+ res2 <- optim(res$par, .obj1, x = P$x, w = w, gf = P$g, ...)
+
+ if (P$optfct == "nlminb")
+ {
+ res2 <- nlminb(res$par, .obj1, x = P$x, w = w, gf = P$g, ...)
+ res2$value <- res2$objective
+ }
+
+ if (P$optfct == "optimize")
+ {
+ res2 <- optimize(.obj1, P$t0, x = P$x, w = w, gf = P$g, ...)
+ res2$par <- res2$minimum
+ res2$value <- res2$objective
+ }
+
+ z = list(coefficients = res2$par, objective = res2$value, k=k, k2=k2, n=n, q=q, df=df)
+ }
+
+ z$x <- P$x
+ z$gt <- P$g(z$coefficients, P$x)
+ z$gradv <- P$gradv
+ z$iid <- P$iid
+ z$g <- P$g
+
+ class(z) <- paste(P$TypeGmm,".res",sep="")
+ return(z)
+ }
+
+momentEstim.baseGmm.twoStep.formula <- function(object, ...)
+ {
+ P <- object
+ g <- P$g
+ dat <- getDat(P$gform, P$x)
+ x <- dat$x
+ k <- dat$k
+ k2 <- k*dat$ny
+ n <- nrow(x)
+ q <- dat$ny*dat$nh
+ df <- q-k*dat$ny
+
+ if (q == k2 | P$wmatrix == "ident")
+ {
+ w <- diag(q)
+ res <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, P$g)
+ z = list(coefficients = res$par, objective = res$value, dat = dat, k = k, k2 = k2, n = n, q = q, df = df)
+ }
+ else
+ {
+ w=diag(rep(1, q))
+ res1 <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, P$g)
+ if (P$vcov == "iid")
+ w <- P$iid(res1$par, x, g)
+ if (P$vcov == "HAC")
+ w <- HAC(g(res1$par, x), kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
+ ar.method = P$ar.method, approx = P$approx, tol = P$tol)
+ res2 <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, g)
+ z = list(coefficients = res2$par, objective = res2$value, dat=dat, k=k, k2=k2, n=n, q=q, df=df)
+ }
+ z$gt <- g(z$coefficients, x)
+ b <- z$coefficients
+ y <- as.matrix(model.response(dat$mf, "numeric"))
+ ny <- dat$ny
+ b <- t(matrix(b, nrow = dat$ny))
+ x <- as.matrix(model.matrix(dat$mt, dat$mf, NULL))
+ yhat <- x %*% b
+ z$dat <- dat
+ z$fitted.values <- yhat
+ z$residuals <- y - yhat
+ z$terms <- dat$mt
+ if(P$model) z$model <- dat$mf
+ if(P$X) z$x <- x
+ if(P$Y) z$y <- y
+ z$gradv <- P$gradv
+ z$iid <- P$iid
+ z$g <- P$g
+
+ namex <- colnames(dat$x[,(dat$ny+1):(dat$ny+dat$k)])
+ nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
+
+ if (dat$ny > 1)
+ {
+ namey <- colnames(dat$x[,1:dat$ny])
+ names(z$coefficients) <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
+ colnames(z$gt) <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny, dat$nh)), sep = "")
+ }
+
+ if (dat$ny == 1)
+ {
+ names(z$coefficients) <- namex
+ colnames(z$gt) <- nameh
+ }
+ class(z) <- paste(P$TypeGmm,".res",sep="")
+ return(z)
+ }
+
+momentEstim.baseGmm.iterative.formula <- function(object, ...)
+ {
+ P <- object
+ g <- P$g
+ dat <- getDat(P$gform, P$x)
+ x <- dat$x
+ k <- dat$k
+ k2 <- k*dat$ny
+ n <- nrow(x)
+ q <- dat$ny*dat$nh
+ df <- q-k*dat$ny
+
+ if (q == k2 | P$wmatrix == "ident")
+ {
+ w <- diag(q)
+ res <- .tetlin(x, w, dat$ny, dat$nh, dat$k, 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(rep(1, q))
+ res <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, g)
+ ch <- 100000
+ j <- 1
+ while(ch > P$crit)
+ {
+ tet <- res$par
+ if (P$vcov == "iid")
+ w <- P$iid(tet, x, g)
+ if (P$vcov == "HAC")
+ w <- HAC(g(tet, x), kernel = P$kernel, bw = P$bw, prewhite = P$prewhite, ar.method = P$ar.method, approx = P$approx, tol = P$tol)
+ res <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, g)
+ ch <- crossprod(abs(tet- res$par)/tet)^.5
+ if (j>P$itermax)
+ {
+ cat("No convergence after ", P$itermax, " iterations")
+ ch <- P$crit
+ }
+ j <- j+1
+ }
+ z = list(coefficients = res$par, objective = res$value, dat=dat, k=k, k2=k2, n=n, q=q, df=df)
+ }
+ z$gt <- g(z$coefficients, x)
+ b <- z$coefficients
+ y <- as.matrix(model.response(dat$mf, "numeric"))
+ ny <- dat$ny
+ b <- t(matrix(b, nrow = dat$ny))
+ x <- as.matrix(model.matrix(dat$mt, dat$mf, NULL))
+ yhat <- x %*% b
+ z$dat <- dat
+ z$fitted.values <- yhat
+ z$residuals <- y - yhat
+ z$terms <- dat$mt
+ if(P$model) z$model <- dat$mf
+ if(P$X) z$x <- x
+ if(P$Y) z$y <- y
+ z$gradv <- P$gradv
+ z$iid <- P$iid
+ z$g <- P$g
+
+ namex <- colnames(dat$x[,(dat$ny+1):(dat$ny+dat$k)])
+ nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
+
+ if (dat$ny > 1)
+ {
+ namey <- colnames(dat$x[,1:dat$ny])
+ names(z$coefficients) <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
+ colnames(z$gt) <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny, dat$nh)), sep = "")
+ }
+
+ if (dat$ny == 1)
+ {
+ names(z$coefficients) <- namex
+ colnames(z$gt) <- nameh
+ }
+ class(z) <- paste(P$TypeGmm,".res",sep="")
+ return(z)
+ }
+
+momentEstim.baseGmm.iterative <- function(object, ...)
+ {
+ P <- object
+ x <- P$x
+ if (P$optfct == "optimize")
+ {
+ n = nrow(P$g(P$t0[1], x))
+ q = ncol(P$g(P$t0[1], x))
+ k = 1
+ }
+ else
+ {
+ n = nrow(P$g(P$t0, x))
+ q = ncol(P$g(P$t0, x))
+ k = length(P$t0)
+ }
+ k2 <- k
+ df <- q - k
+ w=diag(q)
+ if (P$optfct == "optim")
+ res <- optim(P$t0, .obj1, x = P$x, w = w, gf = P$g, ...)
+ if (P$optfct == "nlminb")
+ {
+ res <- nlminb(P$t0, .obj1, x = P$x, w = w, gf = P$g, ...)
+ res$value <- res$objective
+ }
+ if (P$optfct == "optimize")
+ {
+ res <- optimize(.obj1, P$t0, x = P$x, w = w, gf = P$g, ...)
+ res$par <- res$minimum
+ res$value <- res$objective
+ }
+
+ if (q == k2 | P$wmatrix == "ident")
+ {
+ z <- list(coefficients = res$par, objective = res$value, k=k, k2=k2, n=n, q=q, df=df)
+ }
+ else
+ {
+ ch <- 100000
+ j <- 1
+ while(ch > P$crit)
+ {
+ tet <- res$par
+ if (P$vcov == "iid")
+ w <- P$iid(tet, P$x, P$g)
+ if (P$vcov == "HAC")
+ w <- HAC(P$g(tet, P$x), kernel = P$kernel, bw = P$bw, prewhite = P$prewhite,
+ ar.method = P$ar.method, approx = P$approx, tol = P$tol)
+
+ if (P$optfct == "optim")
+ res <- optim(tet, .obj1, x = P$x, w = w, gf = P$g, ...)
+ if (P$optfct == "nlminb")
+ {
+ res <- nlminb(tet, .obj1, x = P$x, w = w, gf = P$g, ...)
+ res$value <- res$objective
+ }
+ if (P$optfct == "optimize")
+ {
+ res <- optimize(.obj1, P$t0, x = P$x, w = w, gf = P$g, ...)
+ res$par <- res$minimum
+ res$value <- res$objective
+ }
+ ch <- crossprod(abs(tet-res$par)/tet)^.5
+ if (j>P$itermax)
+ {
+ cat("No convergence after ", P$itermax, " iterations")
+ ch <- P$crit
+ }
+ j <- j+1
+ }
+ z = list(coefficients = res$par, objective = res$value,k=k, k2=k2, n=n, q=q, df=df)
+ }
+
+ z$x <- P$x
+ z$gt <- P$g(z$coefficients, P$x)
+ z$gradv <- P$gradv
+ z$iid <- P$iid
+ z$g <- P$g
+
+ class(z) <- paste(P$TypeGmm,".res",sep="")
+ return(z)
+ }
+
+momentEstim.baseGmm.cue.formula <- function(object, ...)
+ {
+ P <- object
+ g <- P$g
+ dat <- getDat(P$gform, P$x)
+ x <- dat$x
+ k <- dat$k
+ k2 <- k*dat$ny
+ n <- nrow(x)
+ q <- dat$ny*dat$nh
+ df <- q-k*dat$ny
+
+ if (q == k2 | P$wmatrix == "ident")
+ {
+ w <- diag(q)
+ res <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, g)
+ z = list(coefficients = res$par, objective = res$value, dat = dat, k = k, k2 = k2, n = n, q = q, df = df)
+ }
+ else
+ {
+ if (is.null(P$t0))
+ P$t0 <- .tetlin(x,diag(q), dat$ny, dat$nh, dat$k, P$gradv, g)$par
+ if (P$optfct == "optim")
+ res2 <- optim(P$t0,.objCue, x = x, P = P, ...)
+ if (P$optfct == "nlminb")
+ {
+ res2 <- nlminb(P$t0,.objCue, x = x, P = P, ...)
+ res2$value <- res2$objective
+ }
+ if (P$optfct == "optimize")
+ {
+ res2 <- optimize(.objCue,P$t0, x = x, P = P, ...)
+ res2$par <- res2$minimum
+ res2$value <- res2$objective
+ }
+ z = list(coefficients = res2$par, objective = res2$value, dat = dat, k = k, k2 = k2, n = n, q = q, df = df)
+ }
+
+ z$gt <- g(z$coefficients, x)
+ b <- z$coefficients
+ y <- as.matrix(model.response(dat$mf, "numeric"))
+ ny <- dat$ny
+ b <- t(matrix(b, nrow = dat$ny))
+ x <- as.matrix(model.matrix(dat$mt, dat$mf, NULL))
+ yhat <- x %*% b
+ z$dat <- dat
+ z$fitted.values <- yhat
+ z$residuals <- y - yhat
+ z$terms <- dat$mt
+ if(P$model) z$model <- dat$mf
+ if(P$X) z$x <- x
+ if(P$Y) z$y <- y
+ z$gradv <- P$gradv
+ z$iid <- P$iid
+ z$g <- P$g
+
+ namex <- colnames(dat$x[,(dat$ny+1):(dat$ny+dat$k)])
+ nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
+
+ if (dat$ny > 1)
+ {
+ namey <- colnames(dat$x[,1:dat$ny])
+ names(z$coefficients) <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
+ colnames(z$gt) <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny, dat$nh)), sep = "")
+ }
+
+ if (dat$ny == 1)
+ {
+ names(z$coefficients) <- namex
+ colnames(z$gt) <- nameh
+ }
+ class(z) <- paste(P$TypeGmm,".res",sep="")
+ return(z)
+ }
+
+momentEstim.baseGmm.cue <- function(object, ...)
+ {
+ P <- object
+ x <- P$x
+ if (P$optfct == "optimize")
+ {
+ n = nrow(P$g(P$t0[1], x))
+ q = ncol(P$g(P$t0[1], x))
+ k = 1
+ }
+ else
+ {
+ n = nrow(P$g(P$t0, x))
+ q = ncol(P$g(P$t0, x))
+ k = length(P$t0)
+ }
+ k2 <- k
+ df <- q - k
+ w=diag(q)
+ if (P$optfct == "optim")
+ res <- optim(P$t0, .obj1, x = P$x, w = w, gf = P$g, ...)
+ if (P$optfct == "nlminb")
+ {
+ res <- nlminb(P$t0, .obj1, x = P$x, w = w, gf = P$g, ...)
+ res$value <- res$objective
+ }
+ if (P$optfct == "optimize")
+ {
+ res <- optimize(.obj1, P$t0, x = P$x, w = w, gf = P$g, ...)
+ res$par <- res$minimum
+ res$value <- res$objective
+ }
+
+ if (q == k2 | P$wmatrix == "ident")
+ {
+ z <- list(coefficients = res$par, objective = res$value, k=k, k2=k2, n=n, q=q, df=df)
+ }
+ else
+ {
+ if (P$optfct == "optim")
+ res2 <- optim(P$t0, .objCue, x = x, P = P, ...)
+ if (P$optfct == "nlminb")
+ {
+ res2 <- nlminb(P$t0, .objCue, x = x, P = P, ...)
+ res2$value <- res2$objective
+ }
+ if (P$optfct == "optimize")
+ {
+ res2 <- optimize(.objCue,P$t0, x = x, P = P, ...)
+ res2$par <- res2$minimum
+ res2$value <- res2$objective
+ }
+ z = list(coefficients=res2$par,objective=res2$value, k=k, k2=k2, n=n, q=q, df=df)
+ }
+
+ z$x <- P$x
+ z$gt <- P$g(z$coefficients, P$x)
+ z$gradv <- P$gradv
+ z$iid <- P$iid
+ z$g <- P$g
+
+ class(z) <- paste(P$TypeGmm,".res",sep="")
+ return(z)
+ }
+
+
+
Modified: pkg/gmm/inst/doc/empir.bib
===================================================================
--- pkg/gmm/inst/doc/empir.bib 2009-12-03 06:30:24 UTC (rev 10)
+++ pkg/gmm/inst/doc/empir.bib 2009-12-17 05:48:10 UTC (rev 11)
@@ -18,16 +18,53 @@
Number={}
}
- at article{nolan09,
- AUTHOR={Nolan, J. P.},
- TITLE={Stable Distributions},
- JOURNAL={Math/Stat Department, American University},
- VOLUME={},
- PAGES={},
- YEAR={2009},
- Number={}
- }
+ at book{nolan09,
+author = {J. P. Nolan},
+title = {Stable Distributions - Models for Heavy Tailed Data},
+year = {2010},
+publisher = {Birkh\"auser},
+address = {},
+note = {In progress, Chapter 1 online at academic2.american.edu/$\sim$jpnolan}
+}
+ @Manual{timeDate,
+ title = {timeDate: Rmetrics - Chronological and Calendarical Objects},
+ author = {Diethelm Wuertz and Yohan Chalabi with contributions from Martin Maechler and Joe W. Byers and others},
+ year = {2009},
+ note = {R package version 2100.86},
+ url = {http://CRAN.R-project.org/package=timeDate},
+ }
+
+ at Manual{timeSeries,
+ title = {timeSeries: Rmetrics - Financial Time Series Objects},
+ author = {Diethelm Wuertz and Yohan Chalabi},
+ year = {2009},
+ note = {R package version 2100.84},
+ url = {http://CRAN.R-project.org/package=timeSeries},
+ }
+
+
+
+ at Book{MASS,
+ title = {Modern Applied Statistics with S},
+ author = {W. N. Venables and B. D. Ripley},
+ publisher = {Springer},
+ edition = {Fourth},
+ address = {New York},
+ year = {2002},
+ note = {ISBN 0-387-95457-0},
+ url = {http://www.stats.ox.ac.uk/pub/MASS4}
+ }
+
+ at book{hall05,
+author = {A. R. Hall},
+title = {Generalized Method of Moments (Advanced Texts in Econometrics)},
+year = {2005},
+publisher = {Oxford University Press},
+address = {},
+}
+
+
@Book{wooldridge02,
author = {Wooldridge, J. M.},
title = {Econometric Analysis of Cross Section and Panel Data},
@@ -36,7 +73,7 @@
}
@Book{cochrane01,
-author = {Cochrane, J.H.},
+author = {Cochrane, J. H.},
title = {Asset Pricing},
publisher = {Princeton University Press},
year = {2001},
@@ -76,14 +113,14 @@
}
@Book{campbell-lo-mackinlay96,
-author = {Campbell, J. Y. and Lo, A. W. and Mackinlay, A.C.},
+author = {Campbell, J. Y. and Lo, A. W. and Mackinlay, A. C.},
title = {The Econometrics of Financial Markets},
publisher = {Princeton University Press},
year = {1996},
}
@article{newey-west94,
- AUTHOR={Newey, W.K. and West, K.D.},
+ AUTHOR={Newey, W. K. and West, K. D.},
TITLE={Automatic Lag Selection in Covariance Matrix Estimation},
JOURNAL={Review of Economic Studies},
VOLUME={61},
@@ -94,7 +131,7 @@
@article{andrews-monahan92,
- AUTHOR={Andrews, W.K. and Monahan, J. C.},
+ AUTHOR={Andrews, W. K. and Monahan, J. C.},
TITLE={An Improved Heteroskedasticity and Autocorrelation Consistent Covariance Matrix Estimator},
JOURNAL={Econometrica},
VOLUME={60},
@@ -106,14 +143,14 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 11
More information about the Gmm-commits
mailing list