[Gmm-commits] r40 - in pkg/gmm: . R

[noreply at r-forge.r-project.org]

Author: chaussep
Date: 2011-07-19 16:21:47 +0200 (Tue, 19 Jul 2011)
New Revision: 40

Modified:
   pkg/gmm/DESCRIPTION
   pkg/gmm/NEWS
   pkg/gmm/R/gmm.R
   pkg/gmm/R/momentEstim.R
Log:
A bug with the weighting matrix and linear gmm has been fixed

Modified: pkg/gmm/DESCRIPTION
===================================================================
--- pkg/gmm/DESCRIPTION	2011-06-06 17:20:52 UTC (rev 39)
+++ pkg/gmm/DESCRIPTION	2011-07-19 14:21:47 UTC (rev 40)
@@ -1,6 +1,6 @@
 Package: gmm
-Version: 1.3-7
-Date: 2011-06-03
+Version: 1.3-8
+Date: 2011-07-19
 Title: Generalized Method of Moments and Generalized Empirical
         Likelihood
 Author: Pierre Chausse <pchausse at uwaterloo.ca>

Modified: pkg/gmm/NEWS
===================================================================
--- pkg/gmm/NEWS	2011-06-06 17:20:52 UTC (rev 39)
+++ pkg/gmm/NEWS	2011-07-19 14:21:47 UTC (rev 40)
@@ -1,3 +1,7 @@
+Changes in version 1.3-8
+
+o A bug was found in the computation of linear GMM. The weighting matrix was not use properly. It is now fixed.
+
 Changes in version 1.3-7
 
 o Until now there was no way to know whether optim or nlminb has converged in case of nonlinear GMM. The convergence code

Modified: pkg/gmm/R/gmm.R
===================================================================
--- pkg/gmm/R/gmm.R	2011-06-06 17:20:52 UTC (rev 39)
+++ pkg/gmm/R/gmm.R	2011-07-19 14:21:47 UTC (rev 40)
@@ -116,7 +116,7 @@
 }
 
 
-.tetlin <- function(x, w, ny, nh, k, gradv, g, type=NULL)
+.tetlin <- function(x, w, ny, nh, k, gradv, g, type=NULL, inv=TRUE)
   {
   n <- nrow(x)
   ym <- as.matrix(x[,1:ny])
@@ -148,8 +148,16 @@
   	{
   if (ny>1)
   	{
-     whx <- solve(w, (crossprod(hm, xm) %x% diag(ny)))
-     wvecyh <- solve(w, matrix(crossprod(ym, hm), ncol = 1))
+     if (inv) 
+	{
+	whx <- solve(w, (crossprod(hm, xm) %x% diag(ny)))
+	wvecyh <- solve(w, matrix(crossprod(ym, hm), ncol = 1))	
+	}
+     else
+        {
+	whx <- w%*% (crossprod(hm, xm) %x% diag(ny))
+	wvecyh <- w%*%matrix(crossprod(ym, hm), ncol = 1)
+        }
      dg <- gradv(NULL,x, ny, nh, k)
      xx <- crossprod(dg, whx)
      par <- solve(xx, crossprod(dg, wvecyh))
@@ -158,14 +166,21 @@
   	{   
      if (nh>k)
      	{
-     	xzwz <- crossprod(xm,hm)%*%w%*%t(hm)
+	if(inv)
+           xzwz <- crossprod(xm,hm)%*%solve(w,t(hm))	
+	else
+     	   xzwz <- crossprod(xm,hm)%*%w%*%t(hm)
      	par <- solve(xzwz%*%xm,xzwz%*%ym)	
 	     }
 	else
 		par <- solve(crossprod(hm,xm),crossprod(hm,ym))  	}
 	}
   gb <- matrix(colSums(g(par, x, ny, nh, k))/n, ncol = 1)
-  value <- crossprod(gb, solve(w, gb)) 
+  if(inv)
+	  value <- crossprod(gb, solve(w, gb)) 
+  else
+	  value <- crossprod(gb, w%*%gb) 
+
   res <- list(par = par, value = value)
   return(res)
   }

Modified: pkg/gmm/R/momentEstim.R
===================================================================
--- pkg/gmm/R/momentEstim.R	2011-06-06 17:20:52 UTC (rev 39)
+++ pkg/gmm/R/momentEstim.R	2011-07-19 14:21:47 UTC (rev 40)
@@ -911,7 +911,7 @@
       warning("The matrix of weights is not strictly positive definite")
     }
   
-  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, inv=FALSE)
   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)