[Splm-commits] r150 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 21 22:17:13 CET 2013
Author: gpiras
Date: 2013-03-21 22:17:13 +0100 (Thu, 21 Mar 2013)
New Revision: 150
Modified:
pkg/R/likelihoodsFE.R
pkg/R/spfeml.R
pkg/R/spml.R
pkg/man/spml.Rd
Log:
updated spfeml
Modified: pkg/R/likelihoodsFE.R
===================================================================
--- pkg/R/likelihoodsFE.R 2013-03-21 15:03:55 UTC (rev 149)
+++ pkg/R/likelihoodsFE.R 2013-03-21 21:17:13 UTC (rev 150)
@@ -22,7 +22,7 @@
}
-splaglm<-function(env, zero.policy = zero.policy, interval = interval){
+splaglm<-function(env, zero.policy = zero.policy, interval = interval, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method){
xt <- get("xt", envir = env)
yt <- get("yt", envir = env)
@@ -95,7 +95,7 @@
rho.se <- fdHess[2, 2]
sig.se <- fdHess[1, 1]
- rest.se <- vcov(lm.target)
+ rest.se <- vcov(lm.lag)
}
@@ -195,7 +195,7 @@
-sperrorlm<-function(env, zero.policy = zero.policy, interval = interval){
+sperrorlm<-function(env, zero.policy = zero.policy, interval = interval, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess){
xt <- get("xt", envir = env)
yt <- get("yt", envir = env)
@@ -342,7 +342,7 @@
}
-spsararlm<-function(env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve= 1e-10){
+spsararlm<-function(env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method){
xt <- get("xt", envir = env)
yt <- get("yt", envir = env)
@@ -448,10 +448,11 @@
names(rho) <- "rho"
lambda <- optres$par[1]
names(lambda) <- "lambda"
- nm <- paste(method, "opt", sep = "_")
- timings[[nm]] <- proc.time() - .ptime_start
- .ptime_start <- proc.time()
+ # timings[["coefs"]] <- proc.time() - .ptime_start
+ # .ptime_start <- proc.time()
+ # assign("first_time", TRUE, envir = env)
+
lm.target <- lm(I(yt - lambda * wyt - rho * w2yt + rho * lambda *
w2wyt) ~ I(xt - rho * wxt) - 1)
@@ -464,17 +465,22 @@
betas <- coefficients(lm.target)
names(betas) <- colnames(xt)
coefs <- c(rho, lambda, betas)
- coefsl <- c(s2, rho, lambda, betas)
+ # coefsl <- c(rho, lambda, betas)
###Add the vc matrix exact
-
- fd <- fdHess(coefsl, f_sacpanel_hess, env)
+if(Hess){
+ fd <- fdHess(coefs, f_sacpanel_hess, env)
mat <- fd$Hessian
fdHess<- solve(-(mat), tol.solve = tol.solve)
rownames(fdHess) <- colnames(fdHess) <- c("lambda", "lambda",colnames(xt))
rho.se <- fdHess[2, 2]
lambda.se <- fdHess[1, 1]
rest.se <- vcov(lm.target)
+ }
+
+ else{
+ stop("Asymptotic VC matrix not yet implemented for model SARAR")
+ }
return<-list(coeff=betas,lambda=lambda, rho = rho, s2=s2, asyvar1 = rest.se, lambda.se = lambda.se, rho.se = rho.se)
@@ -484,13 +490,15 @@
{
T<-get("T", envir = env)
NT<-get("NT", envir = env)
- s2 <- coefs[1]
- rho <- coefs[2]
- lambda <- coefs[3]
- beta <- coefs[-(1:3)]
- # SSE <- sar_sac_hess_sse_panel(rho, lambda, beta, env)
+
+ rho <- coefs[1]
+ lambda <- coefs[2]
+ beta <- coefs[-(1:2)]
+ # SSE <- sar_sac_hess_sse_panel(rho, lambda, beta, env)
+ SSE <- sar_sac_hess_sse_panel(rho, lambda, beta, env)
n <- NT/T
- SSE<- s2 *n
+ # SSE<- s2 *n
+ s2<- SSE / n
ldet1 <- do_ldet(rho, env, which = 1)
ldet2 <- do_ldet(lambda, env, which = 2)
ret <- T * ldet1 + T * ldet2 - ((n*T/2) * log(2 * pi)) - (n*T/2) * log(s2) - (1/(2 * s2)) * SSE
@@ -501,16 +509,16 @@
ret
}
-# sar_sac_hess_sse_panel <- function (rho, lambda, beta, env)
-# {
- # yl <- get("yt", envir = env) - lambda * get("wyt", envir = env) -
- # rho * get("w2yt", envir = env) + rho * lambda * get("w2wyt",
- # envir = env)
- # xl <- get("xt", envir = env) - rho * get("wxt", envir = env)
- # res <- yl - (xl %*% beta)
- # SSE <- c(crossprod(res))
- # SSE
-# }
+sar_sac_hess_sse_panel <- function (rho, lambda, beta, env)
+{
+ yl <- get("yt", envir = env) - lambda * get("wyt", envir = env) -
+ rho * get("w2yt", envir = env) + rho * lambda * get("w2wyt",
+ envir = env)
+ xl <- get("xt", envir = env) - rho * get("wxt", envir = env)
+ res <- yl - (xl %*% beta)
+ SSE <- c(crossprod(res))
+ SSE
+}
Modified: pkg/R/spfeml.R
===================================================================
--- pkg/R/spfeml.R 2013-03-21 15:03:55 UTC (rev 149)
+++ pkg/R/spfeml.R 2013-03-21 21:17:13 UTC (rev 150)
@@ -1,4 +1,4 @@
-spfeml<-function(formula, data=list(), index=NULL, listw, listw2 = NULL, na.action, model = c("lag","error", "sarar"),effects = c('pooled','spfe','tpfe','sptpfe'), method="eigen", quiet = TRUE, zero.policy = NULL, interval1 = NULL, interval2 = NULL, trs1 = NULL, trs2 = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE, llprof = NULL, cl = NULL, Hess = FALSE, LeeYu = FALSE, ...){
+spfeml<-function(formula, data=list(), index=NULL, listw, listw2 = NULL, na.action, model = c("lag","error", "sarar"),effects = c('pooled','spfe','tpfe','sptpfe'), method="eigen", quiet = TRUE, zero.policy = NULL, interval1 = NULL, interval2 = NULL, trs1 = NULL, trs2 = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE, llprof = NULL, cl = NULL, Hess = TRUE, LeeYu = FALSE, ...){
timings <- list()
@@ -19,7 +19,7 @@
warning("unknown names in control: ", paste(noNms, collapse = ", "))
if (is.null(quiet))
- quiet <- !get("verbose", envir = .spdepOptions)
+ quiet <- !get("verbose", envir = spdep:::.spdepOptions)
stopifnot(is.logical(quiet))
if (is.null(zero.policy))
@@ -95,6 +95,13 @@
## det. total number of obs. (robust vs. unbalanced panels)
NT<-length(ind)
+
+ mt<-terms(formula,data=data)
+ mf<-lm(formula,data,method="model.frame")#,na.action=na.fail
+
+ na.act<-attr(mf,'na.action')
+
+
##checks on listw
if(is.matrix(listw)) {
if(dim(listw)[[1]] !=N ) stop("Non conformable spatial weights")
@@ -106,7 +113,7 @@
can.sim <- FALSE
if (listw$style %in% c("W", "S"))
- can.sim <- can.be.simmed(listw)
+ can.sim <- spdep:::can.be.simmed(listw)
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
@@ -135,7 +142,7 @@
can.sim2 <- FALSE
if (listw2$style %in% c("W", "S"))
- can.sim2 <- can.be.simmed(listw2)
+ can.sim2 <- spdep:::can.be.simmed(listw2)
if (!is.null(na.act)) {
subset <- !(1:length(listw2$neighbours) %in% na.act)
listw2 <- subset(listw2, subset, zero.policy = zero.policy)
@@ -179,10 +186,6 @@
if(!balanced) stop("Estimation method unavailable for unbalanced panels")
- mt<-terms(formula,data=data)
- mf<-lm(formula,data,method="model.frame")#,na.action=na.fail
- na.act<-attr(mf,'na.action')
-
indic<-seq(1,T)
inde<-as.numeric(rep(indic,each=N)) ####takes the first n observations
indic1<-seq(1,N)
@@ -206,34 +209,35 @@
if(LeeYu){
-if (effects=="spfe" | effects=="sptpfe"){
-IT <- Diagonal(T)
-IN <- Diagonal(n)
-JT <- matrix(1,T,T)
-Jbar <- 1/T * JT
-Qmat <-IT - Jbar
-vec <- eigen(Qmat)
-Fmat <- vec$vectors[,vec$values==1L]
-Ftm <- kronecker(t(Fmat), IN)
-}
+stop("Lee and Yu correction not yet implemented")
+# if (effects=="spfe" | effects=="sptpfe"){
+# IT <- Diagonal(T)
+# IN <- Diagonal(n)
+# JT <- matrix(1,T,T)
+# Jbar <- 1/T * JT
+# Qmat <-IT - Jbar
+# vec <- eigen(Qmat)
+# Fmat <- vec$vectors[,vec$values==1L]
+# Ftm <- kronecker(t(Fmat), IN)
+# }
-if (effects=="tpfe" | effects=="sptpfe"){
-IT <- Diagonal(T)
-IN <- Diagonal(n)
-JT <- matrix(1,T,T)
-Jbar <- 1/T * JT
-Qmat <-IT - Jbar
-vec <- eigen(Qmat)
-Fmat <- matrix(vec$vectors[,vec$values==1L], T, T-1)
-Ftm <- kronecker(t(Fmat), IN)
-iotan <- matrix(1,n,1)
-Jnbar <-1/n * iotan %*% t(iotan)
-Qmat1 <- IN - Jnbar
-vec1 <- eigen(Qmat1)
-Fmat1 <- matrix(vec1$vectors[,vec1$values==1L], n, n-1)
-FFmat<- kronecker(t(Fmat), t(Fmat1))
-}
+# if (effects=="tpfe" | effects=="sptpfe"){
+# IT <- Diagonal(T)
+# IN <- Diagonal(n)
+# JT <- matrix(1,T,T)
+# Jbar <- 1/T * JT
+# Qmat <-IT - Jbar
+# vec <- eigen(Qmat)
+# Fmat <- matrix(vec$vectors[,vec$values==1L], T, T-1)
+# Ftm <- kronecker(t(Fmat), IN)
+# iotan <- matrix(1,n,1)
+# Jnbar <-1/n * iotan %*% t(iotan)
+# Qmat1 <- IN - Jnbar
+# vec1 <- eigen(Qmat1)
+# Fmat1 <- matrix(vec1$vectors[,vec1$values==1L], n, n-1)
+# FFmat<- kronecker(t(Fmat), t(Fmat1))
+# }
@@ -381,29 +385,29 @@
cat(paste("\nSpatial autoregressive error model\n", "Jacobian calculated using "))
if(model == "lag"){
- interval1 <- jacobianSetup(method, env, con, pre_eig = con$pre_eig, trs = trs, interval = interval1)
+ interval1 <- spdep:::jacobianSetup(method, env, con, pre_eig = con$pre_eig, trs = trs1, interval = interval1)
assign("interval1", interval1, envir = env)
nm <- paste(method, "set_up", sep = "_")
timings[[nm]] <- proc.time() - .ptime_start
.ptime_start <- proc.time()
- RES<- splaglm(env = env, zero.policy = zero.policy, interval = interval1)
+ RES<- splaglm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess)
res.eff<-felag(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method =method, rho=RES$rho, legacy = legacy, zero.policy = zero.policy)
}
if(model == "sarar"){
- interval1 <- jacobianSetup(method, env, con, pre_eig = con$pre_eig1, trs = trs1, interval = interval1, which = 1)
+ interval1 <- spdep:::jacobianSetup(method, env, con, pre_eig = con$pre_eig1, trs = trs1, interval = interval1, which = 1)
assign("interval1", interval1, envir = env)
- interval2 <- jacobianSetup(method, env, con, pre_eig = con$pre_eig2, trs = trs2, interval = interval2, which = 2)
+ interval2 <- spdep:::jacobianSetup(method, env, con, pre_eig = con$pre_eig2, trs = trs2, interval = interval2, which = 2)
assign("interval2", interval2, envir = env)
nm <- paste(method, "set_up", sep = "_")
timings[[nm]] <- proc.time() - .ptime_start
.ptime_start <- proc.time()
- RES<- spsararlm(env = env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve = tol.solve)
+ RES<- spsararlm(env = env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve = tol.solve, Hess = Hess)
res.eff<-felag(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method = method, rho = RES$lambda, legacy = legacy, zero.policy = zero.policy)
@@ -414,24 +418,17 @@
if (model=='error'){
- interval1 <- jacobianSetup(method, env, con, pre_eig = con$pre_eig, trs = trs, interval = interval1)
+ interval1 <- spdep:::jacobianSetup(method, env, con, pre_eig = con$pre_eig, trs = trs1, interval = interval1)
assign("interval1", interval1, envir = env)
nm <- paste(method, "set_up", sep = "_")
timings[[nm]] <- proc.time() - .ptime_start
.ptime_start <- proc.time()
- RES<- sperrorlm(env = env, zero.policy = zero.policy, interval = interval1)
+ RES<- sperrorlm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess)
res.eff<-feerror(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method =method, lambda=RES$lambda, legacy = legacy)
}
-if(model=="sarar") {
-
- RES<- spsararlm(env = env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve = tol.solve)
-
-
-res.eff<-felag(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method = method, rho = RES$lambda, legacy = legacy, zero.policy = zero.policy)
- }
@@ -478,6 +475,7 @@
var[((2+1):ncol(var)),((2+1):ncol(var))]<-RES$asyvar1
}
+
spmod <- list(coefficients=Coeff, errcomp=NULL,
vcov = var ,spat.coef=spat.coef,
vcov.errcomp=NULL,
@@ -485,6 +483,10 @@
sigma2=RES$s2, type=type, model=model.data,
call=cl, logLik=RES$ll, method=method, effects=effects,
res.eff=res.eff)
+
+if (!is.null(na.act))
+ spmod$na.action <- na.act
+
class(spmod) <- "splm"
return(spmod)
Modified: pkg/R/spml.R
===================================================================
--- pkg/R/spml.R 2013-03-21 15:03:55 UTC (rev 149)
+++ pkg/R/spml.R 2013-03-21 21:17:13 UTC (rev 150)
@@ -1,4 +1,4 @@
-spml <- function(formula, data, index=NULL, listw, listw2=listw,
+spml <- function(formula, data, index=NULL, listw, listw2=listw, na.action,
model=c("within","random","pooling"),
effect=c("individual","time","twoways"),
lag=FALSE, spatial.error=c("b","kkp","none"),
@@ -35,7 +35,7 @@
effects <- switch(match.arg(effect), individual="spfe",
time="tpfe", twoways="sptpfe")
res <- spfeml(formula=formula, data=data, index=index,
- listw=listw, listw2=listw2,
+ listw=listw, listw2=listw2, na.action,
model=model, effects=effects,
cl=cl, ...)
}, random={
Modified: pkg/man/spml.Rd
===================================================================
--- pkg/man/spml.Rd 2013-03-21 15:03:55 UTC (rev 149)
+++ pkg/man/spml.Rd 2013-03-21 21:17:13 UTC (rev 150)
@@ -6,7 +6,7 @@
\usage{
-spml(formula, data, index = NULL, listw, listw2=listw,
+spml(formula, data, index=NULL, listw, listw2=listw, na.action,
model=c("within","random","pooling"),
effect=c("individual","time","twoways"),
lag=FALSE, spatial.error=c("b","kkp","none"),
@@ -22,6 +22,7 @@
\item{listw2}{an object of class \code{listw} or a
\code{matrix}. Second of set spatial weights for estimation, if
different from the first (e.g., in a 'sarar' model).}
+ \item{na.action}{see \pkg{spdep} for more details.}
\item{model}{one of \code{c("within", "random", "pooling").}}
\item{effect}{one of \code{c("individual","time","twoways")}; the
effects introduced in the model.}
@@ -87,7 +88,7 @@
## the two standard specifications (SEM and SAR) one with FE
## and the other with RE:
## fixed effects panel with spatial errors
-fespaterr <- spml(fm, data = Produc, listw = mat2listw(usaww), model="within", spatial.error="b")
+fespaterr <- spml(fm, data = Produc, listw = mat2listw(usaww), model="within", spatial.error="b", Hess = TRUE)
summary(fespaterr)
## random effects panel with spatial lag
respatlag <- spml(fm, data = Produc, listw = mat2listw(usaww), model="random", spatial.error="none", lag=TRUE)
More information about the Splm-commits
mailing list