[Gmm-commits] r34 - in pkg/gmm: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 26 21:41:06 CET 2010
Author: chaussep
Date: 2010-11-26 21:41:05 +0100 (Fri, 26 Nov 2010)
New Revision: 34
Modified:
pkg/gmm/NEWS
pkg/gmm/R/getModel.R
pkg/gmm/R/gmm.R
pkg/gmm/R/momentEstim.R
pkg/gmm/man/gmm.Rd
Log:
Added 2SLS with iid errors and added the associated example to gmm manual
Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS 2010-10-15 22:19:22 UTC (rev 33)
+++ pkg/gmm/NEWS 2010-11-26 20:41:05 UTC (rev 34)
@@ -4,9 +4,9 @@
(like before) or as a formula. See details and examples in gmm help
o A data argument is added to the gmm function. Therefore, it is no longer required
that the variables in data.frames be attached before using gmm(). You just need to add the option data=your_data.frame.
+o 2SLS is now implemented in a more efficient way for linear models. Just add the option vcov="iid". An example as been included for that case
-
Changes in version 1.3-3
o It is a very small modification to avoid errors in the installation. The function linear.hypothesis of the car package is now deprecated. They have be changed to linearHypothesis.
Modified: pkg/gmm/R/getModel.R
===================================================================
--- pkg/gmm/R/getModel.R 2010-10-15 22:19:22 UTC (rev 33)
+++ pkg/gmm/R/getModel.R 2010-11-26 20:41:05 UTC (rev 34)
@@ -27,7 +27,15 @@
dat <- getDat(object$g, object$x, object$data)
if(is.null(object$weightsMatrix))
- clname <- paste(class(object), ".", object$type, ".formula", sep = "")
+ {
+ if (object$vcov == "iid")
+ {
+ clname <- "baseGmm.twoStep.formula"
+ object$type <- "Linear model with iid errors: Regular IV or 2SLS"
+ }
+ else
+ clname <- paste(class(object), ".", object$type, ".formula", sep = "")
+ }
else
{
clname <- "fixedW.formula"
Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R 2010-10-15 22:19:22 UTC (rev 33)
+++ pkg/gmm/R/gmm.R 2010-11-26 20:41:05 UTC (rev 34)
@@ -98,19 +98,45 @@
}
-.tetlin <- function(x, w, ny, nh, k, gradv, g)
+.tetlin <- function(x, w, ny, nh, k, gradv, g, type=NULL)
{
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))
+ if (!is.null(type))
+ {
+ if(type=="2sls")
+ {
+ fsls <- lm(xm~hm-1)$fitted
+ par <- lm(ym~fsls-1)$coef
+ if (ny == 1)
+ {
+ e2sls <- ym-xm%*%par
+ v2sls <- lm(e2sls~hm-1)$fitted
+ value <- sum(v2sls^2)/sum(e2sls^2)
+ }
+ else
+ {
+ par <- c(t(par))
+ g2sls <- g(par, x, ny, nh, k)
+ w <- crossprod(g2sls)/n
+ gb <- matrix(colMeans(g2sls), ncol = 1)
+ value <- crossprod(gb, solve(w, gb))
+ }
+ }
+ }
+ else
+ {
+ 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)
}
Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R 2010-10-15 22:19:22 UTC (rev 33)
+++ pkg/gmm/R/momentEstim.R 2010-11-26 20:41:05 UTC (rev 34)
@@ -129,12 +129,14 @@
}
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, P$centeredVcov)
+ {
+ res2 <- .tetlin(x, diag(q), dat$ny, dat$nh, dat$k, P$gradv, P$g, type="2sls")
+ }
if (P$vcov == "HAC")
{
+ w=diag(q)
+ res1 <- .tetlin(x, w, dat$ny, dat$nh, dat$k, P$gradv, P$g)
if(P$centeredVcov)
gmat <- lm(g(res1$par, x)~1)
else
@@ -144,9 +146,9 @@
}
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(x, w, dat$ny, dat$nh, dat$k, P$gradv, g)
}
- 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)
Modified: pkg/gmm/man/gmm.Rd
===================================================================
--- pkg/gmm/man/gmm.Rd 2010-10-15 22:19:22 UTC (rev 33)
+++ pkg/gmm/man/gmm.Rd 2010-11-26 20:41:05 UTC (rev 34)
@@ -228,5 +228,32 @@
}
# Now we want to estimate the two parameters using the GMM.
gmm(g, x, c(0, 0), grad = Dg)
+
+# Two-stage-least-squares (2SLS), or IV with iid errors.
+# The model is:
+# Y(t) = b[0] + b[1]C(t) + b[2]Y(t-1) + e(t)
+# e(t) is an MA(1)
+# The instruments are Z(t)={1 C(t) y(t-2) y(t-3) y(t-4)}
+
+getdat <- function(n) {
+e <- arima.sim(n,model=list(ma=.9))
+C <- runif(n,0,5)
+Y <- rep(0,n)
+Y[1] = 1 + 2*C[1] + e[1]
+for (i in 2:n){
+Y[i] = 1 + 2*C[i] + 0.9*Y[i-1] + e[i]
}
+Yt <- Y[5:n]
+X <- cbind(1,C[5:n],Y[4:(n-1)])
+Z <- cbind(1,C[5:n],Y[3:(n-2)],Y[2:(n-3)],Y[1:(n-4)])
+return(list(Y=Yt,X=X,Z=Z))
+}
+d <- getdat(5000)
+res4 <- gmm(d$Y~d$X-1,~d$Z-1,vcov="iid")
+res4
+
+
+
+}
+
More information about the Gmm-commits
mailing list