[Gmm-commits] r171 - in pkg/momentfit: R vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 22 18:58:43 CEST 2020
Author: chaussep
Date: 2020-05-22 18:58:42 +0200 (Fri, 22 May 2020)
New Revision: 171
Modified:
pkg/momentfit/R/rModel-methods.R
pkg/momentfit/R/rsysMomentModel-methods.R
pkg/momentfit/R/sgmmfit-methods.R
pkg/momentfit/R/sysMomentModel-methods.R
pkg/momentfit/vignettes/gmmS4.Rnw
pkg/momentfit/vignettes/gmmS4.pdf
Log:
finalized the system of nonlinear equations
Modified: pkg/momentfit/R/rModel-methods.R
===================================================================
--- pkg/momentfit/R/rModel-methods.R 2020-05-21 19:16:47 UTC (rev 170)
+++ pkg/momentfit/R/rModel-methods.R 2020-05-22 16:58:42 UTC (rev 171)
@@ -103,22 +103,17 @@
stop("The derivative of the constraints at theta0 is either infinite or NAN")
if (qr(dR)$rank < length(R))
stop("The matrix of derivatives of the constraints is not full rank")
- rhs <- as.character(object at fRHS)
- if (!is.null(object at fLHS))
- lhs <- as.character(object at fLHS)
- else
- lhs <- NULL
+ rhs <- object at fRHS
+ lhs <- object at fLHS
+ env <- new.env()
for (r in R)
- {
- rhs <- gsub(as.character(r[2]), paste("(", as.character(r[3]),
- ")", sep=""), rhs)
- if (!is.null(lhs))
- lhs <- gsub(as.character(r[2]),
- paste("(", as.character(r[3]),
- ")", sep=""), lhs)
- }
- rhs <- parse(text=rhs)
- lhs <- parse(text=lhs)
+ {
+ assign(as.character(r[2]),
+ parse(text=paste("(", as.character(r[3]), ")", sep=""))[[1]], env)
+ }
+ rhs <- as.expression(do.call('substitute', list(rhs[[1]], env)))
+ if (!is.null(lhs))
+ lhs <- as.expression(do.call('substitute', list(lhs[[1]], env)))
k <- object at k-length(R)
parNames <- object at parNames[!(object at parNames %in% rest)]
theta0 <- object at theta0[!(object at parNames %in% rest)]
@@ -162,27 +157,20 @@
stop("The derivative of the constraints at theta0 is either infinite or NAN")
if (qr(dR)$rank < length(R))
stop("The matrix of derivatives of the constraints is not full rank")
- rhs <- list()
- lhs <- list()
- for (i in 1:length(object at fRHS))
- {
- rhs[[i]] <- as.character(object at fRHS[[i]])
- if (!is.null(object at fLHS[[i]]))
- lhs[[i]] <- as.character(object at fLHS[[i]])
- else
- lhs[[i]] <- NULL
- for (r in R)
- {
- rhs[[i]] <- gsub(as.character(r[2]), paste("(", as.character(r[3]),
- ")", sep=""), rhs[[i]])
- if (!is.null(lhs[[i]]))
- lhs[[i]] <- gsub(as.character(r[2]),
- paste("(", as.character(r[3]),
- ")", sep=""), lhs[[i]])
- }
- rhs[[i]] <- parse(text=rhs[[i]])
- lhs[[i]] <- parse(text=lhs[[i]])
- }
+ rhs <- object at fRHS
+ lhs <- object at fLHS
+ env <- new.env()
+ for (r in R)
+ {
+ assign(as.character(r[2]),
+ parse(text=paste("(", as.character(r[3]), ")", sep=""))[[1]], env)
+ }
+ for (i in 1:length(rhs))
+ {
+ rhs[[i]] <- as.expression(do.call('substitute', list(rhs[[i]][[1]], env)))
+ if (!is.null(lhs[[i]]))
+ lhs[[i]] <- as.expression(do.call('substitute', list(lhs[[i]][[1]], env)))
+ }
k <- object at k-length(R)
parNames <- object at parNames[!(object at parNames %in% rest)]
if (length(parNames)!=k)
Modified: pkg/momentfit/R/rsysMomentModel-methods.R
===================================================================
--- pkg/momentfit/R/rsysMomentModel-methods.R 2020-05-21 19:16:47 UTC (rev 170)
+++ pkg/momentfit/R/rsysMomentModel-methods.R 2020-05-22 16:58:42 UTC (rev 171)
@@ -32,33 +32,43 @@
stop("The derivative of the constraints at theta0 is either infinite or NAN")
if (qr(dR)$rank < length(R))
stop("The matrix of derivatives of the constraints is not full rank")
- rhs <- lapply(fRHS, function(ri) as.character(ri))
- lhs <- lapply(fLHS, function(ri) ifelse(is.null(ri),NULL,as.character(ri)))
+ rhs <- object at fRHS
+ lhs <- object at fLHS
+ env <- new.env()
+ varNames <- list()
for (r in R)
{
- rhs <- gsub(as.character(r[2]), paste("(", as.character(r[3]),
- ")", sep=""), rhs)
- if (!is.null(lhs))
- lhs <- gsub(as.character(r[2]),
- paste("(", as.character(r[3]),
- ")", sep=""), lhs)
+ assign(as.character(r[2]),
+ parse(text=paste("(", as.character(r[3]), ")", sep=""))[[1]], env)
}
- rhs <- lapply(rhs, function(ri) parse(text=ri))
- lhs <- lapply(lhs, function(ri) parse(text=ri))
+ for (i in 1:length(rhs))
+ {
+ rhs[[i]] <- as.expression(do.call('substitute', list(rhs[[i]][[1]], env)))
+ if (!is.null(lhs[[i]]))
+ lhs[[i]] <- as.expression(do.call('substitute', list(lhs[[i]][[1]], env)))
+ tmp <- c(all.vars(rhs[[i]]))
+ }
k <- sum(spec$k)-length(R)
parNames <- if (is.list(spec$parNames))
- lapply(spec$parNames, function(pi) pi[!(pi %in% rest)])
- else
- spec$parNames[!(spec$parNames %in% rest)]
+ {
+ lapply(spec$parNames, function(pi) pi[!(pi %in% rest)])
+ } else {
+ spec$parNames[!(spec$parNames %in% rest)]
+ }
theta0 <- if (is.list(spec$theta0))
+ {
lapply(spec$theta0, function(ti) ti[!(names(ti) %in% rest)])
- else
+ } else {
spec$theta0[!(names(spec$theta0) %in% rest)]
+ }
+
if (length(rhs) == 1)
{
rhs <- rhs[[1]]
lhs <- lhs[[1]]
}
+ if (is.list(parNames))
+ k <- sapply(parNames, length)
list(rhs=rhs, lhs=lhs, parNames=parNames,
theta0=theta0, k=k, crossEquRest=crossRest)
}
@@ -231,7 +241,8 @@
fLHS <- cst$fLHS
list(k=k, q=object at q, n=object at n, parNames=parNames,
momNames=object at momNames, theta0=theta0,
- fRHS=fRHS, fLHS=fLHS, eqnNames=object at eqnNames)
+ fRHS=fRHS, fLHS=fLHS, eqnNames=object at eqnNames,
+ isEndo=object at isEndo)
})
## printRestrict
@@ -302,7 +313,8 @@
setMethod("print", "rsnonlinearModel",
function (x, ...) {
- print(as(x, "snonlinearModel"))
+ callNextMethod()
+ cat("(The number of endogenous variables is unreliable)\n")
printRestrict(x)
cat("\n")
} )
Modified: pkg/momentfit/R/sgmmfit-methods.R
===================================================================
--- pkg/momentfit/R/sgmmfit-methods.R 2020-05-21 19:16:47 UTC (rev 170)
+++ pkg/momentfit/R/sgmmfit-methods.R 2020-05-22 16:58:42 UTC (rev 171)
@@ -88,7 +88,7 @@
sandwich <- !object at efficientGmm
meat <- meatGmm(object, sandwich)
if (!sandwich) {
- vcov <- solve(meat)/object at model@n
+ vcov <- solve(meat)/spec$n
} else {
if (all(spec$k == spec$q)) {
G <- evalDMoment(object at model, coef(object))
@@ -115,20 +115,22 @@
setMethod("meatGmm", "sgmmfit",
function(object, robust = FALSE)
- {
- G <- evalDMoment(object at model, coef(object))
- G <- .GListToMat(G)
- if (!robust) {
- wObj <- evalWeights(object at model, coef(object), "optimal")
- meat <- quadra(wObj, G)
- } else {
- wObj <- object at wObj
- v <- vcov(object at model, coef(object))
- T1 <- quadra(wObj, v, G)
- meat <- quadra(wObj, G, T1)
- }
- meat
- })
+ {
+ spec <- modelDims(object at model)
+ G <- evalDMoment(object at model, coef(object))
+ full <- all(sapply(1:length(G), function(i) ncol(G[[i]])==sum(spec$k)))
+ G <- .GListToMat(G, full)
+ if (!robust) {
+ wObj <- evalWeights(object at model, coef(object), "optimal")
+ meat <- quadra(wObj, G)
+ } else {
+ wObj <- object at wObj
+ v <- vcov(object at model, coef(object))
+ T1 <- quadra(wObj, v, G)
+ meat <- quadra(wObj, G, T1)
+ }
+ meat
+ })
## bread
@@ -199,8 +201,11 @@
stest <- specTest(object, df.adj = df.adj)
vcovType <- switch(object at model@vcov, HAC = "HAC", iid = "OLS",
MDS = "HC")
- strength <- lapply(1:neqn, function(i)
- momentStrength(object at model[i], par[[i]], vcovType))
+ strength <- lapply(1:neqn, function(i) {
+ if (inherits(object at model, "slinearModel"))
+ momentStrength(object at model[i], par[[i]], vcovType)
+ else
+ NULL})
wSpec <- object at wObj@wSpec
ans <- new("summarySysGmm", coef = coef, type = object at type,
specTest = stest, strength = strength, model = object at model,
Modified: pkg/momentfit/R/sysMomentModel-methods.R
===================================================================
--- pkg/momentfit/R/sysMomentModel-methods.R 2020-05-21 19:16:47 UTC (rev 170)
+++ pkg/momentfit/R/sysMomentModel-methods.R 2020-05-22 16:58:42 UTC (rev 171)
@@ -5,28 +5,28 @@
setMethod("print", "sysModel",
function(x, ...)
- {
- cat("System of Equations Model\n")
- cat("*************************\n")
- type <- gsub("s", "", is(x)[1])
- cat("Moment type: ", strsplit(type, "G")[[1]][1], "\n",
- sep = "")
- cat("Covariance matrix: ", x at vcov, sep = "")
- if (x at vcov == "HAC") {
- option <- x at vcovOptions
- cat(" with ", option$kernel, " kernel and ")
- if (is.numeric(option$bw))
- cat("Fixed bandwidth (", round(option$bw, 3), ")", sep = "")
- else cat(option$bw, " bandwidth", sep = "")
- }
- cat("\n")
- d <- modelDims(x)
- for (i in 1:length(d$eqnNames))
- cat(d$eqnNames[i], ": coefs=", d$k[i],
- ", moments=", d$q[i], ", number of Endogenous: ",
- sum(d$isEndo[[i]]), "\n", sep="")
- cat("Sample size: ", d$n, "\n")
- })
+ {
+ cat("System of Equations Model\n")
+ cat("*************************\n")
+ type <- gsub("s", "", is(x)[1])
+ cat("Moment type: ", strsplit(type, "G")[[1]][1], "\n",
+ sep = "")
+ cat("Covariance matrix: ", x at vcov, sep = "")
+ if (x at vcov == "HAC") {
+ option <- x at vcovOptions
+ cat(" with ", option$kernel, " kernel and ")
+ if (is.numeric(option$bw))
+ cat("Fixed bandwidth (", round(option$bw, 3), ")", sep = "")
+ else cat(option$bw, " bandwidth", sep = "")
+ }
+ cat("\n")
+ d <- modelDims(x)
+ for (i in 1:length(d$eqnNames))
+ cat(d$eqnNames[i], ": coefs=", d$k[i],
+ ", moments=", d$q[i], ", number of Endogenous: ",
+ sum(d$isEndo[[i]]), "\n", sep="")
+ cat("Sample size: ", d$n, "\n")
+ })
setMethod("show", "sysModel", function(object) print(object))
@@ -651,12 +651,12 @@
theta0 <- setCoef(object, theta0)
g <- function(theta, wObj, object){
spec <- modelDims(object)
- theta <- momentfit:::.tetReshape(theta, object at eqnNames, spec$parNames)
+ theta <- .tetReshape(theta, object at eqnNames, spec$parNames)
evalGmmObj(object, theta, wObj)
}
dg <- function(theta, wObj, object) {
spec <- modelDims(object)
- theta <- momentfit:::.tetReshape(theta, object at eqnNames, spec$parNames)
+ theta <- .tetReshape(theta, object at eqnNames, spec$parNames)
gt <- evalMoment(object, theta)
gt <- do.call(cbind, gt)
n <- nrow(gt)
@@ -663,7 +663,7 @@
gt <- colMeans(gt)
G <- evalDMoment(object, theta)
full <- all(sapply(1:length(G), function(i) ncol(G[[i]])==sum(spec$k)))
- G <- momentfit:::.GListToMat(G, full)
+ G <- .GListToMat(G, full)
obj <- 2 * n * quadra(wObj, G, gt)
obj
}
Modified: pkg/momentfit/vignettes/gmmS4.Rnw
===================================================================
--- pkg/momentfit/vignettes/gmmS4.Rnw 2020-05-21 19:16:47 UTC (rev 170)
+++ pkg/momentfit/vignettes/gmmS4.Rnw 2020-05-22 16:58:42 UTC (rev 171)
@@ -472,6 +472,24 @@
\subsection{Methods for momentModel Classes} \label{sec:momentmodel-methods}
\begin{itemize}
+\item \textit{setCoef}: A method to validate and to format the vector
+ of coefficients. It is used by most methods to verify if the vector
+ is correctly specified and to re-format it if needed. Is it
+ particularly useful to create a vector of initial values. For
+ example, if we want to create a named vector with valid names for
+ the nonlinear ``mod2'' model, we can proceed this way:
+
+<<>>=
+setCoef(mod2, 1:3)
+@
+
+The method also reorders the vector to match the order in the model
+object:
+
+<<>>=
+setCoef(mod2, c(theta1=1, theta2=1, theta0=2))
+@
+
\item \textit{residuals}: Only for ``linearModel'' and ``nonlinearModel'',
it returns $\epsilon(\theta)$:
@@ -481,13 +499,16 @@
e2 <- residuals(mod2, theta0)
@
+It returns errors if the names are invalid or the number of
+coefficients is wrong.
+
\item \textit{Dresiduals}: Only for ``linearModel'' and
``nonlinearModel'', it returns the $n\times k$ matrix
$d\epsilon(\theta)/d\theta$:
<<>>=
-theta0 <- c(theta0=1, theta1=1, theta2=2)
e1 <- Dresiduals(mod1)
+theta0 <- setCoef(mod2, c(1,1,2))
e2 <- Dresiduals(mod2, theta0)
@
@@ -564,6 +585,8 @@
above):
<<>>=
theta0 <- c(theta0=.1, theta1=1, theta2=-2)
+## or ##
+theta0 <- setCoef(mod2, c(.1,1,-2))
evalDMoment(mod2, theta0)
@
@@ -706,6 +729,13 @@
coef(R2.mod1, c(1.5,.5))
@
+It is possible to verify that the length or names are valid by using
+the \textit{setCoef} method:
+
+<<>>=
+setCoef(R2.mod1, c(1.5,.5))
+@
+
Notice that any restricted class object contains its unrestricted
version. For example, ``rlinearModel'' is a class that contains a
``linearModel'' class object plus a few additional slots. We can
@@ -736,7 +766,8 @@
R1 <- c("theta1=theta2^2")
R1.mod2 <- restModel(mod2, R1)
evalDMoment(mod2, c(theta0=1, theta1=1, theta2=1))
-evalDMoment(R1.mod2, c(theta0=1, theta2=1))
+## with setCoef:
+evalDMoment(R1.mod2, setCoef(R1.mod2, c(1,1)))
@
Every method uses the method \textit{modelDims} to extract the
@@ -2002,6 +2033,19 @@
``momentModel'' classes. Here, we briefly describe the difference.
\begin{itemize}
+\item \textit{setCoef}: As for single-equation models, it validate and
+ organize the list of coefficients. It is very helpful for large
+ systems. In the ``smod1'' system, we have 11 coefficients. We can
+ create the list using a simple vector:
+
+<<>>=
+setCoef(smod1, 1:11)
+@
+
+If it is a named list of names vectors, the method match the order of
+the model. Or course, it also make sure the dimensions and names are
+valid.
+
\item \textit{[}: The method has two arguments. The first is an vector
of integers to select the equations, and the second is a list of
integers to select the instruments in each of the selected
@@ -2211,12 +2255,30 @@
nsmod at parNames
@
+We can use the \textit{setCoef} method to create valid vectors:
+<<>>=
+setCoef(nsmod, 1:11)
+@
+Creating a restricted model with and without cross-equation
+restrictions is identical. The restricted models are created using the
+\textit{restModel} method, an R is either a vector of characters or a
+list of formulas. There is no need to specify the equation names
+because the coefficient names are unique. The following are two types
+of restrictions, the first being equation by equation and the second
+involving a cross-equation restriction.
+<<>>=
+R1 <- c("theta1=-12*theta2","theta9=0.8", "theta11=0.3")
+R2<- c("theta1=1", "theta6=theta10")
+(rnsmod1 <- restModel(nsmod, R1))
+(rnsmod2 <- restModel(nsmod, R2))
+@
+
\subsection{Genelarized method of moments}\label{sec:sysgmm}
-\subsubsection{A class for GMM weights}\label{sec:smomentmodel-weights}
+\subsubsection{A class for moment weights}\label{sec:smomentmodel-weights}
As for the single equation case, the weighting matrices must have a
particular class in order to work with all model fitting methods. The
@@ -2370,6 +2432,34 @@
coef(rsmod1, theta1)
@
+The method also applies to restricted system of nonlinear equations
+(with and without cross-equation restrictions). It is important to
+provide good starting values to the minimization algorithm if we want
+the method to converge to the global minimum. In the following, a
+vector of 0's is used to get the first-step estimate, but in practice
+it is recommended to find a better strategy. The starting values for
+the second-step estimate is the first-step estimate. It is the ideal
+starting values provided that the first-step method converged.
+
+<<>>=
+### Without cross-equation restrictions
+wObj1 <- evalWeights(rnsmod1, w="ident")
+theta0 <- solveGmm(rnsmod1, wObj1, theta0=rep(0, 8))$theta
+wObj2 <- evalWeights(rnsmod1, theta=theta0)
+theta1 <- solveGmm(rnsmod1, wObj2, theta0=theta0)$theta
+### Verify that the restrictions are correctly imposed:
+printRestrict(rnsmod1)
+coef(rnsmod1, theta1)
+### With cross-equation restrictions
+wObj1 <- evalWeights(rnsmod2, w="ident")
+theta0 <- solveGmm(rnsmod2, wObj1, theta0=rep(0, 9))$theta
+wObj2 <- evalWeights(rnsmod2, theta=theta0)
+theta1 <- solveGmm(rnsmod2, wObj2, theta0=theta0)$theta
+### Verify that the restrictions are correctly imposed:
+printRestrict(rnsmod2)
+coef(rnsmod2, theta1)
+@
+
\subsubsection{The \textit{gmmFit} method for system of equations}\label{sec:smomentmodel-gmmfit}
This is the main algorithm to obtain GMM estimates of systems of
@@ -2451,6 +2541,18 @@
gmmFit(rsmod1)@theta
@
+The following are the estimation of the two above restricted systems
+of nonlinear equations (the vector of 0's is used again as starting
+values because the initial values included in the model object does a
+poor job):
+
+<<>>=
+theta0 <- setCoef(rnsmod1, rep(0,8))
+gmmFit(rnsmod1, theta0=theta0)@theta
+theta0 <- setCoef(rnsmod2, rep(0,9))
+gmmFit(rnsmod2, theta0=theta0)@theta
+@
+
\subsubsection{The \textit{tsls} and \textit{ThreeSLS} methods}
A system of equation can be estimated by 2SLS equation by equation
@@ -2603,7 +2705,33 @@
hypothesisTest(res.u, res2.r, type="LR")
@
+For the nonlinear model, it works in a very similar way. First we
+estimate the unrestricted model and the two restricted ones.
+<<>>=
+R1 <- c("theta1=-12*theta2","theta9=0.8", "theta11=0.3")
+R2<- c("theta1=1", "theta6=theta10")
+rnsmod1 <- restModel(nsmod, R1)
+rnsmod2 <- restModel(nsmod, R2)
+theta0 <- setCoef(nsmod, rep(0,11))
+fit <- gmmFit(nsmod, theta0=theta0)
+theta0 <- setCoef(rnsmod1, rep(0,8))
+rfit1 <- gmmFit(rnsmod1, theta0=theta0)
+theta0 <- setCoef(rnsmod2, rep(0,9))
+rfit2 <- gmmFit(rnsmod2, theta0=theta0)
+@
+
+Then, we test the two restrictions using the different options:
+
+<<>>=
+hypothesisTest(object.u=fit, R=R1)
+hypothesisTest(object.u=fit, object.r=rfit1, type="LR")
+hypothesisTest(object.u=fit, object.r=rfit1, type="LM")
+hypothesisTest(object.u=fit, R=R2)
+hypothesisTest(object.u=fit, object.r=rfit2, type="LR")
+hypothesisTest(object.u=fit, object.r=rfit2, type="LM")
+@
+
\end{itemize}
\subsubsection{Direct estimation with \textit{gmm4}}\label{sec:smomentmodel-gmm4}
@@ -2643,6 +2771,25 @@
res
@
+It is the same for nonlinear systems. Notice that theta0 that needs to
+be provided is for the unrestricted model even if impose
+restriction. \textit{gmm4} first creates the unrestricted model with
+theta0 and use \textit{restModel} after to created the restricted
+model.
+<<>>=
+h <- list(~z1+z2+z3, ~x3+z1+z2+z3+z4, ~x3+x4+z1+z2+z3)
+nlg <- list(Supply=y1~theta0+theta1*x1+theta2*z2,
+ Demand1=y2~alpha0+alpha1*x1+alpha2*x2+alpha3*x3,
+ Demand2=y3~beta0+beta1*x3+beta2*x4+beta3*z1)
+theta0 <- list(c(theta0=0,theta1=0,theta2=0),
+ c(alpha0=0,alpha1=0,alpha2=0, alpha3=0),
+ c(beta0=0,beta1=0,beta2=0,beta3=0))
+fit <- gmm4(nlg, h, theta0,data=simData)
+## the restricted estimation (:
+R2<- c("theta1=1", "alpha1=beta2")
+fit2 <- gmm4(nlg, h, theta0,data=simData, cstLHS=R2)
+@
+
\subsection{Textbooks Applications} \label{sec:smomentmodel-app}
\subsubsection{Greene}
Modified: pkg/momentfit/vignettes/gmmS4.pdf
===================================================================
(Binary files differ)
More information about the Gmm-commits
mailing list