[Gmm-commits] r240 - in pkg/momentfit: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 4 23:34:54 CEST 2024
Author: chaussep
Date: 2024-06-04 23:34:54 +0200 (Tue, 04 Jun 2024)
New Revision: 240
Added:
pkg/momentfit/R/utils.R
pkg/momentfit/src/utils.f
Modified:
pkg/momentfit/DESCRIPTION
pkg/momentfit/NEWS
pkg/momentfit/R/momentModel-methods.R
pkg/momentfit/R/momentWeights-methods.R
pkg/momentfit/man/quadra-methods.Rd
pkg/momentfit/src/momentfit.h
pkg/momentfit/src/src.c
Log:
starting to add the option of using a generalized inverse to invert the weighting matrix
Modified: pkg/momentfit/DESCRIPTION
===================================================================
--- pkg/momentfit/DESCRIPTION 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/DESCRIPTION 2024-06-04 21:34:54 UTC (rev 240)
@@ -14,7 +14,7 @@
'hypothesisTest-methods.R' 'sysMomentModel.R'
'sysMomentModel-methods.R' 'rsysMomentModel-methods.R'
'sgmmfit-methods.R' 'gmm4.R' 'gel.R' 'gelfit-methods.R' 'gel4.R' 'weak.R'
- 'minAlgo.R'
+ 'minAlgo.R' 'utils.R'
License: GPL (>= 2)
NeedsCompilation: yes
VignetteBuilder: knitr
Modified: pkg/momentfit/NEWS
===================================================================
--- pkg/momentfit/NEWS 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/NEWS 2024-06-04 21:34:54 UTC (rev 240)
@@ -8,6 +8,8 @@
o Related to the previous point, there is now an option to select a different optimization algorithm, even from different packages. It is explained in the vignette. For that matter, the argument algo is no longer defined as vector of character with the option "optim" and "nlminb". Instead, it is a new object. See the vignette.
+o The gmm estimation does not work well when the weighting matrix is singular or nearly singular. It was shown that using a generalized inverse instead of the regular inverse is preferred. It is now optional to use a Moore-Penrose inverse when ever it is needed.
+
Changes in version 0.5
o The Makevars was modified. The previous version was hard to install on most MAC because of an unnecessary flag. It is no longer the case.
Modified: pkg/momentfit/R/momentModel-methods.R
===================================================================
--- pkg/momentfit/R/momentModel-methods.R 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/R/momentModel-methods.R 2024-06-04 21:34:54 UTC (rev 240)
@@ -629,12 +629,12 @@
class(gt) <- "momentFct"
opt <- object at vcovOptions
opt$x <- gt
- w <- chol(do.call(meatCL, opt))
+ w <- chol(do.call(meatCL, opt), pivot=TRUE)
type <- "chol"
} else {
w <- vcovHAC(object, theta)
wSpec <- attr(w,"Spec")
- w <- chol(w[,])
+ w <- chol(w[,], pivot=TRUE)
type <- "chol"
}
}
Modified: pkg/momentfit/R/momentWeights-methods.R
===================================================================
--- pkg/momentfit/R/momentWeights-methods.R 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/R/momentWeights-methods.R 2024-06-04 21:34:54 UTC (rev 240)
@@ -1,37 +1,54 @@
####### Methods for gmmWeights objects
#####################################
-
################ method to compute X'WX where W is the object weights #############
setGeneric("quadra", function(w, x, y, ...) standardGeneric("quadra"))
setMethod("quadra", c("momentWeights", "matrixORnumeric", "missing"),
- function(w, x, y) {
+ function(w, x, y, genInv=FALSE) {
x <- as.matrix(x)
if (is.character(w at w))
+ {
+ obj <- crossprod(x)
+ } else {
+ if (w at type == "weights")
{
- obj <- crossprod(x)
+ obj <- crossprod(x,w at w)%*%x
} else {
- if (w at type == "weights")
- {
- obj <- crossprod(x,w at w)%*%x
- } else if (w at type == "vcov") {
+ piv <- switch(w at type,
+ chol = attr(w, "pivot"),
+ qr = w at w$pivot,
+ NULL)
+ if (genInv)
+ {
+ w2 <- switch(w at type,
+ chol = crossprod(w at w),
+ qr = crossprod(qr.R(w at w)),
+ w at w)
+ iw <- MPinv(w2)
+ if (!is.null(piv))
+ iw[piv,piv] <- iw
+ obj <- crossprod(x, iw) %*% x
+ } else {
+ if (w at type == "vcov") {
obj <- crossprod(x, solve(w at w,x))
- } else if (w at type == "chol") {
- obj <- crossprod(forwardsolve(t(w at w),x))
} else {
- R <- qr.R(w at w)
- if (!all(w at w$pivot==(1:ncol(R))))
- {
- iw <- matrix(NA, ncol(R), ncol(R))
- iw[w at w$pivot, w at w$pivot] <- chol2inv(R)
- obj <- crossprod(x,iw)%*%x
- } else {
- obj <- crossprod(forwardsolve(t(R),x))
- }
+ w2 <- switch(w at type,
+ chol = w at w,
+ qr = qr.R(w at w))
+ if (!all(piv==(1:ncol(w2))))
+ {
+ iw <- chol2inv(w2)
+ iw[piv, piv] <- iw
+ obj <- crossprod(x,iw)%*%x
+ } else {
+ obj <- crossprod(forwardsolve(t(w2),x))
+ }
}
+ }
}
+ }
if (all(dim(obj)==1))
c(obj)
else
@@ -39,66 +56,95 @@
})
setMethod("quadra", c("momentWeights", "matrixORnumeric", "matrixORnumeric"),
- function(w, x, y) {
+ function(w, x, y, genInv=FALSE) {
x <- as.matrix(x)
y <- as.matrix(y)
if (is.character(w at w))
- {
+ {
obj <- crossprod(x,y)
} else {
if (w at type == "weights")
+ {
+ obj <- crossprod(x,w at w)%*%y
+ } else {
+ piv <- switch(w at type,
+ chol = attr(w, "pivot"),
+ qr = w at w$pivot,
+ NULL)
+ if (genInv)
{
- obj <- crossprod(x,w at w)%*%y
- } else if (w at type == "vcov") {
- obj <- crossprod(x, solve(w at w,y))
- } else if (w at type == "chol") {
- T1 <- forwardsolve(t(w at w), x)
- T2 <- forwardsolve(t(w at w), y)
- obj <- crossprod(T1,T2)
+ w2 <- switch(w at type,
+ chol = crossprod(w at w),
+ qr = crossprod(qr.R(w at w)),
+ w at w)
+ iw <- MPinv(w2)
+ if (!is.null(piv))
+ iw[piv,piv] <- iw
+ obj <- crossprod(x, iw) %*% y
} else {
- R <- qr.R(w at w)
- if (!all(w at w$pivot==(1:ncol(R))))
+ if (w at type == "vcov") {
+ obj <- crossprod(x, solve(w at w,y))
+ } else {
+ w2 <- switch(w at type,
+ chol = w at w,
+ qr = qr.R(w at w))
+ if (!all(piv==(1:ncol(w2))))
{
- iw <- matrix(NA, ncol(R), ncol(R))
- iw[w at w$pivot, w at w$pivot] <- chol2inv(R)
+ iw <- chol2inv(w2)
+ iw[piv, piv] <- iw
obj <- crossprod(x,iw)%*%y
} else {
- T1 <- forwardsolve(t(R), x)
- T2 <- forwardsolve(t(R), y)
+ T1 <- forwardsolve(t(w2), x)
+ T2 <- forwardsolve(t(w2), y)
obj <- crossprod(T1,T2)
}
+ }
}
+ }
}
if (all(dim(obj)==1))
c(obj)
else
- obj
+ obj
})
setMethod("quadra", c("momentWeights", "missing", "missing"),
- function(w, x, y) {
+ function(w, x, y, genInv=FALSE) {
if (is.character(w at w))
{
obj <- "Identity"
} else {
if (w at type == "weights")
+ {
+ obj <- w at w
+ } else {
+ piv <- switch(w at type,
+ chol = attr(w, "pivot"),
+ qr = w at w$pivot,
+ NULL)
+ if (genInv)
{
- obj <- w at w
- } else if (w at type == "vcov") {
- obj <- solve(w at w)
- } else if (w at type == "chol") {
- obj <- chol2inv(w at w)
+ w2 <- switch(w at type,
+ chol = crossprod(w at w),
+ qr = crossprod(qr.R(w at w)),
+ w at w)
+ obj <- MPinv(w2)
+ if (!is.null(piv))
+ obj[piv,piv] <- obj
} else {
- R <- qr.R(w at w)
- if (!all(w at w$pivot==(1:ncol(R))))
- {
- iw <- matrix(NA, ncol(R), ncol(R))
- iw[w at w$pivot, w at w$pivot] <- chol2inv(R)
- obj <- iw
- } else {
- obj <- chol2inv(R)
- }
+ w2 <- switch(w at type,
+ chol = w at w,
+ qr = qr.R(w at w),
+ w at w)
+ if (w at type == "vcov")
+ {
+ obj <- solve(w2)
+ } else {
+ obj <- chol2inv(w2)
+ obj[piv, piv] <- obj
+ }
}
+ }
}
if (all(dim(obj)==1))
c(obj)
Added: pkg/momentfit/R/utils.R
===================================================================
--- pkg/momentfit/R/utils.R (rev 0)
+++ pkg/momentfit/R/utils.R 2024-06-04 21:34:54 UTC (rev 240)
@@ -0,0 +1,24 @@
+MPinv <- function(x)
+{
+ d <- dim(x)
+ tr <- FALSE
+ if (d[1]<d[2])
+ {
+ x <- t(x)
+ d <- d[c(2,1)]
+ tr <- TRUE
+ }
+ res <- .Fortran(F_pinv, as.double(x), as.integer(d[1]), as.integer(d[2]),
+ r=integer(1), L=double(d[2]^2))
+ L <- matrix(res$L, d[2], d[2])[,1:res$r]
+ if (dim(L)[1] != dim(L)[2])
+ {
+ M <- solve(crossprod(L))
+ ix <- crossprod(M%*%t(L))%*%t(x)
+ } else {
+ ix <- crossprod(forwardsolve(L, diag(dim(L)[1])))%*%t(x)
+ }
+ if (tr)
+ ix = t(ix)
+ ix
+}
Modified: pkg/momentfit/man/quadra-methods.Rd
===================================================================
--- pkg/momentfit/man/quadra-methods.Rd 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/man/quadra-methods.Rd 2024-06-04 21:34:54 UTC (rev 240)
@@ -15,16 +15,18 @@
\code{momentWeights} object ~~
}
\usage{
-\S4method{quadra}{momentWeights,missing,missing}(w, x, y)
+\S4method{quadra}{momentWeights,missing,missing}(w, x, y, genInv=FALSE)
-\S4method{quadra}{momentWeights,matrixORnumeric,missing}(w, x, y)
+\S4method{quadra}{momentWeights,matrixORnumeric,missing}(w, x, y,
+genInv=FALSE)
\S4method{quadra}{momentWeights,matrixORnumeric,matrixORnumeric}(w, x,
-y)
+y, genInv=FALSE)
-\S4method{quadra}{sysMomentWeights,matrixORnumeric,matrixORnumeric}(w, x,
-y)
+\S4method{quadra}{sysMomentWeights,matrixORnumeric,matrixORnumeric}(w,
+x, y)
+
\S4method{quadra}{sysMomentWeights,matrixORnumeric,missing}(w, x, y)
\S4method{quadra}{sysMomentWeights,missing,missing}(w, x, y)
@@ -35,7 +37,9 @@
\item{w}{An object of class \code{"momentWeights"}}
\item{x}{A matrix or numeric vector}
\item{y}{A matrix or numeric vector}
+ \item{genInv}{Should we invert the center matrix using a generalized inverse?}
}
+
\section{Methods}{
\describe{
\item{\code{signature(w = "momentWeights", x = "matrixORnumeric", y =
@@ -55,8 +59,17 @@
object. The \code{quadra} method with no \code{y} and \code{x} is
therefore a way to invert it. The same applies to system of equations
+}}}
+
+\value{
+ It returns a single numeric value.
}
-}}
+
+\references{
+ Courrieu P (2005), Fast Computation of Moore-Penrose Inverse Matrices.
+ \emph{Neural Information Processing - Letters and Reviews}, \bold{8}(2), 25--29.
+}
+
\examples{
data(simData)
@@ -66,13 +79,30 @@
gbar <- evalMoment(model1, theta)
gbar <- colMeans(gbar)
-### Onjective function of GMM with identity matrix
+### Objective function of GMM with identity matrix
wObj <- evalWeights(model1, w="ident")
quadra(wObj, gbar)
-### Onjective function of GMM with efficient weights
+### Objective function of GMM with efficient weights
wObj <- evalWeights(model1, theta)
quadra(wObj, gbar)
+### Linearly dependent instruments
+
+simData$z3 <- simData$z1+simData$z2
+model2 <- momentModel(y~x1, ~z1+z2+z3, data=simData)
+gbar2 <- evalMoment(model2, theta)
+gbar2 <- colMeans(gbar2)
+
+## A warning is printed about the singularity of the weighting matrix
+wObj <- evalWeights(model2, theta)
+
+## The regular inverse using the QR decomposition:
+quadra(wObj)
+
+## The regular inverse using the generalized inverse:
+quadra(wObj, genInv=TRUE)
+
+
}
\keyword{methods}
Modified: pkg/momentfit/src/momentfit.h
===================================================================
--- pkg/momentfit/src/momentfit.h 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/src/momentfit.h 2024-06-04 21:34:54 UTC (rev 240)
@@ -12,4 +12,7 @@
int *maxit, int *conv, double *lam,
double *pt, double *obj);
+void F77_SUB(pinv)(double *x, int *m, int *n, int *r,
+ double *l);
+
#endif
Modified: pkg/momentfit/src/src.c
===================================================================
--- pkg/momentfit/src/src.c 2024-05-31 21:00:05 UTC (rev 239)
+++ pkg/momentfit/src/src.c 2024-06-04 21:34:54 UTC (rev 240)
@@ -5,6 +5,7 @@
static const R_FortranMethodDef fortranMethods[] = {
{"wu", (DL_FUNC) &F77_SUB(wu), 9},
{"lamcuep", (DL_FUNC) &F77_SUB(lamcuep), 9},
+ {"pinv", (DL_FUNC) &F77_SUB(pinv), 5},
{NULL, NULL, 0}
};
Added: pkg/momentfit/src/utils.f
===================================================================
--- pkg/momentfit/src/utils.f (rev 0)
+++ pkg/momentfit/src/utils.f 2024-06-04 21:34:54 UTC (rev 240)
@@ -0,0 +1,46 @@
+c Generalized inverse Courrieu (2005)
+c The GINV of mxn X is X'LMML'G', with m>n
+c where M = inv(L'L) and L'L is full rank rxr
+c The subroutine only returns L and r
+c when m=n, L is non singular and the solution is
+c inv(L')inv(L)G'
+
+ subroutine pinv(x, m, n, r, l)
+ integer n, m, k, r
+ double precision x(m,n), ix(n,m), a(n,n), tol
+ double precision l(n,n), tmp(n,1)
+
+ a = matmul(transpose(x), x)
+ tol = 0.0d0
+ do i=1,n
+ if (a(i,i) .gt. 0) then
+ if (tol .eq. 0.0d0) then
+ tol = a(i,i)
+ else
+ if (a(i,i) .lt. tol) then
+ tol = a(i,i)
+ end if
+ endif
+ end if
+ end do
+ tol = tol*1.0d-9
+ l = 0.0d0
+ r = 0
+ do k=1,n
+ r = r+1
+ if (r .eq. 1) then
+ l(k:n,r) = a(k:n,k)
+ else
+ l(k:n,r) = a(k:n,k)-matmul(l(k:n,1:(r-1)), l(k,1:(r-1)))
+ end if
+ if (l(k,r) .gt. tol) then
+ l(k,r) = sqrt(l(k,r))
+ if (k .lt. n) then
+ l((k+1):n,r) = l((k+1):n,r)/l(k,r)
+ end if
+ else
+ r = r-1
+ end if
+ end do
+ end
+
More information about the Gmm-commits
mailing list