Nothing
spsepgm<-function(formula, data = list(), index = NULL, listw, model = c("within","random"), lag = NULL, spatial.error = NULL, moments = c("initial", "weights", "fullweights"), endog = NULL, instruments = NULL, verbose = FALSE, method = c("w2sls", "b2sls", "g2sls", "ec2sls"), control = list()) {
#source("spregm.R")
effects<-match.arg(model)
moments<-match.arg(moments)
iindex<-index
if (!is.null(index)) {
require(plm)
data <- plm.data(data, index)
}
index <- data[, 1]
tindex <- data[, 2]
cl <- match.call()
if (dim(data)[[1]] != length(index)) stop("Non conformable arguments")
mt <- terms(formula[[1]], data = data)
mf <- lm(formula[[1]], data, na.action = na.fail, method = "model.frame")
y <- model.extract(mf, "response")
x <- model.matrix(mt, mf)
names(index) <- row.names(data)
ind <- index[which(names(index) %in% row.names(x))]
tind <- tindex[which(names(index) %in% row.names(x))]
N <- length(unique(ind))
T <- max(tapply(x[, 1], ind, length))
NT <- length(ind)
balanced <- N * T == NT
if (!balanced)
stop("Estimation method unavailable for unbalanced panels")
index<-iindex
if (!inherits(formula, "list"))
stop("formula should be a list in simultaneous equation models")
eq <- length(formula)
if (!inherits(listw, "list")) {
if (!inherits(listw, c("listw", "matrix")))
stop("listw has to be of class 'listw' or 'matrix' if it is not a list")
else {
w.list <- vector("list", length = eq )
for ( i in 1:eq) w.list[[i]]<-listw
}
}
else{
w.list<-listw
}
n <- dim(data)[[1]]
if(!is.null(lag) & length(lag) != eq) stop("Argument lag incorrectly specified")
if(is.null(lag)) lag<-rep(FALSE, eq)
if(!is.null(endog) && is.null(instruments)) stop("Instruments not specified")
est.res <- vector("list", length = eq)
bigk<-0
expv<-vector("numeric", length=0)
rhos<-vector("numeric", length=0)
namesx<-vector("list", length=0)
for (i in 1:eq) {
est.res[[i]]<-spgm(formula=formula[[i]], data=data, index = index, listw = w.list[[i]],
model = model, lag = lag[[i]], spatial.error = spatial.error[[i]], moments = moments, endog = endog, instruments = instruments,
verbose = verbose, method = method, control = list())
#print(summary(est.res[[i]]))
bigk<-bigk + est.res[[i]]$rhs
expv<-c(expv,length(est.res[[i]]$coefficients))
namesx[[i]]<-names(est.res[[i]]$coefficients)
rhos<-c(rhos,est.res[[i]]$rho[1])
}
Umat <- matrix(,NT,eq)
bigy<-matrix(0,NT*eq,1)
bigx<-matrix(0,NT*eq,bigk)
model.data<-data.frame(matrix(0,NT,0))
counter<-0
for (i in 1:eq) {
Umat[,i]<-est.res[[i]]$residuals
upl<-NT*i
lowl<-upl-NT+1
counter<-counter+est.res[[i]]$rhs
bigy[lowl:upl,]<-est.res[[i]]$coy
bigx[lowl:upl,((counter-est.res[[i]]$rhs + 1):counter)]<-est.res[[i]]$cox
model.data<-data.frame(cbind(model.data,est.res[[i]]$model))
}
#pos<-c(1,10,13,22)
#print(bigx[1:10,])
ind <- rep(seq(1, N),T)
avctp <- function(x) rep(unlist(tapply(x, ind, mean, simplify = TRUE)), T)
Q1Umat <- apply(Umat, 2, avctp)
Q2Umat <- Umat - Q1Umat
Omv<-crossprod(Umat,Q2Umat)/(N*(T-1))
Om1<-crossprod(Umat,Q1Umat)/N
#print(Omv)
#print(Om1)
Omvinv<-invfrac(Omv, -0.5)
Om1inv<-invfrac(Om1, -0.5)
#print(Omvinv)
#print(Om1inv)
JT<-matrix(1/T,T,T)
IN<-Diagonal(N)
Q1<-kronecker(JT, IN)
IT<-Diagonal(T)
FP<-matrix((1-(1/T)),T,T)
Q2<-kronecker(FP, IN)
Omuinv<-kronecker(Om1inv, Q1) + kronecker(Omvinv, Q2)
#print(Omuinv)
finy<-Omuinv %*% bigy
finx<-Omuinv %*% bigx
finx<-as.matrix(finx)
#print(finx[1:10,])
#print(finx)
#finxpfinx<-crossprod(as.matrix(finx))
#print(cor(as.matrix(finxpfinx)))
#finxpfinxinv<-solve(as.matrix(finxpfinx), tol = 2e-20)
#pos<-c(1,10,13,22)
#print(finx[,pos])
modfin<-lm(as.matrix(finy)~as.matrix(finx)-1)
#betaFGLS<-finxpfinxinv %*% crossprod(finx,finy)
betaFGLS<-coefficients(modfin)
#print(betaFGLS)
# fv <- as.vector(finx %*% betaFGLS)
fv<-fitted.values(modfin)
egls <- finy - fv
SGLS <- sum(egls^2)/(N - 1)
# xfpxfNT <- (1/NT) * crossprod(finx)
# PSI <- solve(xfpxfNT, tol=2e-20)
covbeta <- vcov(modfin)
type<- "Seemingly unrelated regressions with spatial error components"
#print(covbeta)
spmod<- list(coefficients= betaFGLS, vcov = covbeta, residuals = as.vector(egls), fitted.values = fv,
sigma2 = SGLS, type = type, rho= rhos, model.data = model.data, call = cl, expv = expv, namesx = namesx, NT = NT, N = N, T = T, eq = eq)
class(spmod)<-"spsugm"
return(spmod)
}
summary.spsugm<-function(object,...){
coeff<-object$coefficients
eq <- object$eq
expv <- object$expv
var<-object$vcov
residuals<-object$residuals
serr<-sqrt(diag(as.matrix(var)))
tval<-coeff/serr
pval<-pnorm(abs(as.matrix(tval)), lower.tail = FALSE) * 2
Xnam<-object$namesx
NT<-object$NT
N<-object$N
T<-object$T
NT<-object$NT
rho<-object$rho
tmp<-1:eq
tmp2<-rep(tmp,expv)
beta<-split(coeff,tmp2)
#print(coeff)
#print(beta)
se<-split(serr,tmp2)
tv<-split(tval,tmp2)
pv<-split(pval,tmp2)
namesx<-split(unlist(object$namesx),tmp2)
tmp3<-rep(tmp,each=NT)
resid<-split(residuals, tmp3)
object$beta<-beta
object$se<-se
object$tv<-tv
object$pv<-pv
object$eq<-eq
object$namesx<-namesx
object$Xnam<-Xnam
object$resid<-resid
class(object)<- c("summary.spsugm","spsugm")
object
}
print.summary.spsugm<-function(x, digits = max(3, getOption("digits") - 2), width = getOption("width"),
...){
cat("\nSeemingly unrelated regressions with spatial error components\n")
cat("\nCall:\n")
print(x$call)
eq <- x$eq
save.digits <- unlist(options(digits = digits))
on.exit(options(digits = save.digits))
namesx <- x$namesx
for (i in 1:eq) {
cat(" \n")
cat(paste("Equation", i, sep = " ", collapse = ""), "\n", sep = "")
cat("\nResiduals:\n")
save.digits <- unlist(options(digits=digits))
on.exit(options(digits=save.digits))
sr <- summary(x$resid[[i]])
srm <- mean(x$resid[[i]])
if (abs(srm) < 1e-10) sr <- sr[c(1:3,5:6)]
print(sr)
tables <- cbind(as.matrix(x$beta[[i]]), as.matrix(x$se[[i]]), as.matrix(x$tv[[i]]), as.matrix(x$pv[[i]]))
dimnames(tables) <- list(namesx[[i]], c("Estimate", "Std.Error", "t value", "Pr(>|t|)"))
# rownames(tables)
if (i == eq) {
legend = TRUE
}
else {
legend = FALSE
}
printCoefmat(tables, digits = digits, signif.stars = TRUE,
signif.legend = legend)
cat(" \n")
if (!is.null(x$rho[i])) {
cat(paste("Spatial autoregressive parameter:",
round(x$rho[i], 4), sep = " ", collapse = ""),
"\n", sep = "")
}
cat("\n _______________________________________________________ \n")
}
invisible(x)
}
invfrac<-function(Rmat, p){
#if(p<0) Rmat<-solve(Rmat)
#print(p)
#print(solve(Rmat))
ee<-eigen(Rmat, symmetric = TRUE)
Vmat<-ee$vectors
Dmat<-diag(ee$values)
Dmatp<-Dmat^abs(p)
Rmatoh<-Vmat %*% Dmatp %*% solve(Vmat)
if(p<0) Rmatoh<-solve(Rmatoh)
return(Rmatoh)
}
`print.spsugm` <-
function(x, digits = max(3, getOption("digits") - 2), ...) {
cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
if (length(coef(x))) {
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2,
quote = FALSE)
} else {
cat("No coefficients\n")
}
## add printing of error variance parameters
cat("\n")
ec <- x$errcomp
if (length(ec)) {
cat("Error covariance parameters:\n")
print.default(format(ec, digits = digits), print.gap = 2,
quote = FALSE)
}
else cat("No error covariance parameters\n")
cat("\n")
invisible(x)
}
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.