Nothing
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
five <- function(g, h, commonCoef = FALSE, data = NULL)
{
res <- sysGmm(g, h, wmatrix = "optimal", vcov = "CondHom", commonCoef=commonCoef, data=data)
attr(res$dat, "sysInfo")$typeDesc <- "Full-Information IV (FIVE)"
res$call <- match.call()
res
}
threeSLS <- function(g, h, commonCoef = FALSE, data = NULL)
{
if (!is(h, "formula"))
if (length(h) != 1)
stop("In 3SLS, h is a single equation since the same instruments are used in each equation")
res <- sysGmm(g, h, vcov = "CondHom", commonCoef=commonCoef, data=data)
attr(res$dat, "sysInfo")$typeDesc <- "3-stage Least Squares"
res$call <- match.call()
res
}
sur <- function(g, commonCoef = FALSE, data = NULL)
{
if (!is.list(g))
stop("g must be list of formulas")
if (length(g) == 1)
stop("For single equation GMM, use the function gmm()")
if (!all(sapply(1:length(g), function(i) is(g[[i]], "formula"))))
stop("g must be a list of formulas")
tm <- lapply(1:length(g), function(i) terms(g[[i]]))
reg <- lapply(1:length(g), function(i) attr(tm[[i]], "term.labels"))
reg <- unique(do.call(c,reg))
h <- paste(reg, collapse = "+")
if (all(sapply(1:length(g), function(i) attr(tm[[i]], "intercept") == 0)))
h <- paste(h, "-1")
h <- as.formula(paste("~", h))
res <- sysGmm(g, h, vcov = "CondHom", commonCoef=commonCoef, data=data)
attr(res$dat, "sysInfo")$typeDesc <- "Seemingly Unrelated Regression (SUR)"
res$call <- match.call()
res
}
randEffect <- function(g, data = NULL)
{
res <- sur(g, commonCoef = TRUE, data = data)
attr(res$dat, "sysInfo")$typeDesc <- "Random Effect Estimator"
res$call <- match.call()
res
}
sysGmm <- function(g, h, wmatrix = c("optimal","ident"),
vcov=c("MDS", "HAC", "CondHom", "TrueFixed"),
kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen", "Tukey-Hanning"),
crit=10e-7,bw = bwAndrews, prewhite = FALSE, ar.method = "ols", approx="AR(1)",tol = 1e-7,
model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, weightsMatrix = NULL,
data, crossEquConst = NULL, commonCoef = FALSE)
{
kernel <- match.arg(kernel)
vcov <- match.arg(vcov)
wmatrix <- match.arg(wmatrix)
TypeGmm = "sysGmm"
if(vcov=="TrueFixed" & is.null(weightsMatrix))
stop("TrueFixed vcov only for fixed weighting matrix")
if(!is.null(weightsMatrix))
wmatrix <- "optimal"
if(missing(data))
data<-NULL
all_args<-list(data = data, g = g, h = h, wmatrix = wmatrix, vcov = vcov, kernel = kernel,
crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method, approx = approx,
weightsMatrix = weightsMatrix, centeredVcov = centeredVcov, tol = tol,
model = model, X = X, Y = Y, call = match.call(), commonCoef=commonCoef, crossEquConst = crossEquConst)
class(all_args)<-TypeGmm
Model_info<-getModel(all_args)
z <- momentEstim(Model_info)
z <- FinRes(z, Model_info)
return(z)
}
.momentFct_Sys <- function(tet, dat)
{
q <- length(dat)
f <- function(i, dat)
{
d <- dat[[i]]
attr(d, "eqConst") <- attr(dat, "eqConst")
attr(d, "ModelType") <- attr(dat, "ModelType")
attr(d,"momentfct") <- attr(dat,"momentfct")
attr(d, "smooth") <- attr(dat, "smooth")
.momentFct(tet[[i]], d)
}
mom <- lapply(1:q, function(i) f(i, dat))
do.call(cbind, mom)
}
.DmomentFct_Sys <- function(tet, dat, pt=NULL)
{
q <- length(dat)
f <- function(i, dat)
{
d <- dat[[i]]
attr(d, "eqConst") <- attr(dat, "eqConst")
attr(d, "ModelType") <- attr(dat, "ModelType")
attr(d,"momentfct") <- attr(dat,"momentfct")
attr(d, "smooth") <- attr(dat, "smooth")
.DmomentFct(tet[[i]], d, pt)
}
dmom <- lapply(1:q, function(i) f(i, dat))
if (attr(dat, "sysInfo")$commonCoef)
do.call(rbind,dmom)
else
.diagMatrix(dmom)
}
.weightFct_Sys<- function(tet, dat, type=c("MDS", "HAC", "CondHom", "ident", "fct", "fixed"))
{
type <- match.arg(type)
if (type == "fixed")
{
w <- attr(dat, "weight")$w
attr(w, "inv") <- FALSE
} else if (type == "ident") {
w <- diag(attr(dat, "q"))
attr(w, "inv") <- FALSE
} else {
if (type == "HAC")
{
gt <- .momentFct_Sys(tet,dat)
if(attr(dat, "weight")$centeredVcov)
gt <- residuals(lm(gt~1))
n <- NROW(gt)
obj <- attr(dat, "weight")
obj$centeredVcov <- FALSE
w <- .myKernHAC(gt, obj)
attr(w, "inv") <- TRUE
}
if (type == "MDS")
{
gt <- .momentFct_Sys(tet,dat)
n <- NROW(gt)
if(attr(dat, "weight")$centeredVcov)
gt <- scale(gt, scale=FALSE)
w <- crossprod(gt)/n
attr(w, "inv") <- TRUE
}
if (type == "CondHom")
{
e <- lapply(1:length(dat), function(i) .residuals(tet[[i]], dat[[i]])$residuals)
e <- do.call(cbind, e)
Sig <- crossprod(scale(e, scale=FALSE))/nrow(e)
Z <- lapply(1:length(dat), function(i) dat[[i]]$x[,(2+dat[[i]]$k):ncol(dat[[i]]$x)])
Z <- do.call(cbind, Z)
w <- crossprod(Z)/nrow(e)
for (i in 1:length(dat))
for (j in 1:length(dat))
{
s1 <- 1+(i-1)*dat[[i]]$nh
e1 <- i*dat[[i]]$nh
s2 <- 1+(j-1)*dat[[j]]$nh
e2 <- j*dat[[j]]$nh
w[s1:e1, s2:e2] <- w[s1:e1, s2:e2]*Sig[i,j]
}
attr(w, "inv") <- TRUE
}
}
return(w)
}
.diagMatrix <- function(xlist, which=NULL)
{
# Create block diagonal matrix from matrices with the same number of rows.
m <- length(xlist)
n <- NROW(xlist[[1]])
l <- sapply(1:m, function(i) dim(as.matrix(xlist[[i]])))
if (!is.null(which))
{
if (any(l[2,1]!=l[2,]))
stop("diagMatrix with which given is for X with the same number of columns")
which <- sort(which)
if (length(which) == l[2,1])
{
dimX <- rowSums(l)
which <- NULL
} else {
dimX <- c(sum(l[1,]), length(which)*m + l[2,1]-length(which))
}
} else {
dimX <- rowSums(l)
}
X <- matrix(0, dimX[1], dimX[2])
if (is.null(which))
{
for (i in 1:m)
{
s1 <- 1 + (i-1)*n
e1 <- n*i
s2 <- 1 + sum(l[2,][-(i:m)])
e2 <- sum(l[2,][1:i])
X[s1:e1, s2:e2] <- xlist[[i]]
}
} else {
q <- match(1:l[2,1],which)
q2 <- which(is.na(q))
wai <- 1
for (j in q2)
{
if (wai < j)
{
xlist2 <- lapply(1:length(xlist), function(i) xlist[[i]][,wai:(j-1), drop=FALSE])
k <- j-wai
Xtmp <- .diagMatrix(xlist2)
X[,wai:(wai+length(xlist)*k-1)] <- Xtmp
wai <- wai + length(xlist)*k
}
X[,wai] <- do.call(c, lapply(1:length(xlist), function(i) xlist[[i]][,j]))
wai <- wai+1
}
if (max(q2) < l[2,1])
{
xlist2 <- lapply(1:length(xlist), function(i) xlist[[i]][,-(1:max(q2))])
Xtmp <- .diagMatrix(xlist2)
X[,wai:ncol(X)] <- Xtmp
}
}
X
}
.getThetaList <- function(tet, dat)
{
neq <- length(dat)
if (attr(dat, "sysInfo")$commonCoef)
{
k <- attr(dat, "k")[[1]]
if (length(tet) != k)
stop("Wrong length of tet")
tet2 <- rep(list(c(tet)), neq)
} else if (!is.null(attr(dat, "sysInfo")$crossEquConst)) {
k <- attr(dat, "k")[[1]]
cst <- attr(dat, "sysInfo")$crossEquConst
nk <- (k-length(cst))*neq + length(cst)
if (nk != length(tet))
stop("Wrong length of tet")
tet2 <- matrix(NA, k, neq)
wai <- 1
for (j in cst)
{
if (wai < j)
{
ind <- wai:(j-1)
tet2[ind,] <- tet[1:(neq*length(ind))]
tet <- tet[-(1:(neq*(j-wai)))]
wai <- wai + length(ind)
}
tet2[wai,] <- tet[1]
tet <- tet[-1]
wai <- wai+1
}
if (max(cst) < k)
tet2[(max(cst)+1):k,] <- tet
attr(dat, "sysInfo")$crossEquConst <- NULL
tet2 <- .getThetaList(c(tet2), dat)
} else {
k2 <- do.call(c, attr(dat, "k"))
if (sum(k2) != length(tet))
stop("Wrong length of tet")
tet2 <- list()
for (j in 1:length(k2))
{
tet2[[j]] <- tet[1:k2[j]]
tet <- tet[-(1:k2[j])]
}
}
tet2
}
.chkPerfectFit <- function(obj)
{
r <- as.matrix(obj$residuals)
f <- as.matrix(obj$fitted.values)
rdf <- obj$df.residual
rss <- colSums(r^2)
mf <- colMeans(f)
mv <- apply(f, 2, var)
resvar <- rss/rdf
resvar < (mf^2 + mv) * 1e-30
}
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.