[Splm-commits] r151 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 22 16:50:20 CET 2013
Author: gpiras
Date: 2013-03-22 16:50:20 +0100 (Fri, 22 Mar 2013)
New Revision: 151
Modified:
pkg/R/fixed_effects.R
pkg/R/likelihoodsFE.R
pkg/R/spfeml.R
Log:
updated spfeml
Modified: pkg/R/fixed_effects.R
===================================================================
--- pkg/R/fixed_effects.R 2013-03-21 21:17:13 UTC (rev 150)
+++ pkg/R/fixed_effects.R 2013-03-22 15:50:20 UTC (rev 151)
@@ -96,7 +96,7 @@
#felag<-function(y,x,wy,ysms,xsms,ytms, xtms, wytms, wysms, beta,sige,yt,xt,N,T,NT,k,effects,method, rho,listw,inde){
-felag<-function(env, beta,sige, effects, method, rho, legacy, zero.policy){
+felag<-function(env, beta,sige, effects, method, lambda, legacy, zero.policy){
y<-get("y", envir = env)
x<-get("x", envir = env)
wy<-get("wy", envir = env)
@@ -112,14 +112,14 @@
inde<-get("inde", envir = env)
mx<-apply(x,2,mean)
- intercept <- mean(y)- mean(wy)*rho - mx%*%beta
+ intercept <- mean(y)- mean(wy)*lambda - mx%*%beta
if (effects=="spfe"){
ysms<-get("ysms", envir = env)
xsms<-get("xsms", envir = env)
wysms<-get("wysms", envir = env)
- res.sfe <- as.matrix(ysms) - as.matrix(wysms) *rho - xsms %*% as.matrix(beta) - as.numeric(intercept)
+ res.sfe <- as.matrix(ysms) - as.matrix(wysms) *lambda - xsms %*% as.matrix(beta) - as.numeric(intercept)
xhat <- x %*% as.matrix(beta) + rep(res.sfe,T) + as.numeric(intercept)
res.var.sfe<- (sige / T) + diag((as.numeric(sige)*(xsms%*% solve(crossprod(xt)) %*% t(xsms))))
res.se.sfe<-sqrt(res.var.sfe)
@@ -127,7 +127,7 @@
res.se.con<-sqrt(as.numeric(sige) / NT + as.numeric(sige) * t(as.matrix(mx)) %*% solve(crossprod(xt)) %*% as.matrix(mx))
res.t.con <- as.numeric(intercept) / res.se.con
N.vars <- k + N
- res.e <- y - xhat - rho* wy
+ res.e <- y - xhat - lambda* wy
FE.out<-list(res.sfe=res.sfe, res.se.sfe=res.se.sfe, intercept=intercept, res.se.con=res.se.con,xhat=xhat,N.vars=N.vars,res.e=res.e)
}
@@ -136,7 +136,7 @@
xtms<-get("xtms", envir = env)
wytms<-get("wytms", envir = env)
- res.tfe <- as.matrix(ytms) - as.matrix(wytms)* rho - xtms %*% as.matrix(beta) - as.numeric(intercept)
+ res.tfe <- as.matrix(ytms) - as.matrix(wytms)* lambda - xtms %*% as.matrix(beta) - as.numeric(intercept)
xhat <- x %*% as.matrix(beta) + rep(res.tfe,each=N) + as.numeric(intercept)
res.var.tfe <- (sige / N) + (as.numeric(sige)*(xtms%*% solve(crossprod(xt)) %*% t(xtms)))
res.se.tfe <-sqrt(diag(res.var.tfe))
@@ -144,7 +144,7 @@
res.se.con<-sqrt(as.numeric(sige) / NT + as.numeric(sige) * t(as.matrix(mx)) %*% solve(crossprod(xt)) %*% as.matrix(mx))
res.t.con <- as.numeric(intercept) / res.se.con
N.vars <- k + T
- res.e <- y - xhat - rho* wy
+ res.e <- y - xhat - lambda * wy
FE.out<-list(res.tfe=res.tfe, res.se.tfe=res.se.tfe, intercept=intercept, res.se.con=res.se.con,xhat=xhat,N.vars=N.vars,res.e=res.e)
}
if (effects== "sptpfe"){
@@ -156,8 +156,8 @@
xtms<-get("xtms", envir = env)
wytms<-get("wytms", envir = env)
- res.sfe <- as.matrix(ysms) - as.matrix(wysms) * rho - xsms %*% as.matrix(beta) - as.numeric(intercept)
- res.tfe <- as.matrix(ytms) - as.matrix(wytms) * rho - xtms %*% as.matrix(beta) - as.numeric(intercept)
+ res.sfe <- as.matrix(ysms) - as.matrix(wysms) * lambda - xsms %*% as.matrix(beta) - as.numeric(intercept)
+ res.tfe <- as.matrix(ytms) - as.matrix(wytms) * lambda - xtms %*% as.matrix(beta) - as.numeric(intercept)
res.var.sfe<- (sige / T) + (as.numeric(sige)*(xsms%*% solve(crossprod(xt)) %*% t(xsms)))
res.se.sfe <-sqrt(diag(res.var.sfe))
res.var.tfe <- (sige / N) + (as.numeric(sige)*(xtms%*% solve(crossprod(xt)) %*% t(xtms)))
@@ -168,18 +168,18 @@
res.t.con <- as.numeric(intercept) / res.se.con
xhat<- x %*% as.matrix(beta) + rep(res.sfe,T) + rep(res.tfe,each=N) + as.numeric(intercept)
N.vars <- k + N + T - 1
- res.e <- y - xhat - rho* wy
+ res.e <- y - xhat - lambda * wy
FE.out<-list(res.tfe=res.tfe, res.se.tfe=res.se.tfe, res.sfe=res.sfe, res.se.sfe=res.se.sfe, intercept=intercept, res.se.con=res.se.con,xhat=xhat,N.vars=N.vars,res.e=res.e)
}
if (effects=="pooled") {
xhat <- x %*% as.matrix(beta)
- res.e <- y - xhat - rho* wy
+ res.e <- y - xhat - lambda* wy
FE.out<-list(xhat=xhat,N.vars=k,res.e=res.e)
}
if(legacy){
W <- listw2dgCMatrix(listw, zero.policy = zero.policy)
- IrW<-Diagonal(N)-rho*W
+ IrW<- sparseMatrix(i=1:N, j=1:N, x=1) -lambda*W
IrWi<-solve(IrW)
xtb <- xt %*% beta
yhat <- unlist(tapply(xhat,inde, function(u) IrWi %*% u))
Modified: pkg/R/likelihoodsFE.R
===================================================================
--- pkg/R/likelihoodsFE.R 2013-03-21 21:17:13 UTC (rev 150)
+++ pkg/R/likelihoodsFE.R 2013-03-22 15:50:20 UTC (rev 151)
@@ -34,7 +34,7 @@
inde<-get("inde", envir = env)
interval1 <- get("interval1", envir = env)
- XpX<-crossprod(xt)
+ XpX<-crossprod(xt)
b0<-solve(XpX,crossprod(xt,yt)) ####y on X
b1<-solve(XpX,crossprod(xt,wyt)) ####Wy on x
e0<-yt - xt%*% b0
@@ -53,48 +53,47 @@
#opt <- nlminb(0.02138744, conclikpan, lower = interval[1], upper= interval[2], env = env)
- rho <- opt$maximum
+ lambda <- opt$maximum
- if (isTRUE(all.equal(rho, interval[1])) || isTRUE(all.equal(rho,
- interval[2])))
- warning("rho on interval bound - results should not be used")
+ if (isTRUE(all.equal(lambda, interval[1])) || isTRUE(all.equal(lambda,interval[2])))
+ warning("lambda on interval bound - results should not be used")
- names(rho) <- "lambda"
+ names(lambda) <- "lambda"
LL <- opt$objective
optres <- opt
- nm <- paste(method, "opt", sep = "_")
- timings[[nm]] <- proc.time() - .ptime_start
- .ptime_start <- proc.time()
+ # nm <- paste(method, "opt", sep = "_")
+ # timings[[nm]] <- proc.time() - .ptime_start
+ # .ptime_start <- proc.time()
- lm.lag <- lm((yt - rho * wyt) ~ xt - 1)
+ lm.lag <- lm((yt - lambda * wyt) ~ xt - 1)
p <- lm.lag$rank
r <- residuals(lm.lag)
fit <- yt - r
names(r) <- names(fit)
betas <- coefficients(lm.lag)
names(betas) <- colnames(xt)
- coefs <- c(rho, betas)
+ coefs <- c(lambda, betas)
SSE <- deviance(lm.lag)
s2 <- SSE/NT
- coefsl <- c(s2, rho, betas)
+ # coefsl <- c(s2, rho, betas)
- timings[["coefs"]] <- proc.time() - .ptime_start
- .ptime_start <- proc.time()
- assign("first_time", TRUE, envir = env)
+ # timings[["coefs"]] <- proc.time() - .ptime_start
+ # .ptime_start <- proc.time()
+ # assign("first_time", TRUE, envir = env)
### add numerical hessian
-if(!Hess){
+if(Hess){
- fd <- fdHess(coefsl, f_sarpanel_hess, env)
+ fd <- fdHess(coefs, f_sarpanel_hess, env)
mat <- fd$Hessian
fdHess<- solve(-(mat), tol.solve = tol.solve)
- rownames(fdHess) <- colnames(fdHess) <- c("s2", "lambda", colnames(xt))
+ rownames(fdHess) <- colnames(fdHess) <- c("lambda", colnames(xt))
- rho.se <- fdHess[2, 2]
- sig.se <- fdHess[1, 1]
+ lambda.se <- fdHess[1, 1]
+ sig.se <- NULL
rest.se <- vcov(lm.lag)
@@ -104,7 +103,7 @@
tr <- function(A) sum(diag(A))
W <-listw2dgCMatrix(listw, zero.policy = zero.policy)
- A <- solve(diag(NT/T) - rho * W)
+ A <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - lambda * W)
WA <- W %*% A
one <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
@@ -128,7 +127,7 @@
asyv <- solve(asyvar, tol = con$tol.solve)
rownames(asyv) <- colnames(asyv) <- c("sigma","lambda", colnames(xt))
- rho.se <- sqrt(asyv[2, 2])
+ lambda.se <- sqrt(asyv[2, 2])
rest.se <- sqrt(diag(asyv))[-c(1:2)]
sig.se <- sqrt(asyv[1, 1])
asyvar1 <- asyv[-1,-1]
@@ -137,7 +136,7 @@
}
- return<-list(coeff=betas,rho=rho,s2=s2, rest.se=rest.se, rho.se=rho.se, sig.se = sig.se, asyvar1=asyvar1)
+ return<-list(coeff = betas, lambda = lambda, s2 = s2, rest.se = rest.se, lambda.se = lambda.se, sig.se = sig.se, asyvar1 = asyvar1)
}
@@ -147,22 +146,33 @@
T<-get("T", envir = env)
NT<-get("NT", envir = env)
- s2<- coefs[1]
- lambda <- coefs[2]
- beta <- coefs[-c(1,2)]
- SSE <- s2*n
+
+ lambda <- coefs[1]
+ beta <- coefs[-1]
+
n <- NT/T
+ SSE <- sar_hess_sse_panel(lambda, beta, env)
+ s2 <- SSE /n
+
ldet <- do_ldet(lambda, env, which = 1)
ret <- (T * ldet - ((n*T/2) * log(2 * pi)) - (n*T/2) * log(s2) -
(1/(2 * s2)) * SSE)
if (get("verbose", envir = env))
- cat("lambda:", lambda, " function:", ret,
- " Jacobian:", ldet," SSE:", SSE, "\n")
+ cat("lambda:", lambda, " function:", ret, " Jacobian:", ldet," SSE:", SSE, "\n")
ret
}
+sar_hess_sse_panel <- function (lambda, beta, env)
+{
+ yl <- get("yt", envir = env) - lambda * get("wyt", envir = env)
+ res <- yl - (get("xt", envir = env) %*% beta)
+ SSE <- c(crossprod(res))
+ SSE
+}
+
+
####### ERROR MODEL
-sarpanelerror<-function (lambda, env=env)
+sarpanelerror<-function (rho, env=env)
{
yt<- get("yt", envir = env)
xt<- get("xt", envir = env)
@@ -176,18 +186,18 @@
inde<- get("inde", envir = env)
T<- get("T", envir = env)
- yco <- yt - lambda * wyt
- xco <- xt - lambda * wxt
+ yco <- yt - rho * wyt
+ xco <- xt - rho * wxt
bb<- solve(crossprod(xco),crossprod(xco, yco) )
ehat<- yco - xco %*% bb
SSE <- crossprod(ehat)
- ldet <- do_ldet(lambda, env)
+ ldet <- do_ldet(rho, env)
ret <- T*ldet - (NT/2) * log(SSE)
if (get("verbose", envir = env))
- cat("rho:", lambda, " function:", ret, " Jacobian:", ldet, " SSE:", SSE, "\n")
+ cat("rho:", rho, " function:", ret, " Jacobian:", ldet, " SSE:", SSE, "\n")
ret
}
@@ -472,7 +482,8 @@
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))
+ rownames(fdHess) <- colnames(fdHess) <- c("lambda", "rho",colnames(xt))
+
rho.se <- fdHess[2, 2]
lambda.se <- fdHess[1, 1]
rest.se <- vcov(lm.target)
Modified: pkg/R/spfeml.R
===================================================================
--- pkg/R/spfeml.R 2013-03-21 21:17:13 UTC (rev 150)
+++ pkg/R/spfeml.R 2013-03-22 15:50:20 UTC (rev 151)
@@ -1,8 +1,8 @@
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()
- .ptime_start <- proc.time()
+ # timings <- list()
+ # .ptime_start <- proc.time()
model<-match.arg(model)
@@ -369,7 +369,7 @@
}
assign("verbose", !quiet, envir = env)
-assign("first_time", TRUE, envir = env)
+# assign("first_time", TRUE, envir = env)
assign("LAPACK", con$LAPACK, envir = env)
assign("can.sim", can.sim, envir=env)
assign("similar", FALSE, envir = env)
@@ -378,8 +378,8 @@
assign("con", con, envir=env)
-timings[["set_up"]] <- proc.time() - .ptime_start
-.ptime_start <- proc.time()
+# timings[["set_up"]] <- proc.time() - .ptime_start
+# .ptime_start <- proc.time()
if (!quiet)
cat(paste("\nSpatial autoregressive error model\n", "Jacobian calculated using "))
@@ -387,13 +387,14 @@
if(model == "lag"){
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()
+ # 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, 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)
+
+ res.eff<-felag(env = env, beta = RES$coeff, sige = RES$s2, effects = effects, method = method, lambda = RES$lambda, legacy = legacy, zero.policy = zero.policy)
}
@@ -403,14 +404,14 @@
assign("interval1", interval1, envir = env)
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()
+ # 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, 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)
+res.eff<-felag(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method = method, lambda = RES$lambda, legacy = legacy, zero.policy = zero.policy)
}
@@ -420,9 +421,9 @@
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()
+ # 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, Hess = Hess)
res.eff<-feerror(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method =method, lambda=RES$lambda, legacy = legacy)
More information about the Splm-commits
mailing list