Nothing
####### All methods with sysGmmModels (and its subclasses) signature
#####################################################################
## print and show
setMethod("print", "sysGmmModels",
function(x, ...)
{
cat("System GMM Model\n")
cat("****************\n")
type <- gsub("s", "", is(x)[1])
cat("Moment type: ", strsplit(type, "G")[[1]][1], "\n",
sep = "")
cat("Covariance matrix: ", x@vcov, sep = "")
if (x@vcov == "HAC") {
option <- x@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", "sysGmmModels", function(object) print(object))
## modelDims
setMethod("modelDims", "slinearGmm",
function(object) {
list(q=object@q, k=object@k, n=object@n, parNames=object@parNames,
momNames=object@momNames, eqnNames=object@eqnNames,
isEndo=object@isEndo)
})
setMethod("modelDims", "snonlinearGmm",
function(object) {
list(k=object@k, q=object@q, n=object@n, parNames=object@parNames,
momNames=object@momNames, theta0=object@theta0,
fRHS=object@fRHS, fLHS=object@fLHS, eqnNames=object@eqnNames,
isEndo=object@isEndo)
})
## Subsetting '['
setMethod("[", c("sysGmmModels", "missing", "list"),
function(x, i, j){
if (length(j) != length(x@q))
stop("j must be a list with a length equals to the number of equations")
spec <- modelDims(x)
x@SUR <- FALSE
if (x@sameMom)
{
chk <- sapply(j[-1], function(ji) identical(ji, j[[1]]))
if (!all(chk))
x@sameMom <- FALSE
}
for (s in 1:length(j))
{
if (length(j[[s]]) > 0)
{
q <- spec$q[s]
if (!all(abs(j[[s]]) %in% (1:q)))
stop("SubMoment must be between 1 and q")
momNames <- spec$momNames[[s]][j[[s]]]
if (length(momNames)<spec$k[s])
{
error <- paste("Equation", s, "is under-identified")
stop(error)
}
if (momNames[1] == "(Intercept)")
{
f <- reformulate(momNames[-1], NULL, TRUE)
} else {
f <- reformulate(momNames, NULL, FALSE)
}
attr(f, ".Environment")<- .GlobalEnv
x@q[s] <- length(momNames)
x@instT[[s]] <- terms(f)
x@momNames[[s]] <- momNames
}
}
x
})
setMethod("[", c("snonlinearGmm", "numeric", "missing"),
function(x, i, j){
i <- unique(as.integer(i))
spec <- modelDims(x)
neqn <- length(spec$k)
if (!all(abs(i) %in% (1:neqn)))
stop("Selected equations out of range")
x@fLHS <- x@fLHS[i]
x@fRHS <- x@fRHS[i]
if (length(x@fLHS) == 0)
stop("Removed too many equations; the model is empty")
x@instT <- x@instT[i]
x@k=x@k[i]
x@q <- x@q[i]
x@parNames <- x@parNames[i]
x@momNames <- x@momNames[i]
x@eqnNames <- x@eqnNames[i]
x@theta0 <- x@theta0[i]
x@varNames <- x@varNames[i]
x@isEndo <- x@isEndo[i]
x@SUR <- FALSE
if (length(x@q) > 1)
return(x)
instF <- model.frame(x@instT[[1]], x@data)
varN <- c(all.vars(x@fLHS[[1]]), all.vars(x@fRHS))
varN <- varN[!(varN%in%names(x@theta0[[1]]))]
modelF <- x@data[,varN, drop=FALSE]
new("nonlinearGmm", instF=instF, modelF=modelF, q=x@q[[1]],
fLHS=x@fLHS[[1]], fRHS=x@fRHS[[1]], theta0=x@theta0[[1]],
k=x@k[[1]], parNames=x@parNames[[1]], momNames=x@momNames[[1]],
vcov=x@vcov, n=spec$n,vcovOptions=x@vcovOptions,
centeredVcov=x@centeredVcov, varNames=x@varNames[[1]],
isEndo=x@isEndo[[1]], survOptions=x@survOptions)
})
setMethod("[", c("slinearGmm", "numeric", "missing"),
function(x, i, j){
i <- unique(as.integer(i))
spec <- modelDims(x)
neqn <- length(spec$k)
if (!all(abs(i) %in% (1:neqn)))
stop("Selected equations out of range")
x@modelT <- x@modelT[i]
if (length(x@modelT) == 0)
stop("Removed too many equations; the model is empty")
x@instT <- x@instT[i]
x@k=x@k[i]
x@q <- x@q[i]
x@parNames <- x@parNames[i]
x@momNames <- x@momNames[i]
x@eqnNames <- x@eqnNames[i]
x@varNames <- x@varNames[i]
x@isEndo <- x@isEndo[i]
x@SUR <- FALSE
if (length(x@q) > 1)
return(x)
instF <- model.frame(x@instT[[1]], x@data)
modelF <- model.frame(x@modelT[[1]], x@data)
new("linearGmm", instF=instF, modelF=modelF, q=x@q[[1]],
k=x@k[[1]], parNames=x@parNames[[1]], momNames=x@momNames[[1]],
vcov=x@vcov, n=spec$n, vcovOptions=x@vcovOptions,
centeredVcov=x@centeredVcov, varNames=x@varNames[[1]],
isEndo=x@isEndo[[1]],survOptions=x@survOptions)
})
setMethod("[", c("sysGmmModels", "numeric", "list"),
function(x, i, j){
x <- x[i]
if (!inherits(x, "sysGmmModels"))
{
if (length(j)>1)
warning("length(j)>1, only the first element used")
x[,j[[1]]]
} else {
x[,j]
}
})
setMethod("[", c("sysGmmModels", "missing", "missing"),
function(x, i, j) x)
## Observation subset
setMethod("subset", "sysGmmModels",
function(x, i) {
x@data <- x@data[i,,drop=FALSE]
if (!is.null(x@vcovOptions$cluster))
{
if (!is.null(dim(x@vcovOptions$cluster)))
x@vcovOptions$cluster <- x@vcovOptions$cluster[i]
else
x@vcovOptions$cluster <- x@vcovOptions$cluster[i,,drop=FALSE]
}
x@n <- nrow(x@data)
x})
### merge
setGeneric("merge")
setMethod("merge", c("linearGmm", "linearGmm"),
function(x, y, ...) {
all <- c(list(x), y, list(...))
cl <- sapply(all, class)
if (any(cl != "linearGmm"))
stop("Can only merge linearGmm with other linearGmm models")
n <- sapply(all, function(s) modelDims(s)$n)
if (any(n[-1L] != n[1L]))
stop("You can only merge models with the same number of observations")
k <- sapply(all, function(s) modelDims(s)$k)
q <- sapply(all, function(s) modelDims(s)$q)
parNames <- lapply(all, function(s) modelDims(s)$parNames)
momNames <- lapply(all, function(s) modelDims(s)$momNames)
varNames <- lapply(all, function(s) s@varNames)
isEndo <- lapply(all, function(s) s@isEndo)
instT <- lapply(all, function(s) terms(s@instF))
modelT <- lapply(all, function(s) terms(s@modelF))
dat <- do.call(cbind, lapply(all, function(s) cbind(s@modelF, s@instF)))
dat <- dat[,!duplicated(colnames(dat))]
eqnNames <- paste("Eqn", 1:length(all), sep="")
new("slinearGmm", data=dat, instT=instT, modelT=modelT,
eqnNames=eqnNames, vcov=x@vcov, vcovOptions=x@vcovOptions,
centeredVcov = x@centeredVcov, k=k, q=q, n=n[1], parNames=parNames,
momNames=momNames, sameMom=FALSE, isEndo=isEndo, varNames=varNames,
SUR=FALSE,survOptions=x@survOptions)
})
setMethod("merge", c("nonlinearGmm", "nonlinearGmm"),
function(x, y, ...) {
all <- c(list(x), y, list(...))
cl <- sapply(all, class)
if (any(cl != "nonlinearGmm"))
stop("Can only merge nonlinearGmm with oter nonlinearGmm models")
n <- sapply(all, function(s) modelDims(s)$n)
if (any(n[-1L] != n[1L]))
stop("You can only merge models with the same number of observations")
fRHS <- lapply(all, function(s) s@fRHS)
fLHS <- lapply(all, function(s) s@fLHS)
k <- sapply(all, function(s) modelDims(s)$k)
q <- sapply(all, function(s) modelDims(s)$q)
parNames <- lapply(all, function(s) modelDims(s)$parNames)
momNames <- lapply(all, function(s) modelDims(s)$momNames)
varNames <- lapply(all, function(s) s@varNames)
isEndo <- lapply(all, function(s) s@isEndo)
instT <- lapply(all, function(s) terms(s@instF))
theta0 <- lapply(all, function(s) s@theta0)
eqnNames <- paste("Eqn", 1:length(all), sep="")
dat <- do.call(cbind, lapply(all, function(s) cbind(s@modelF, s@instF)))
dat <- dat[,!duplicated(colnames(dat))]
new("snonlinearGmm", data=dat, instT=instT,
theta0=theta0,fRHS=fRHS,eqnNames=eqnNames,
fLHS=fLHS, vcov=x@vcov, vcovOptions=x@vcovOptions,
centeredVcov = x@centeredVcov, k=k, q=q,
n=n[1], parNames=parNames, isEndo=isEndo, varNames=varNames,
momNames=momNames, sameMom=FALSE, SUR=FALSE,
survOptions=x@survOptions)
})
setMethod("merge", c("snonlinearGmm", "nonlinearGmm"),
function(x, y, ...) {
all <- c(list(y), list(...))
cl <- sapply(all, class)
if (any(cl != "nonlinearGmm"))
stop("Can only merge nonlinearGmm with oter nonlinearGmm models")
n <- sapply(all, function(s) modelDims(s)$n)
spec <- modelDims(x)
if (any(n != spec$n))
stop("You can only merge models with the same number of observations")
fRHS <- c(spec$fRHS, lapply(all, function(s) modelDims(s)$fRHS))
fLHS <- c(spec$fLHS, lapply(all, function(s) modelDims(s)$fLHS))
k <- c(spec$k, sapply(all, function(s) modelDims(s)$k))
q <- c(spec$q, sapply(all, function(s) modelDims(s)$q))
parNames <- c(spec$parNames, lapply(all, function(s) modelDims(s)$parNames))
momNames <- c(spec$momNames, lapply(all, function(s) modelDims(s)$momNames))
varNames <- c(x@varNames, lapply(all, function(s) s@varNames))
isEndo <- c(x@isEndo, lapply(all, function(s) s@isEndo))
instT <- c(x@instT, lapply(all, function(s) terms(s@instF)))
theta0 <- c(spec$theta0, lapply(all, function(s) modelDims(s)$theta0))
eqNames <- x@eqnNames
eqnNames <- c(eqNames, paste("Eqn",
(length(eqNames)+1):length(fRHS), sep=""))
dat <- do.call(cbind, lapply(all, function(s) cbind(s@modelF, s@instF)))
dat <- dat[,!duplicated(colnames(dat))]
new("snonlinearGmm", data=dat, instT=instT,
theta0=theta0,fRHS=fRHS,eqnNames=eqnNames,
fLHS=fLHS, vcov=x@vcov,vcovOptions=x@vcovOptions,
centeredVcov = x@centeredVcov, k=k, q=q,
n=n[1], parNames=parNames, isEndo=isEndo, varNames=varNames,
momNames=momNames, sameMom=FALSE, SUR=FALSE,
survOptions=x@survOptions)
})
setMethod("merge", c("slinearGmm", "linearGmm"),
function(x, y, ...) {
all <- c(list(y), list(...))
cl <- sapply(all, class)
if (any(cl != "linearGmm"))
stop("Can only merge linearGmm with other linearGmm models")
n <- sapply(all, function(s) modelDims(s)$n)
spec <- modelDims(x)
if (any(n != spec$n))
stop("You can only merge models with the same number of observations")
k <- c(spec$k, sapply(all, function(s) modelDims(s)$k))
q <- c(spec$q, sapply(all, function(s) modelDims(s)$q))
parNames <- c(spec$parNames, lapply(all, function(s) modelDims(s)$parNames))
momNames <- c(spec$momNames, lapply(all, function(s) modelDims(s)$momNames))
varNames <- c(x@varNames, lapply(all, function(s) s@varNames))
isEndo <- c(x@isEndo, lapply(all, function(s) s@isEndo))
instT <- c(x@instT, lapply(all, function(s) terms(s@instF)))
modelT <- c(x@modelT, lapply(all, function(s) terms(s@modelF)))
dat <- do.call(cbind, lapply(all, function(s) cbind(s@modelF, s@instF)))
dat <- dat[,!duplicated(colnames(dat))]
eqNames <- x@eqnNames
eqnNames <- c(eqNames, paste("Eqn",
(length(eqNames)+1):length(instT), sep=""))
new("slinearGmm", data=dat, instT=instT, modelT=modelT,
eqnNames=eqnNames, vcov=x@vcov, vcovOptions=x@vcovOptions,
centeredVcov = x@centeredVcov, k=k, q=q, n=n[1], parNames=parNames,
momNames=momNames, sameMom=FALSE, isEndo=isEndo, varNames=varNames,
SUR=FALSE, survOptions=x@survOptions)
})
## residuals
setMethod("residuals", "sysGmmModels",
function(object, theta) {
neqn <- length(object@eqnNames)
r <- sapply(1:neqn, function(i) residuals(object[i], theta[[i]]))
colnames(r) <- object@eqnNames
r
})
## Dresiduals
setMethod("Dresiduals", "sysGmmModels",
function(object, theta) {
neqn <- length(object@eqnNames)
if (missing(theta))
r <- sapply(1:neqn, function(i) Dresiduals(object[i]))
else
r <- sapply(1:neqn, function(i) Dresiduals(object[i], theta[[i]]))
names(r) <- object@eqnNames
r
})
## evalMoment
setMethod("evalMoment", "sysGmmModels",
function(object, theta) {
neqn <- length(object@eqnNames)
gt <- lapply(1:neqn, function(i) evalMoment(object[i], theta[[i]]))
names(gt) <- object@eqnNames
gt
})
## evalDmoments
setMethod("evalDMoment", "sysGmmModels",
function(object, theta) {
neqn <- length(object@eqnNames)
if (missing(theta))
dgt <- lapply(1:neqn, function(i) evalDMoment(object[i]))
else
dgt <- lapply(1:neqn, function(i) evalDMoment(object[i], theta[[i]]))
names(dgt) <- object@eqnNames
dgt
})
## model.matrix
setMethod("model.matrix", "slinearGmm",
function(object, type = c("regressors", "instruments")) {
type <- match.arg(type)
mm <- lapply(1:length(object@eqnNames), function(i)
model.matrix(object[i], type))
names(mm) <- object@eqnNames
mm
})
setMethod("model.matrix", "snonlinearGmm",
function(object, type = c("regressors", "instruments"))
{
type <- match.arg(type)
if (type == "regressors")
stop("no model.matrix of type regressors for nonlinear Gmm. set type to 'instruments' to get the matrix of instruments")
mm <- lapply(1:length(object@eqnNames), function(i)
model.matrix(object[i], type))
names(mm) <- object@eqnNames
mm
})
## evalWeights
### The following multiplies each block matrix of ZZ, when the dimensions
### are defined by dimr and dimc (dim row and dim column) by each element
### of Sigma. Sigma must be length(dimr) by length(dimc)
.SigmaZZ <- function(ZZ, Sigma, dimr, dimc=NULL, lowerTri=FALSE, isSym=TRUE)
{
r1 <- 1
c1 <- 1
if (is.null(dimc))
dimc <- dimr
for (i in 1:length(dimr))
{
r2 <- sum(dimr[1:i])
start <- ifelse(isSym, i, 1)
for (j in start:length(dimc))
{
c2 <- sum(dimc[1:j])
ZZ[r1:r2, c1:c2] <- ZZ[r1:r2, c1:c2]*Sigma[i,j]
c1 <- c1+dimc[j]
}
r1 <- r1+dimr[i]
c1 <- ifelse(isSym, sum(dimc[1:i])+1, 1)
}
if (lowerTri && isSym)
ZZ[lower.tri(ZZ)] <- t(ZZ)[lower.tri(ZZ)]
ZZ
}
setMethod("evalWeights", "sysGmmModels",
function(object, theta = NULL, w="optimal", wObj=NULL)
{
spec <- modelDims(object)
sameMom <- object@sameMom
n <- object@n
neqn <- length(spec$eqnNames)
if (is.matrix(w))
{
if (!all(dim(w) == sum(spec$q)))
stop("The weights matrix has the wrong dimension")
}
if (is.matrix(w) || (w == "ident"))
{
return(new("sysGmmWeights", type="weights",momNames=object@momNames,
wSpec=list(), w=w, Sigma=NULL, sameMom=sameMom,
eqnNames=object@eqnNames))
}
if (w != "optimal")
stop("w is either 'ident', 'optimal' or a matrix")
if (object@vcov == "iid")
{
e <- residuals(object, theta)
type <- "iid"
Sigma <- chol(crossprod(e)/n)
if (!is.null(wObj))
{
if (wObj@type != "iid")
stop("wObj must come from a model with iid errors")
if (ncol(wObj@w$qr) != spec$q[1])
stop("The qr decomposition has the wrong dimension")
w <- wObj@w
} else {
Z <- model.matrix(object, type="instruments")
if (sameMom)
{
w <- qr(Z[[1]]/sqrt(n))
} else {
w <- crossprod(do.call(cbind,Z))/n
}
}
} else if (object@vcov == "MDS") {
type <- "MDS"
gt <- evalMoment(object, theta)
w <- qr(do.call(cbind, gt)/sqrt(n))
Sigma <- NULL
} else {
stop("Only identity, iid and MDS if allowed for now")
}
return(new("sysGmmWeights", type=type,momNames=object@momNames,
wSpec=list(), w=w, Sigma=Sigma, sameMom=sameMom,
eqnNames=object@eqnNames))
})
## evalObjective
setMethod("evalObjective", signature("sysGmmModels", "list", "sysGmmWeights"),
function(object, theta, wObj, ...)
{
gt <- evalMoment(object, theta)
gt <- lapply(gt, function(g) colMeans(g))
gt <- do.call("c", gt)
n <- object@n
obj <- quadra(wObj, gt)
n*obj
})
## modelResponse
setMethod("modelResponse", signature("slinearGmm"),
function(object)
{
neqn <- length(object@eqnNames)
Y <- lapply(1:neqn, function(i) modelResponse(object[i]))
Y
})
## solveGMM
.GListToMat <- function(G)
{
dimG <- sapply(G, dim)
Gmat <- matrix(0, sum(dimG[1,]), sum(dimG[2,]))
r1 <- 1
c1 <- 1
for (i in 1:length(G)) {
r2 <- sum(dimG[1,1:i])
c2 <- sum(dimG[2,1:i])
Gmat[r1:r2, c1:c2] <- as.matrix(G[[i]])
r1 <- r1 + dimG[1,i]
c1 <- sum(dimG[2,1:i]) + 1
}
Gmat
}
.tetReshape <- function(theta, eqnNames, parNames)
{
if (is.list(theta))
{
theta2 <- do.call("c", theta)
k <- sapply(parNames, length)
tn <- paste(rep(eqnNames, k), ".", do.call("c", parNames),
sep = "")
names(theta2) <- tn
} else {
k <- cumsum(sapply(parNames, length))
names(theta) <- do.call("c", parNames)
theta2 <- list(theta[1:k[1]])
if (length(k)>1)
theta2[2:length(k)] <- lapply(2:length(k),
function(i) theta[(k[i-1]+1):k[i]])
names(theta2) <- eqnNames
}
theta2
}
setMethod("solveGmm", c("slinearGmm", "sysGmmWeights"),
function(object, wObj, theta0 = NULL) {
if (wObj@type=="iid" && object@sameMom)
return(ThreeSLS(object, Sigma=wObj@Sigma, qrZ=wObj@w, coefOnly=TRUE))
spec <- modelDims(object)
Y <- modelResponse(object)
Z <- model.matrix(object, type="instruments")
Syz <- lapply(1:length(Y),
function(i) colMeans(Y[[i]]*Z[[i]]))
Syz <- do.call("c", Syz)
G <- evalDMoment(object)
G <- .GListToMat(G)
T1 <- quadra(wObj, G)
T2 <- quadra(wObj, G, Syz)
theta <- -solve(T1, T2)
theta <- .tetReshape(theta, object@eqnNames, object@parNames)
list(theta=theta, convergence=NULL)
})
setMethod("solveGmm", signature("snonlinearGmm", "sysGmmWeights"),
function (object, wObj, theta0 = NULL, ...)
{
if (is.null(theta0))
theta0 <- modelDims(object)$theta0
if (!is.list(theta0))
stop("theta0 must be a list of vectors")
g <- function(theta, wObj, object){
spec <- modelDims(object)
theta <- .tetReshape(theta, object@eqnNames, spec$parNames)
evalObjective(object, theta, wObj)
}
dg <- function(theta, wObj, object) {
spec <- modelDims(object)
theta <- .tetReshape(theta, object@eqnNames, spec$parNames)
gt <- evalMoment(object, theta)
gt <- do.call(cbind, gt)
n <- nrow(gt)
gt <- colMeans(gt)
G <- evalDMoment(object, theta)
G <- .GListToMat(G)
obj <- 2 * n * quadra(wObj, G, gt)
obj
}
spec <- modelDims(object)
theta0 <- .tetReshape(theta0, object@eqnNames, spec$parNames)
res <- optim(par = theta0, fn = g, gr = dg, method = "BFGS",
object = object, wObj = wObj, ...)
theta <- .tetReshape(res$par, object@eqnNames, object@parNames)
list(theta = theta, convergence = res$convergence)
})
## momentVcov
setMethod("momentVcov", signature("sysGmmModels"),
function(object, theta, ...){
spec <- modelDims(object)
q <- spec$q
if (object@vcov == "MDS")
{
gt <- evalMoment(object, theta)
gt <- do.call(cbind, gt)
if (object@centeredVcov)
gt <- scale(gt, scale=FALSE)
w <- crossprod(gt)/nrow(gt)
} else if (object@vcov == "iid") {
e <- residuals(object, theta)
Sigma <- crossprod(e)/nrow(e)
Z <- model.matrix(object, "instrument")
if (object@sameMom)
{
w <- kronecker(Sigma, crossprod(Z[[1]])/nrow(e))
} else {
Z <- crossprod(Z)/nrow(e)
w <- .SigmaZZ(Z, Sigma, q)
}
} else {
stop("not yet implemented for HAC")
}
wn <- paste(rep(spec$eqnNames, q), ".", do.call("c", spec$momNames),
sep = "")
dimnames(w) <- list(wn,wn)
w
})
### tsls
setMethod("tsls", "slinearGmm",
function(model)
{
Call <- try(match.call(call=sys.call(sys.parent())), silent=TRUE)
if (inherits(Call,"try-error"))
Call <- NULL
neqn <- length(model@eqnNames)
res <- lapply(1:neqn, function(i) tsls(model[i]))
w <- lapply(1:neqn, function(i) quadra(res[[i]]@wObj))
w <- .GListToMat(w)
wObj <- evalWeights(model, w=w)
theta <- lapply(res, coef)
names(theta) <- model@eqnNames
new("stsls", theta=theta, convergence=NULL, convIter=NULL,
call=Call, type="tsls", wObj=wObj, niter=1L,
efficientGmm=FALSE, model=model)
})
### 3SLS and SUR
setGeneric("ThreeSLS", function(model, ...) standardGeneric("ThreeSLS"))
setMethod("ThreeSLS", "slinearGmm",
function(model, coefOnly=FALSE, qrZ=NULL, Sigma=NULL) {
Call <- try(match.call(call=sys.call(sys.parent())), silent=TRUE)
if (inherits(Call,"try-error"))
Call <- NULL
if (!inherits(model, "slinearGmm"))
stop("3SLS is for slinearGmm classes")
if (!model@sameMom)
stop("For 3SLS, the instruments must be the same in each equation")
efficientGmm <- model@vcov == "iid"
spec <- modelDims(model)
n <- spec$n
neqn <- length(model@eqnNames)
x <- model.matrix(model)
y <- modelResponse(model)
if (is.null(qrZ))
{
z <- model.matrix(model[1], "instruments")
qrZ <- qr(z/sqrt(n))
} else {
if (!inherits(qrZ,"qr"))
stop("qrZ must be the qr decomposition of Z")
if (ncol(qrZ$qr) != length(spec$momNames[[1]]))
stop("The qr decomposition has the wrong dimension")
}
if (is.null(Sigma))
{
theta <- lapply(1:neqn, function(i) lm.fit(qr.fitted(qrZ, x[[i]]),
y[[i]])$coefficients)
e <- residuals(model, theta)
Sigma <- chol(crossprod(e)/n)
} else {
if (!all(Sigma[lower.tri(Sigma)] == 0))
stop("Sigma must be a cholesky and therefore upper tiangular")
}
if (model@SUR)
{
xhat <- do.call(cbind, x)
type <- "SUR"
} else {
xhat <- qr.fitted(qrZ, do.call(cbind, x))
type <- "3SLS"
}
iSigma <- chol2inv(Sigma)
y <- do.call(cbind,y)
A <- crossprod(xhat)
A <- .SigmaZZ(A, iSigma, spec$k, spec$k, TRUE)
C <- crossprod(xhat, y)
C <- .SigmaZZ(C, iSigma, spec$k, rep(1, neqn), FALSE, FALSE)
C <- rowSums(C)
theta <- .tetReshape(solve(A, C), model@eqnNames, spec$parNames)
if (coefOnly)
return(list(theta=theta, convergence=NULL))
wObj <- new("sysGmmWeights", w=qrZ, Sigma=Sigma, type="iid",
momNames=spec$momNames,
wSpec=list(), sameMom=TRUE, eqnNames=model@eqnNames)
new("sgmmfit", theta=theta, convergence=NULL,
convIter=rep(NULL, neqn), call=Call, type=type, wObj=wObj,
niter=2L, efficientGmm=efficientGmm, model=model)
})
## modelFit
setMethod("modelFit", signature("sysGmmModels"), valueClass="sgmmfit",
function(model, type=c("twostep", "iter","cue", "onestep"),
itertol=1e-7, initW=c("ident", "tsls", "EbyE"), weights="optimal",
itermaxit=100, efficientWeights=FALSE, theta0=NULL,
EbyE=FALSE, ...)
{
Call <- try(match.call(call=sys.call(sys.parent())), silent=TRUE)
if (inherits(Call,"try-error"))
Call <- NULL
chk <- validObject(model)
type <- match.arg(type)
initW <- match.arg(initW)
i <- 1L
chk <- validObject(model, TRUE)
if (!chk)
stop("model is not a valid gmmModels object")
if (is.character(weights) && !(weights%in%c("optimal","ident")))
stop("weights is a matrix or one of 'optimal' or 'ident'")
spec <- modelDims(model)
if (all(spec$q==spec$k))
{
weights <- "ident"
type <- "onestep"
EbyE <- TRUE
} else if (type == "onestep" && !is.matrix(weights)) {
weights <- "ident"
EbyE <- TRUE
} else if (is.matrix(weights) || inherits(weights,"sysGmmWeights")) {
type <- "onestep"
EbyE <- FALSE
} else if (weights == "ident") {
type <- "onestep"
EbyE <- TRUE
}
if (EbyE)
{
neqn <- length(model@eqnNames)
res <- lapply(1:neqn, function(i)
modelFit(model[i], type=type, weights=weights,itertol=itertol,
initW=initW, itermaxit=itermaxit,
efficientWeights=efficientWeights, theta0=theta0, ...))
theta <- lapply(res, coef)
convergence <- sapply(res, function(r) r@convergence)
if (is.list(convergence))
convergence <- do.call("c", convergence)
convIter <- sapply(res, function(r) r@convIter)
niter <- sapply(res, function(r) r@niter)
if (is.list(convIter))
convIter <- do.call("c", convIter)
type <- paste("EBE", type, sep="")
efficientGmm <- FALSE
if (is.character(weights) && weights=="ident")
{
wObj <- evalWeights(model, NULL, "ident")
} else {
wObj <- lapply(res, function(r) quadra(r@wObj))
wObj <- .GListToMat(wObj)
wObj <- evalWeights(model, w=wObj)
}
ans <- new("sgmmfit", theta=theta, convergence=convergence,
convIter=convIter, call=Call, type=type, wObj=wObj,
niter=niter, efficientGmm=efficientGmm, model=model)
return(ans)
}
if (type == "onestep")
{
if (inherits(weights,"sysGmmWeights"))
wObj <- weights
else
wObj <- evalWeights(model, w=weights)
res <- solveGmm(model, wObj, theta0, ...)
convergence <- res$convergence
efficientGmm <- efficientWeights
ans <- new("sgmmfit", theta=res$theta,
convergence=convergence, convIter=NULL, type=type,
wObj=wObj, model=model, call=Call, niter=i,
efficientGmm=efficientGmm)
return(ans)
}
if (type == "twostep")
{
itermaxit <- 1
if (inherits(model, "slinearGmm"))
{
if (model@vcov=="iid" && !model@sameMom && !model@SUR)
type <- "FIVE"
if (initW=="tsls" && model@vcov=="iid" && model@sameMom)
{
obj <- ThreeSLS(model)
obj@call <- Call
}
}
}
if (initW=="tsls")
{
theta0 <- try(coef(tsls(model)), silent=TRUE)
if (inherits(theta0,"try-error"))
stop("Cannot get the initial weights using 2SLS")
} else if (initW == "EbyE") {
neqn <- length(model@eqnNames)
res <- lapply(1:neqn, function(i)
modelFit(model[i], type=type, weights=weights,itertol=itertol,
itermaxit=itermaxit,
efficientWeights=efficientWeights, theta0=theta0, ...))
theta0 <- lapply(res, coef)
} else {
wObj <- evalWeights(model, NULL, "ident")
theta0 <- solveGmm(model, wObj, theta0, ...)$theta
}
bw <- model@vcovOptions$bw
if (type != "cue")
{
while(TRUE)
{
if (i>1 && model@vcov=="iid")
wObj0 <- wObj
else
wObj0 <- NULL
wObj <- evalWeights(model, theta0, "optimal", wObj0)
if (model@vcov=="HAC" && is.character(model@bw))
model@vcovOptions$bw <- wObj@wSpec$bw
res <- solveGmm(model, wObj, theta0, ...)
theta1 <- res$theta
convergence <- res$convergence
tet0 <- do.call("c", theta0)
dif1 <- do.call("c", theta1)-tet0
crit <- sqrt(sum(dif1^2))/(1+sqrt(sum(tet0^2)))
if (crit < itertol & type=="iter")
{
convIter <- 0
break
}
i <- i + 1L
theta0 <- theta1
if (i>itermaxit)
{
if (type %in% c("twostep", "FIVE"))
convIter <- NULL
else
convIter <- 1
break
}
}
} else {
convIter <- NULL
if (model@vcov=="HAC" && is.character(bw))
{
w <- momentVcov(model, theta0)
model@vcovOptions$bw <- attr(w, "Spec")$bw
}
if (model@vcov == "iid")
wObj0 <- evalWeights(model, theta0)
else
wObj0 <- NULL
obj <- function(theta, model, wObj0, spec)
{
theta <- .tetReshape(theta, model@eqnNames,
spec$parNames)
wObj <- evalWeights(model, theta, "optimal", wObj0)
evalObjective(model, theta, wObj)
}
res <- optim(do.call("c",theta0), obj, model=model, wObj0=wObj0,
spec=spec, ...)
theta1 <- .tetReshape(res$par, model@eqnNames,spec$parNames)
convergence <- res$convergence
wObj <- evalWeights(model, theta1, "optimal", wObj0)
}
model@vcovOptions$bw <- bw
new("sgmmfit", theta=theta1, convergence=convergence,
convIter=convIter, call=Call, type=type, wObj=wObj,
niter=i, efficientGmm=TRUE, model=model)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.