[Gmm-commits] r9 - in pkg/gmm: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 3 04:47:05 CET 2009
Author: chaussep
Date: 2009-12-03 04:47:05 +0100 (Thu, 03 Dec 2009)
New Revision: 9
Added:
pkg/gmm/R/HAC.R
pkg/gmm/R/Methods.gel.R
pkg/gmm/R/Methods.gmm.R
pkg/gmm/R/PlotMethods.R
pkg/gmm/R/charStable.R
pkg/gmm/man/coef.Rd
pkg/gmm/man/confint.Rd
pkg/gmm/man/fitted.Rd
pkg/gmm/man/formula.Rd
pkg/gmm/man/getDat.Rd
pkg/gmm/man/getLamb.Rd
pkg/gmm/man/plot.Rd
pkg/gmm/man/print.Rd
pkg/gmm/man/residuals.Rd
pkg/gmm/man/smoothG.Rd
pkg/gmm/man/summary.Rd
pkg/gmm/man/vcov.Rd
Log:
Stucture modifications
Added: pkg/gmm/R/HAC.R
===================================================================
--- pkg/gmm/R/HAC.R (rev 0)
+++ pkg/gmm/R/HAC.R 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,204 @@
+# 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/
+
+kweights2 <- function(x, kernel = c("Truncated", "Bartlett", "Parzen",
+ "Tukey-Hanning", "Quadratic Spectral"), normalize = FALSE)
+{
+ kernel <- match.arg(kernel)
+ if(normalize) {
+ ca <- switch(kernel,
+ "Truncated" = 2,
+ "Bartlett" = 2/3,
+ "Parzen" = .539285,
+ "Tukey-Hanning" = 3/4,
+ "Quadratic Spectral" = 1)
+ } else ca <- 1
+
+ switch(kernel,
+ "Truncated" = { ifelse(ca * x > 1, 0, 1) },
+ "Bartlett" = { ifelse(ca * x > 1, 0, 1 - abs(ca * x)) },
+ "Parzen" = {
+ ifelse(ca * x > 1, 0, ifelse(ca * x < 0.5,
+ 1 - 6 * (ca * x)^2 + 6 * abs(ca * x)^3, 2 * (1 - abs(ca * x))^3))
+ },
+ "Tukey-Hanning" = {
+ ifelse(ca * x > 1, 0, (1 + cos(pi * ca * x))/2)
+ },
+ "Quadratic Spectral" = {
+ y <- 6 * pi * x/5
+ ifelse(x < 1e-4, 1, 3 * (1/y)^2 * (sin(y)/y - cos(y)))
+ })
+}
+
+
+
+bwNeweyWest2 <- function (x, kernel = c("Bartlett", "Parzen",
+ "Quadratic Spectral", "Truncated", "Tukey-Hanning"),
+ prewhite = 1, ar.method = "ols",...)
+{
+ kernel <- match.arg(kernel)
+ if (kernel %in% c("Truncated", "Tukey-Hanning"))
+ stop(paste("Automatic bandwidth selection only available for ",
+ dQuote("Bartlett"), ", ", dQuote("Parzen"), " and ",
+ dQuote("Quadratic Spectral"), " kernel. Use ", sQuote("bwAndrews2"),
+ " instead.", sep = ""))
+ prewhite <- as.integer(prewhite)
+ n <- nrow(x)
+ k <- ncol(x)
+ weights <- rep(1, k)
+ if (length(weights) < 2)
+ weights <- 1
+ mrate <- switch(kernel, Bartlett = 2/9, Parzen = 4/25, "Quadratic Spectral" = 2/25)
+ m <- floor(ifelse(prewhite > 0, 3, 4) * (n/100)^mrate)
+ if (prewhite > 0) {
+ x <- as.matrix(na.omit(ar(x, order.max = prewhite,
+ demean = FALSE, aic = FALSE, method = ar.method)$resid))
+ n <- n - prewhite
+ }
+ hw <- x %*% weights
+ sigmaj <- function(j) sum(hw[1:(n - j)] * hw[(j + 1):n])/n
+ sigma <- sapply(0:m, sigmaj)
+ s0 <- sigma[1] + 2 * sum(sigma[-1])
+ s1 <- 2 * sum(1:m * sigma[-1])
+ s2 <- 2 * sum((1:m)^2 * sigma[-1])
+ qrate <- 1/(2 * ifelse(kernel == "Bartlett", 1, 2) + 1)
+ rval <- switch(kernel, Bartlett = {
+ 1.1447 * ((s1/s0)^2)^qrate
+ }, Parzen = {
+ 2.6614 * ((s2/s0)^2)^qrate
+ }, "Quadratic Spectral" = {
+ 1.3221 * ((s2/s0)^2)^qrate
+ })
+ rval <- rval * (n + prewhite)^qrate
+ return(rval)
+}
+
+
+HAC <- function (x, weights = weightsAndrews2, bw = bwAndrews2, prewhite = FALSE, ar.method = "ols", kernel=c("Quadratic Spectral", "Truncated", "Bartlett", "Parzen", "Tukey-Hanning"), approx="AR(1)",tol = 1e-7)
+{
+ n.orig <- n <- nrow(x)
+ k <- ncol(x)
+ kernel=match.arg(kernel)
+ if(prewhite > 0)
+ {
+ var.fit <- ar(x, order.max = prewhite, demean = FALSE, aic = FALSE, method = ar.method)
+ if(k > 1) D <- solve(diag(ncol(x)) - apply(var.fit$ar, 2:3, sum))
+ else D <- as.matrix(1/(1 - sum(var.fit$ar)))
+ x <- as.matrix(na.omit(var.fit$resid))
+ n <- n - prewhite
+ }
+ weights <- weights(x, ar.method = ar.method,kernel=kernel,bw=bw, approx = approx, prewhite = 1, tol = tol)
+ if (length(weights) > n)
+ {
+ warning("more weights than observations, only first n used")
+ weights <- weights[1:n]
+ }
+ utu <- 0.5 * crossprod(x) * weights[1]
+ wsum <- n * weights[1]/2
+ w2sum <- n * weights[1]^2/2
+ if (length(weights) > 1) {
+ for (ii in 2:length(weights)) {
+ utu <- utu + weights[ii] * crossprod(x[1:(n -
+ ii + 1), , drop = FALSE], x[ii:n, , drop = FALSE])
+ wsum <- wsum + (n - ii + 1) * weights[ii]
+ w2sum <- w2sum + (n - ii + 1) * weights[ii]^2
+ }
+ }
+ utu <- utu + t(utu)
+
+ if(prewhite > 0) {
+ utu <- crossprod(t(D), utu) %*% t(D)
+ }
+ wsum <- 2 * wsum
+ w2sum <- 2 * w2sum
+ bc <- n^2/(n^2 - wsum)
+ df <- n^2/w2sum
+ rval <- utu/n.orig
+
+ return(rval)
+}
+
+weightsAndrews2 <- function (x, bw = bwAndrews2, kernel = c("Quadratic Spectral",
+ "Truncated", "Bartlett", "Parzen", "Tukey-Hanning"), approx = c("AR(1)",
+ "ARMA(1,1)"), prewhite = 1, ar.method = "ols", tol = 1e-7, verbose = FALSE)
+{
+ kernel <- match.arg(kernel)
+ approx=match.arg(approx)
+
+ if (is.function(bw))
+ bw <- bw(x, kernel = kernel, prewhite = prewhite, ar.method = ar.method, approx=approx)
+ n <- NROW(x)
+ weights <- kweights2(0:(n - 1)/bw, kernel = kernel)
+ weights <- weights[1:max(which(abs(weights) > tol))]
+ return(weights)
+}
+
+
+bwAndrews2 <- function (x, kernel = c("Quadratic Spectral",
+ "Truncated", "Bartlett", "Parzen", "Tukey-Hanning"), approx = c("AR(1)",
+ "ARMA(1,1)"), prewhite = 1, ar.method = "ols")
+{
+ kernel <- match.arg(kernel)
+ approx <- match.arg(approx)
+ n <- nrow(x)
+ k <- ncol(x)
+
+ if (approx == "AR(1)") {
+ fitAR1 <- function(x) {
+ rval <- ar(x, order.max = 1, aic = FALSE, method = "ols")
+ rval <- c(rval$ar, sqrt(rval$var.pred))
+ names(rval) <- c("rho", "sigma")
+ return(rval)
+ }
+ ar.coef <- apply(x, 2, fitAR1)
+ denum <- sum((ar.coef["sigma", ]/(1 - ar.coef["rho",
+ ]))^4)
+ alpha2 <- sum(4 * ar.coef["rho", ]^2 * ar.coef["sigma",
+ ]^4/(1 - ar.coef["rho", ])^8)/denum
+ alpha1 <- sum(4 * ar.coef["rho", ]^2 * ar.coef["sigma",
+ ]^4/((1 - ar.coef["rho", ])^6 * (1 + ar.coef["rho",
+ ])^2))/denum
+ }
+ else {
+ fitARMA11 <- function(x) {
+ rval <- arima(x, order = c(1, 0, 1), include.mean = FALSE)
+ rval <- c(rval$coef, sqrt(rval$sigma2))
+ names(rval) <- c("rho", "psi", "sigma")
+ return(rval)
+ }
+ arma.coef <- apply(x, 2, fitARMA11)
+ denum <- sum(((1 + arma.coef["psi", ]) * arma.coef["sigma",
+ ]/(1 - arma.coef["rho", ]))^4)
+ alpha2 <- sum(4 * ((1 + arma.coef["rho", ] *
+ arma.coef["psi", ]) * (arma.coef["rho", ] + arma.coef["psi",
+ ]))^2 * arma.coef["sigma", ]^4/(1 - arma.coef["rho",
+ ])^8)/denum
+ alpha1 <- sum(4 * ((1 + arma.coef["rho", ] *
+ arma.coef["psi", ]) * (arma.coef["rho", ] + arma.coef["psi",
+ ]))^2 * arma.coef["sigma", ]^4/((1 - arma.coef["rho",
+ ])^6 * (1 + arma.coef["rho", ])^2))/denum
+ }
+ rval <- switch(kernel, Truncated = {
+ 0.6611 * (n * alpha2)^(1/5)
+ }, Bartlett = {
+ 1.1447 * (n * alpha1)^(1/3)
+ }, Parzen = {
+ 2.6614 * (n * alpha2)^(1/5)
+ }, "Tukey-Hanning" = {
+ 1.7462 * (n * alpha2)^(1/5)
+ }, "Quadratic Spectral" = {
+ 1.3221 * (n * alpha2)^(1/5)
+ })
+ return(rval)
+}
+
Added: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R (rev 0)
+++ pkg/gmm/R/Methods.gel.R 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,170 @@
+# 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/
+
+confint.gel <- function(object, parm, level=0.95, lambda=FALSE, ...)
+ {
+ z <- object
+ n <- nrow(z$gt)
+
+ se_par <- sqrt(diag(z$vcov_par))
+ par <- z$coefficients
+ tval <- par/se_par
+
+ se_parl <- sqrt(diag(z$vcov_lambda))
+ lamb <- z$lambda
+
+ zs <- qnorm((1-level)/2,lower.tail=FALSE)
+ ch <- zs*se_par
+
+ if(!lambda)
+ {
+ ans <- cbind(par-ch,par+ch)
+ dimnames(ans) <- list(names(par),c((1-level)/2,0.5+level/2))
+ }
+ if(lambda)
+ {
+ chl <- zs*se_parl
+ ans <- cbind(lamb-chl,lamb+chl)
+ dimnames(ans) <- list(names(lamb),c((1-level)/2,0.5+level/2))
+ }
+ if(!missing(parm))
+ ans <- ans[parm,]
+ ans
+ }
+
+coef.gel <- function(object, lambda=FALSE, ...)
+ {
+ if(!lambda)
+ object$coefficients
+ else
+ object$lambda
+ }
+
+vcov.gel <- function(object, lambda=FALSE, ...)
+ {
+ if(!lambda)
+ object$vcov_par
+ else
+ object$vcov_lambda
+ }
+
+print.gel <- function(x, digits=5, ...)
+ {
+ cat("Type de GEL: ", x$type,"\n\n")
+ cat("Coefficients:\n")
+ print.default(format(coef(x), digits=digits),
+ print.gap = 2, quote = FALSE)
+ cat("\n")
+ cat("Lambdas:\n")
+ print.default(format(coef(x,lambda=TRUE), digits=digits),
+ print.gap = 2, quote = FALSE)
+ invisible(x)
+ }
+
+print.summary.gel <- function(x, digits = 5, ...)
+ {
+ cat("\nCall:\n")
+ cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
+ cat("\nType of GEL: ", x$type,"\n\n")
+ cat("Kernel: ", x$kernel,"\n\n")
+ cat("Coefficients:\n")
+ print.default(format(x$coefficients, digits=digits),
+ print.gap = 2, quote = FALSE)
+
+ cat("\nLambdas:\n")
+ print.default(format(x$lambda, digits=digits),
+ print.gap = 2, quote = FALSE)
+
+ cat("\nTests of overidentifying restrictions:\n")
+ print.default(format(x$test, digits=digits),
+ print.gap = 2, quote = FALSE)
+ cat("\nConvergence code for the coefficients: ",x$conv_par,"\n")
+ cat("\nConvergence code for the lambdas: ",x$conv_lambda,"\n")
+
+ invisible(x)
+ }
+
+summary.gel <- function(object, ...)
+ {
+ z <- object
+ n <- nrow(z$gt)
+ khat <- crossprod(z$gt)/n
+ gbar <- colMeans(z$gt)
+
+ se_par <- sqrt(diag(z$vcov_par))
+ par <- z$coefficients
+ tval <- par/se_par
+
+ se_parl <- sqrt(diag(z$vcov_lambda))
+ lamb <- z$lambda
+ tvall <- lamb/se_parl
+
+ LR_test <- 2*z$objective*n
+ LM_test <- n*crossprod(z$lambda,crossprod(khat,z$lambda))
+ J_test <- n*crossprod(gbar,solve(khat,gbar))
+ test <- c(LR_test,LM_test,J_test)
+ vptest <- pchisq(test,(ncol(z$gt)-length(z$par)),lower.tail=FALSE)
+ ans <- list(type=z$type,call=z$call)
+ names(ans$type) <-"Type of GEL"
+
+ ans$coefficients <- round(cbind(par,se_par, tval, 2 * pnorm(abs(tval), lower.tail = FALSE)),5)
+ ans$lambda <- round(cbind(lamb,se_parl, tvall, 2 * pnorm(abs(tvall), lower.tail = FALSE)),5)
+
+ dimnames(ans$coefficients) <- list(names(z$coefficients),
+ c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
+ dimnames(ans$lambda) <- list(names(z$lambda),
+ c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
+
+ ans$test <- cbind(test,vptest)
+ dimnames(ans$test) <- list(c("LR test","LM test","J test"),c("statistics","p-value"))
+
+ if (z$type == "EL")
+ ans$badrho <- z$badrho
+ if (!is.null(z$weights))
+ {
+ ans$weights <- z$weights
+ }
+ ans$conv_par <- z$conv_par
+ ans$conv_pt <- z$conv_pt
+ ans$conv_moment <- cbind(z$conv_moment)
+ ans$conv_lambda <- z$conv_lambda
+ names(ans$conv_par) <- "Convergence_code_theta"
+ names(ans$conv_pt) <- "Sum_of_pt"
+ names(ans$conv_lambda) <- "Convergence_code_for_lambda"
+ dimnames(ans$conv_moment) <- list(names(z$gt),"Sample_moment_with_pt")
+ class(ans) <- "summary.gel"
+ ans
+}
+
+residuals.gel <- function(object,...)
+ {
+ if(is.null(object$model))
+ stop("The residuals method is valid only for g=formula")
+ object$residuals
+ }
+
+fitted.gel <- function(object,...)
+ {
+ if(is.null(object$model))
+ stop("The residuals method is valid only for g=formula")
+ object$fitted.value
+ }
+
+formula.gel <- function(x, ...)
+{
+ if(is.null(x$terms))
+ stop("The gel object was not created by a formula")
+ else
+ formula(x$terms)
+}
+
Added: pkg/gmm/R/Methods.gmm.R
===================================================================
--- pkg/gmm/R/Methods.gmm.R (rev 0)
+++ pkg/gmm/R/Methods.gmm.R 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,120 @@
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# A copy of the GNU General Public License is available at
+# http://www.r-project.org/Licenses/
+
+
+summary.gmm <- function(object, ...)
+ {
+ z <- object
+ se <- sqrt(diag(z$vcov))
+ par <- z$coefficients
+ tval <- par/se
+ j <- z$objective*z$n
+ ans <- list(met=z$met,kernel=z$kernel,algo=z$algo,call=z$call)
+ names(ans$met) <- "GMM method"
+ names(ans$kernel) <- "kernel for cov matrix"
+
+ ans$coefficients <- round(cbind(par,se, tval, 2 * pnorm(abs(tval), lower.tail = FALSE)),5)
+
+ dimnames(ans$coefficients) <- list(names(z$coefficients),
+ c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
+
+ ans$J_test <- noquote(paste("Test-J degrees of freedom is ",z$df,sep=""))
+ ans$j <- noquote(cbind(j,ifelse(z$df>0,pchisq(j,z$df,lower.tail = FALSE),"*******")))
+ dimnames(ans$j) <- list("Test E(g)=0: ",c("J-test","Pz(>j)"))
+ class(ans) <- "summary.gmm"
+ ans
+ }
+
+
+formula.gmm <- function(x, ...)
+{
+ if(is.null(x$terms))
+ stop("The gmm object was not created by a formula")
+ else
+ formula(x$terms)
+}
+
+confint.gmm <- function(object, parm, level=0.95, ...)
+ {
+ z <- object
+ se <- sqrt(diag(z$vcov))
+ par <- z$coefficients
+
+ zs <- qnorm((1-level)/2,lower.tail=FALSE)
+ ch <- zs*se
+ ans <- cbind(par-ch,par+ch)
+ dimnames(ans) <- list(names(par),c((1-level)/2,0.5+level/2))
+ if(!missing(parm))
+ ans <- ans[parm,]
+ ans
+ }
+
+
+
+residuals.gmm <- function(object,...)
+ {
+ if(is.null(object$model))
+ stop("The residuals method is valid only for g=formula")
+ object$residuals
+ }
+
+fitted.gmm <- function(object,...)
+ {
+ if(is.null(object$model))
+ stop("The residuals method is valid only for g=formula")
+ object$fitted.value
+ }
+
+print.gmm <- function(x, digits=5, ...)
+ {
+ cat("Method\n", x$met,"\n\n")
+ cat("Objective function value: ",x$objective,"\n\n")
+ print.default(format(coef(x), digits=digits),
+ print.gap = 2, quote = FALSE)
+ cat("\n")
+ invisible(x)
+ }
+
+coef.gmm <- function(object,...) object$coefficients
+
+vcov.gmm <- function(object,...) object$vcov
+
+
+print.summary.gmm <- function(x, digits = 5, ...)
+ {
+ cat("\nCall:\n")
+ cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
+ cat("\nMethod: ", x$met,"\n\n")
+ cat("Kernel: ", x$kernel,"\n\n")
+ cat("Coefficients:\n")
+ print.default(format(x$coefficients, digits=digits),
+ print.gap = 2, quote = FALSE)
+
+ cat("\nJ-test:\n")
+ print.default(format(x$j, digits=digits),
+ print.gap = 2, quote = FALSE)
+ cat("\n")
+ invisible(x)
+ }
+
+
+
+
+
+
+
+
+
+
+
+
Added: pkg/gmm/R/PlotMethods.R
===================================================================
--- pkg/gmm/R/PlotMethods.R (rev 0)
+++ pkg/gmm/R/PlotMethods.R 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,102 @@
+# 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/
+
+plot.gmm <- function (x, which = c(1L:3),
+ main = list("Residuals vs Fitted values", "Normal Q-Q",
+ "Response variable and fitted values"),
+ panel = if(add.smooth) panel.smooth else points,
+ ask = prod(par("mfcol")) < length(which) && dev.interactive(), ...,
+ add.smooth = getOption("add.smooth"))
+{
+ if (!inherits(x, "gmm"))
+ stop("use only with \"gmm\" objects")
+ if (!inherits(x, "gel"))
+ {
+ if(!is.numeric(which) || any(which < 1) || any(which > 3))
+ stop("'which' must be in 1L:3")
+ show <- rep(FALSE, 3)
+ }
+ else
+ {
+ if(!is.numeric(which) || any(which < 1) || any(which > 4))
+ stop("'which' must be in 1L:4")
+ show <- rep(FALSE, 4)
+ }
+
+ show[which] <- TRUE
+ r <- residuals(x)
+ if(ncol(r)>1)
+ stop("plot.gmm is not yet implemented for system of equations")
+
+ yh <- fitted(x)
+ n <- length(r)
+
+ if (ask) {
+ oask <- devAskNewPage(TRUE)
+ on.exit(devAskNewPage(oask))
+ }
+
+ ##---------- Do the individual plots : ----------
+
+ if (show[1L]) {
+ ylim <- range(r, na.rm=TRUE)
+ ylim <- extendrange(r= ylim, f = 0.08)
+ plot(yh, r, xlab = "Fitted", ylab = "Residuals", main = main[1L],
+ ylim = ylim, type = "n", ...)
+ panel(yh, r, ...)
+ abline(h = 0, lty = 3, col = "gray")
+ }
+ if (show[2L]) { ## Normal
+ rs <- (r-mean(r))/sd(r)
+ ylim <- range(rs, na.rm=TRUE)
+ ylim[2L] <- ylim[2L] + diff(ylim) * 0.075
+ qq <- qqnorm(rs, main = main[2L], ylab = "stand. residuals", ylim = ylim, ...)
+ qqline(rs, lty = 3, col = "gray50")
+ }
+ if (show[3L]) {
+ y <- as.matrix(model.response(x$model, "numeric"))
+ ylim <- range(yh, na.rm=TRUE)
+ ylim <- extendrange(r= ylim, f = 0.08)
+ plot(y, main = main[3L],
+ ylim = ylim, ...)
+ lines(yh,col=2)
+ }
+ if (inherits(x, "gel"))
+ {
+ if (show[4L]) {
+ pt <- x$pt
+ plot(pt, type='l',main = main[4L],ylab="Implied Prob.", ...)
+ emp_pt <- rep(1/length(pt),length(pt))
+ lines(emp_pt,col=2)
+ legend("topright",c("Imp. Prob.","Empirical (1/T)"),col=1:2,lty=c(1,1))
+ }
+ }
+
+
+ invisible()
+}
+
+plot.gel <- function (x, which = c(1L:4),
+ main = list("Residuals vs Fitted values", "Normal Q-Q",
+ "Response variable and fitted values","Implied probabilities"),
+ panel = if(add.smooth) panel.smooth else points,
+ ask = prod(par("mfcol")) < length(which) && dev.interactive(), ...,
+ add.smooth = getOption("add.smooth"))
+ {
+ if (!inherits(x, "gel"))
+ stop("use only with \"gel\" objects")
+ class(x) <- c("gmm","gel")
+ plot(x,which=which,main=main,panel=panel, ask=ask, ..., add.smooth=add.smooth)
+
+ invisible()
+ }
Added: pkg/gmm/R/charStable.R
===================================================================
--- pkg/gmm/R/charStable.R (rev 0)
+++ pkg/gmm/R/charStable.R 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,81 @@
+# 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/
+
+
+charStable <- function(theta,tau,pm=0)
+ {
+ # pm is the type parametrization as described by Nolan(2009)
+ # It takes the value 0 or 1
+
+ # const can fixe parameters. It is NULL for no constraint or
+ # a matrix in which case the constraint is theta[const[,1]]=const[,2]
+
+ a <- theta[1]
+ b <- theta[2]
+ g <- theta[3]
+ d <- theta[4]
+ if(pm == 0)
+ {
+ if(a == 1)
+ {
+ if(g == 0)
+ {
+ the_car <- exp(complex(ima=d*tau))
+ }
+ else
+ {
+ re_p <- -g*abs(tau)
+ im_p <- d*tau
+ im_p[tau!=0] <- im_p[tau!=0] + re_p[tau!=0]*2/pi*b*sign(tau[tau!=0])*log(g*abs(tau[tau!=0]))
+ the_car <- exp(complex(re=re_p,ima=im_p))
+ }
+ }
+ else
+ {
+ if(g == 0)
+ {
+ the_car <- exp(complex(ima=d*tau))
+ }
+ else
+ {
+ phi <- tan(pi*a/2)
+ re_p <- -g^a*abs(tau)^a
+ im_p <- d*tau*1i
+ im_p[tau!=0] <- im_p[tau!=0] + re_p*( b*phi*sign(tau[tau!=0])*(abs(g*tau[tau!=0])^(1-a)-1) )
+ the_car <- exp(complex(re=re_p,ima=im_p))
+ }
+ }
+ }
+
+ if(pm == 1)
+ {
+ if(a == 1)
+ {
+ re_p <- -g*abs(tau)
+ im_p <- d*tau
+ im_p[tau!=0] <- im_p[tau!=0]+re_p*(b*2/pi*sign(tau[tau!=0])*log(abs(tau[tau!=0])))
+ the_car <- exp(complex(re=re_p,ima=im_p))
+ }
+ else
+ {
+ phi <- tan(pi*a/2)
+ re_p <- -g^a*abs(tau)^a
+ im_p <- re_p*(-phi*b*sign(tau))+d*tau
+ the_car <- exp(complex(re=re_p,ima=im_p))
+ }
+ }
+ return(the_car)
+ }
+
+
+
Added: pkg/gmm/man/coef.Rd
===================================================================
--- pkg/gmm/man/coef.Rd (rev 0)
+++ pkg/gmm/man/coef.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,49 @@
+\name{coef}
+\alias{coef.gel}
+\alias{coef.gmm}
+\title{Coefficients of GEL or GMM}
+\description{
+ It extracts the coefficients from \code{gel} or \code{gmm} objects.
+}
+\usage{
+\method{coef}{gmm}(object, ...)
+\method{coef}{gel}(object, lambda = FALSE, ...)
+}
+\arguments{
+ \item{object}{An object of class \code{gel} or \code{gmm} returned by the function \code{\link{gel}} or \code{\link{gmm}}}
+\item{lambda}{If set to TRUE, the lagrange multipliers are extracted instead of the vector of coefficients}
+\item{...}{Other arguments when \code{coef} is applied to an other class object}
+}
+
+\value{
+Vector of coefficients
+}
+
+
+\examples{
+
+#################
+n = 500
+phi<-c(.2,.7)
+thet <- 0
+sd <- .2
+x <- matrix(arima.sim(n=n,list(order=c(2,0,1),ar=phi,ma=thet,sd=sd)),ncol=1)
+y <- x[7:n]
+ym1 <- x[6:(n-1)]
+ym2 <- x[5:(n-2)]
+
+H <- cbind(x[4:(n-3)], x[3:(n-4)], x[2:(n-5)], x[1:(n-6)])
+g <- y ~ ym1 + ym2
+x <- H
+t0 <- c(0,.5,.5)
+
+res <- gel(g, x, t0)
+
+coef(res)
+coef(res, lambda = TRUE)
+###################
+res <- gmm(g, x)
+coef(res)
+
+}
+
Added: pkg/gmm/man/confint.Rd
===================================================================
--- pkg/gmm/man/confint.Rd (rev 0)
+++ pkg/gmm/man/confint.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,68 @@
+\name{confint}
+\alias{confint.gel}
+\alias{confint.gmm}
+\title{Confidence intervals for GMM or GEL}
+\description{
+It produces confidence intervals for the coefficients from \code{gel} or \code{gmm} estimation.
+}
+\usage{
+\method{confint}{gel}(object, parm, level = 0.95, lambda = FALSE, ...)
+\method{confint}{gmm}(object, parm, level = 0.95, ...)
+}
+\arguments{
+ \item{object}{An object of class \code{gel} or \code{gmm} returned by the function \code{\link{gel}} or \code{\link{gmm}}}
+\item{parm}{A specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector
+ of names. If missing, all parameters are considered.}
+\item{level}{The confidence level}
+\item{lambda}{If set to TRUE, the confidence intervals for the Lagrange multipliers are produced.}
+\item{...}{Other arguments when \code{confint} is applied to another classe object}
+}
+
+\value{
+It returns a matrix with the first column being the lower bound and the second the upper bound.}
+
+
+\references{
+ Hansen, L.P. (1982),
+ Large Sample Properties of Generalized Method of Moments Estimators.
+ \emph{Econometrica}, \bold{50},
+ 1029-1054,
+
+ Hansen, L.P. and Heaton, J. and Yaron, A.(1996),
+ Finit-Sample Properties of Some Alternative GMM Estimators.
+ \emph{Journal of Business and Economic Statistics}, \bold{14}
+ 262-280.
+}
+
+
+\examples{
+#################
+n = 500
+phi<-c(.2,.7)
+thet <- 0
+sd <- .2
+x <- matrix(arima.sim(n = n, list(order = c(2,0,1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+y <- x[7:n]
+ym1 <- x[6:(n-1)]
+ym2 <- x[5:(n-2)]
+
+H <- cbind(x[4:(n-3)], x[3:(n-4)], x[2:(n-5)], x[1:(n-6)])
+g <- y ~ ym1 + ym2
+x <- H
+t0 <- c(0,.5,.5)
+
+res <- gel(g, x, t0)
+
+confint(res)
+confint(res, level = 0.90)
+confint(res, lambda = TRUE)
+
+########################
+
+res <- gmm(g, x)
+
+confint(res)
+confint(res, level = 0.90)
+
+}
+
Added: pkg/gmm/man/fitted.Rd
===================================================================
--- pkg/gmm/man/fitted.Rd (rev 0)
+++ pkg/gmm/man/fitted.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,58 @@
+\name{fitted}
+\alias{fitted.gel}
+\alias{fitted.gmm}
+\title{Fitted values of GEL and GMM}
+\description{
+ Method to extract the fitted values of the model estimated by \code{\link{gel}} or \code{\link{gmm}}.
+}
+\usage{
+\method{fitted}{gel}(object, ...)
+\method{fitted}{gmm}(object, ...)
+}
+\arguments{
+ \item{object}{An object of class \code{gel} or \code{gel} returned by the function \code{\link{gel}} or \code{\link{gmm}}}
+\item{...}{Other arguments when \code{fitted} is applied to an other class object}
+}
+
+\value{
+It returns a matrix of the estimated mean \eqn{\hat{y}} in \code{g=y~x} as it is done by \code{fitted.lm}.
+}
+
+\examples{
+
+# GEL can deal with endogeneity problems
+
+n = 200
+phi<-c(.2,.7)
+thet <- 0.2
+sd <- .2
+set.seed(123)
+x <- matrix(arima.sim(n = n, list(order = c(2,0,1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+
+y <- x[7:n]
+ym1 <- x[6:(n-1)]
+ym2 <- x[5:(n-2)]
+H <- cbind(x[4:(n-3)], x[3:(n-4)], x[2:(n-5)], x[1:(n-6)])
+g <- y ~ ym1 + ym2
+x <- H
+
+res <- gel(g, x, c(0,.3,.6))
+plot(y, main = "Fitted ARMA with GEL")
+lines(fitted(res), col = 2)
+
+# GMM is like GLS for linear models without endogeneity problems
+
+set.seed(345)
+n = 200
+phi<-c(.2,.7)
+thet <- 0
+sd <- .2
+x <- matrix(arima.sim(n = n, list(order = c(2,0,1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+y <- 10 + 5*rnorm(n) + x
+
+res <- gmm(y ~ x, x)
+plot(x, y, main = "Fitted model with GMM")
+lines(x, fitted(res), col = 2)
+legend("topright", c("Y","Yhat"), col = 1:2, lty = c(1,1))
+}
+
Added: pkg/gmm/man/formula.Rd
===================================================================
--- pkg/gmm/man/formula.Rd (rev 0)
+++ pkg/gmm/man/formula.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,50 @@
+\name{formula}
+\alias{formula.gel}
+\alias{formula.gmm}
+\title{Formula method for gel and gmm objects}
+\description{
+ Method to extract the formula from \code{gel} or \code{gmm} objects.
+}
+\usage{
+\method{formula}{gel}(x, ...)
+\method{formula}{gmm}(x, ...)
+}
+\arguments{
+ \item{x}{An object of class \code{gel} or \code{gmm} returned by the function \code{\link{gel}} or \code{\link{gmm}}}
+\item{...}{Other arguments to pass to other methods}
+}
+\examples{
+
+## GEL ##
+n = 200
+phi<-c(.2,.7)
+thet <- 0.2
+sd <- .2
+set.seed(123)
+x <- matrix(arima.sim(n = n, list(order = c(2,0,1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+
+y <- x[7:n]
+ym1 <- x[6:(n-1)]
+ym2 <- x[5:(n-2)]
+H <- cbind(x[4:(n-3)], x[3:(n-4)], x[2:(n-5)], x[1:(n-6)])
+g <- y ~ ym1 + ym2
+x <- H
+
+res <- gel(g, x, c(0,.3,.6))
+formula(res)
+
+# GMM is like GLS for linear models without endogeneity problems
+
+set.seed(345)
+n = 200
+phi<-c(.2,.7)
+thet <- 0
+sd <- .2
+x <- matrix(arima.sim(n = n, list(order = c(2,0,1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+y <- 10 + 5*rnorm(n) + x
+
+res <- gmm(y ~ x, x)
+formula(res)
+
+
+}
Added: pkg/gmm/man/getDat.Rd
===================================================================
--- pkg/gmm/man/getDat.Rd (rev 0)
+++ pkg/gmm/man/getDat.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,39 @@
+\name{getDat}
+\alias{getDat}
+\title{Extracting data from a formula}
+\description{
+It extract the data from a formula y~z with instrument h and put everything in a matrix. It helps redefine the function \eqn{g(\theta,x)} that is required by \code{\link{gmm}} and \code{\link{gel}}.
+}
+\usage{
+getDat(formula, h)
+}
+\arguments{
+\item{formula}{A formula that defines the linear model to be estimated (see details).}
+\item{h}{A \eqn{n\times nh} matrix of intruments(see details).}
+}
+\details{The model to be estimated is based on the moment conditions \eqn{<h,(y-z\theta)>=0}. It adds a column of ones to z and h by default. They are removed if -1 is added to the formula}
+\value{
+x: A \eqn{n \times l} matrix, where \eqn{l = ncol(y)+ncol(z)+ncol(h)+2} if "intercept" is TRUE and \eqn{ncol(y)+ncol(z)+xcol(h)} if "intercept" is FALSE.
+
+nh: dimension of h
+
+k: dimension of z
+
+ny: dimension of y
+}
+
+\examples{
+n = 500
+phi<-c(.2, .7)
+thet <- 0.2
+sd <- .2
+x <- matrix(arima.sim(n = n, list(order = c(2, 0, 1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+y <- x[7:n]
+ym1 <- x[6:(n-1)]
+ym2 <- x[5:(n-2)]
+H <- cbind(x[4:(n-3)], x[3:(n-4)], x[2:(n-5)], x[1:(n-6)])
+
+x <- getDat(y ~ ym1 + ym2, H)
+
+}
+
Added: pkg/gmm/man/getLamb.Rd
===================================================================
--- pkg/gmm/man/getLamb.Rd (rev 0)
+++ pkg/gmm/man/getLamb.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,60 @@
+\name{getLamb}
+
+\alias{getLamb}
+
+\title{Solving for the Lagrange multipliers of Generalized Empirical Likelihood (GEL) }
+
+\description{
+ It computes the vector of Lagrange multipliers, which maximizes the GEL objective function, using an iterative Newton method.
+}
+\usage{
+getLamb(g, tet, x, type = c('EL','ET','CUE'), tol_lam = 1e-12, maxiterlam = 1000, tol_obj = 1e-7)
+}
+\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.}
+
+\item{tet}{A \eqn{k \times 1} vector of parameters at which the function \eqn{g(\theta,x)} has to be evaluated}
+
+\item{x}{The matrix or vector of data from which the function \eqn{g(\theta,x)} is computed.}
+
+\item{type}{"EL" for empirical likelihood, "ET" for exponential tilting and "CUE" for continuous updated estimator.}
+
+\item{tol_lam}{Tolerance for \eqn{\lambda} between two iterations. The algorithm stops when \eqn{\|\lambda_i -\lambda_{i-1}\|} reaches \code{tol_lam} }
+
+\item{maxiterlam}{The algorithm stops if there is no convergence after "maxiterlam" iterations.}
+
+\item{tol_obj}{Tolerance for the gradiant of the objective function. The algorithm returns a non-convergence message if \eqn{\max(|gradiant|)} does not reach \code{tol_obj}. It helps the \code{gel} algorithm to select the right space to look for \eqn{\theta}}
+
+}
+\details{
+It solves the problem \eqn{\frac{1}{n}\sum_{t=1}^n D\rho(<g(\theta,x_t),\lambda>)g(\theta,x_t)=0}.
+}
+
+\value{
+lambda: A \eqn{q\times 1} vector of Lagrange multipliers which solve the system of equations given above.
+singular: 0 for a normal solution, 1 if the algorithm does not converge and 2 if the algorithm produces a singular system, NaN or Inf values.
+\code{conv_mes}: A message with details about the convergence.
+ }
+
+\references{
+Newey, W.K. and Smith, R.J. (2004), Higher Order Properties of GMM and
+Generalized Empirical Likelihood Estimators. \emph{Econometrica}, \bold{72}, 219-255.
+}
+
+\examples{
+g <- function(tet,x)
+ {
+ n <- nrow(x)
+ u <- (x[7:n] - tet[1] - tet[2]*x[6:(n-1)] - tet[3]*x[5:(n-2)])
+ f <- cbind(u, u*x[4:(n-3)], u*x[3:(n-4)], u*x[2:(n-5)], u*x[1:(n-6)])
+ return(f)
+ }
+n = 500
+phi<-c(.2, .7)
+thet <- 0.2
+sd <- .2
+x <- matrix(arima.sim(n = n, list(order = c(2, 0, 1), ar = phi, ma = thet, sd = sd)), ncol = 1)
+getLamb(g, c(0, phi), x, type = "EL")
+}
+
+
Added: pkg/gmm/man/plot.Rd
===================================================================
--- pkg/gmm/man/plot.Rd (rev 0)
+++ pkg/gmm/man/plot.Rd 2009-12-03 03:47:05 UTC (rev 9)
@@ -0,0 +1,71 @@
+\name{plot}
+\alias{plot.gel}
+\alias{plot.gmm}
+\title{Plot Diagnostics for gel and gmm objects}
+\description{
+ It is a plot method for \code{gel} or \code{gmm} objects.
+}
+\usage{
+\method{plot}{gel}(x, which = c(1L:4),
+ main = list("Residuals vs Fitted values", "Normal Q-Q",
+ "Response variable and fitted values","Implied probabilities"),
+ panel = if(add.smooth) panel.smooth else points,
+ ask = prod(par("mfcol")) < length(which) && dev.interactive(), ...,
+ add.smooth = getOption("add.smooth"))
+
+\method{plot}{gmm}(x, which = c(1L:3),
+ main = list("Residuals vs Fitted values", "Normal Q-Q",
+ "Response variable and fitted values"),
+ panel = if(add.smooth) panel.smooth else points,
+ ask = prod(par("mfcol")) < length(which) && dev.interactive(), ...,
+ add.smooth = getOption("add.smooth"))
+
+}
+\arguments{
+ \item{x}{\code{gel} or \code{gmm} object, typically result of \code{\link{gel}} or \code{\link{gmm}}.}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/gmm -r 9
More information about the Gmm-commits
mailing list