[Gmm-commits] r232 - pkg/momentfit/R

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

Author: chaussep
Date: 2024-03-25 22:02:23 +0100 (Mon, 25 Mar 2024)
New Revision: 232

Modified:
   pkg/momentfit/R/gel.R
   pkg/momentfit/R/gelfit-methods.R
   pkg/momentfit/R/gmmfit-methods.R
Log:
working on lambda algotihm for GEL

Modified: pkg/momentfit/R/gel.R
===================================================================
--- pkg/momentfit/R/gel.R	2024-03-22 21:10:17 UTC (rev 231)
+++ pkg/momentfit/R/gel.R	2024-03-25 21:02:23 UTC (rev 232)
@@ -136,7 +136,16 @@
             algo <- "Wu"
     }
     algo <- match.arg(algo)
-    gmat <- as.matrix(gmat)    
+    gmat <- as.matrix(gmat)
+    chk1 <- any(apply(gmat, 2, function(x) all(x>0) | all(x<0)))
+    chk2 <- any(is.na(gmat))
+    chk3 <- any(!is.finite(gmat))
+    if (chk1 | chk2 | chk3)
+    {
+        return(list(lambda = as.numeric(rep(NA, ncol(gmat))),
+                    convergence = list(convergence=1, message='gt has some wrong values'),
+                    obj= NA))
+    }
     if (length(restrictedLam))
     {
         if (length(restrictedLam) > ncol(gmat))

Modified: pkg/momentfit/R/gelfit-methods.R
===================================================================
--- pkg/momentfit/R/gelfit-methods.R	2024-03-22 21:10:17 UTC (rev 231)
+++ pkg/momentfit/R/gelfit-methods.R	2024-03-25 21:02:23 UTC (rev 232)
@@ -207,6 +207,15 @@
 setMethod("vcov", "gelfit",
           function(object, withImpProb=FALSE, tol=1e-10, robToMiss=FALSE) {
               spec <- modelDims(object at model)
+              if (any(is.na(object at lambda)))
+              {
+                  Sigma <- matrix(NA, spec$k, spec$k)
+                  if (spec$k>0)
+                      dimnames(Sigma) <- list(spec$parNames, spec$parNames)
+                  SigmaLam <- matrix(NA, spec$q, spec$q)
+                  dimnames(SigmaLam) <- list(spec$momNames, spec$momNames)
+                  return(list(vcov_par = Sigma, vcov_lambda = SigmaLam))
+              }
               if (robToMiss)
               {
                   eta <- c(coef(object), object at lambda)
@@ -215,8 +224,13 @@
                   fit <- evalGmm(mod, theta=eta)
                   v <- vcov(fit)
                   spec <- modelDims(object at model)
-                  Sigma <- v[1:spec$k, 1:spec$k]
-                  dimnames(Sigma) <- list(spec$parNames, spec$parNames)
+                  if (spec$k>0)
+                  {
+                      Sigma <- v[1:spec$k, 1:spec$k]
+                      dimnames(Sigma) <- list(spec$parNames, spec$parNames)
+                  } else {
+                      Sigma <- matrix(nrow=0, ncol=0)
+                  }
                   SigmaLam <- v[(spec$k+1):nrow(v), (spec$k+1):ncol(v)]
                   dimnames(SigmaLam) <- list(spec$momNames, spec$momNames)
                   return(list(vcov_par = Sigma, vcov_lambda = SigmaLam))
@@ -237,22 +251,47 @@
                   G <- G/k[1]
                   gt <- gt * sqrt(bw/k[2]/n)
               }
-              qrGt <- qr(gt)
-              piv <- qrGt$pivot
-              R <- qr.R(qrGt)
-              X <- forwardsolve(t(R), G[piv,])
-              Y <- forwardsolve(t(R), diag(q)[piv,])
-              res <- lm.fit(as.matrix(X), Y)
-              u <- res$residuals
-              Sigma <- chol2inv(res$qr$qr)/n
-              diag(Sigma)[diag(Sigma) < 0] <- tol
-              if (q == ncol(G)) {
+              if (ncol(G) > 0)
+              {
+                  qrGt <- qr(gt)
+                  piv <- qrGt$pivot
+                  R <- qr.R(qrGt)
+                  X <- forwardsolve(t(R), G[piv,])
+                  Y <- forwardsolve(t(R), diag(q)[piv,])
+                  res <- lm.fit(as.matrix(X), Y)
+                  u <- res$residuals
+                  Sigma <- try(chol2inv(res$qr$qr)/n, silent=TRUE)
+                  if (inherits(Sigma, "try-error"))
+                  {
+                      Sigma <- matrix(NA, ncol(G), ncol(G))
+                      warning("Failed to compute the variance of the coefficients")
+                  } else {
+                      diag(Sigma)[diag(Sigma) < 0] <- tol
+                  }
+                  dimnames(Sigma) <- list(spec$parNames,spec$parNames)                  
+              } else {
+                  Sigma <- matrix(nrow=0, ncol=0)
+              }
+              if (q == ncol(G))
+              {
                   SigmaLam <- matrix(0, q, q)
               } else {
-                  SigmaLam <- crossprod(Y, u)/n * bw^2
-                  diag(SigmaLam)[diag(SigmaLam) < 0] <- tol
+                  if (ncol(G) > 0)
+                  {
+                      SigmaLam <- crossprod(Y, u)/n * bw^2
+                      diag(SigmaLam)[diag(SigmaLam) < 0] <- tol
+                  } else {
+                      SigmaLam <- try(solve(crossprod(gt))/n * bw^2, silent=TRUE)
+                      if (inherits(SigmaLam, "try-error"))
+                      {
+                          SigmaLam <- matrix(NA, q, q)
+                          warning("Failed to compute the variance of the lambdas")
+                      } else {
+                          diag(SigmaLam)[diag(SigmaLam) < 0] <- tol
+                      }
+                  }
               }
-              piv <- sort.int(piv, index.return = TRUE)$ix
+              dimnames(SigmaLam) <- list(spec$momNames, spec$momNames)
               list(vcov_par = Sigma, vcov_lambda = SigmaLam)
           })
 

Modified: pkg/momentfit/R/gmmfit-methods.R
===================================================================
--- pkg/momentfit/R/gmmfit-methods.R	2024-03-22 21:10:17 UTC (rev 231)
+++ pkg/momentfit/R/gmmfit-methods.R	2024-03-25 21:02:23 UTC (rev 232)
@@ -47,6 +47,13 @@
 setMethod("vcov", "gmmfit",
           function(object, sandwich=NULL, df.adj=FALSE, breadOnly=FALSE,
                    modelVcov=NULL) {
+              if (length(coef(object)) == 0)
+              {
+                  vcov <- matrix(nrow=0, ncol=0)
+                  attr(vcov, "type") <- list(sandwich=sandwich, df.adj=df.adj,
+                                             breadOnly=breadOnly)
+                  return(vcov)
+              }
               if (!is.null(modelVcov))
                   {
                       if (modelVcov != object at model@vcov)