[Splm-commits] r183 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 11 02:43:42 CET 2013
Author: gpiras
Date: 2013-12-11 02:43:41 +0100 (Wed, 11 Dec 2013)
New Revision: 183
Modified:
pkg/R/likelihoodsFE.R
pkg/R/spfeml.R
Log:
some changes to spfeml
Modified: pkg/R/likelihoodsFE.R
===================================================================
--- pkg/R/likelihoodsFE.R 2013-12-09 22:28:29 UTC (rev 182)
+++ pkg/R/likelihoodsFE.R 2013-12-11 01:43:41 UTC (rev 183)
@@ -44,6 +44,7 @@
e1e1<-crossprod(e1)
e0e1<-t(e1)%*%e0
+
assign("e0e0", e0e0, envir = env)
assign("e1e1", e1e1, envir = env)
assign("e0e1", e0e1, envir = env)
@@ -123,7 +124,7 @@
}
-
+
### add numerical hessian
if(Hess){
@@ -183,13 +184,19 @@
lambda.se <- asyv[2, 2]
rest.se <- sqrt(diag(asyv))[-c(1:2)]
sig.se <- sqrt(asyv[1, 1])
- asyvar1 <- asyv[-c(1,2),-c(1,2)]
-
- rownames(asyvar1) <- colnames(asyvar1) <- c(colnames(xt))
+ asyvar1 <- as.matrix(asyv[-c(1,2),-c(1,2)])
+ rownames(asyvar1) <- colnames(asyvar1) <- colnames(xt)
+
+
}
+
+if(Hess) asyv <- NULL
+else asyv <- asyv
+
+
- return<-list(coeff = betas, lambda = lambda, s2 = s2, rest.se = rest.se, lambda.se = lambda.se, sig.se = sig.se, asyvar1 = asyvar1, residuals = r)
+ return<-list(coeff = betas, lambda = lambda, s2 = s2, rest.se = rest.se, lambda.se = lambda.se, sig.se = sig.se, asyvar1 = asyvar1, residuals = r, asyv = asyv)
}
@@ -411,9 +418,15 @@
asyvar1 <- asyv[-c(1,2),-c(1,2)]
# rownames(asyvar1) <- colnames(asyvar1) <- c(colnames(xt))
+
+
}
- return<-list(coeff=betas, rho = rho, s2 = s2, rest.se = rest.se, rho.se = rho.se, s2.se = s2.se, asyvar1=asyvar1, residuals = r)
+if(Hess) asyv <- NULL
+else asyv <- asyv
+
+
+ return<-list(coeff=betas, rho = rho, s2 = s2, rest.se = rest.se, rho.se = rho.se, s2.se = s2.se, asyvar1=asyvar1, residuals = r, asyv = asyv)
}
@@ -802,10 +815,14 @@
rest.se <- sqrt(diag(asyv))[-((p+1):(p+3))]
asyvar1 <- asyv[-((p+1):(p+3)),-((p+1):(p+3))]
+
}
+if(Hess) asyv <- NULL
+else asyv <- asyv
+
-return<-list(coeff = betas, lambda = lambda, rho = rho, s2 = s2, asyvar1 = asyvar1, lambda.se = lambda.se, rho.se = rho.se, s2.se = s2.se, residuals = r)
+return<-list(coeff = betas, lambda = lambda, rho = rho, s2 = s2, asyvar1 = asyvar1, lambda.se = lambda.se, rho.se = rho.se, s2.se = s2.se, residuals = r, asyv = asyv)
}
f_sacpanel_hess <- function (coefs, env, LeeYu = LeeYu, effects = effects)
Modified: pkg/R/spfeml.R
===================================================================
--- pkg/R/spfeml.R 2013-12-09 22:28:29 UTC (rev 182)
+++ pkg/R/spfeml.R 2013-12-11 01:43:41 UTC (rev 183)
@@ -55,6 +55,9 @@
## reduce X,y to model matrix values (no NAs)
x<-model.matrix(formula,data=data)
+ clnames <- colnames(x)
+ rwnames <- rownames(x)
+
y<-model.response(model.frame(formula,data=data))
## reduce index accordingly
names(index)<-row.names(data)
@@ -68,32 +71,25 @@
ind<-ind[oo]
tind<-tind[oo]
- clnames<-colnames(x)
- rwnames<-rownames(x)
+
#make sure that the model has no intercept if effects !=pooled
- if (colnames(x)[1]=="(Intercept)") {
- x<-x[,-1]
- #cat('\n Warning: x may not contain an intercept if fixed effects are specified \n')
- }
-
-if(is.vector(x)){
- k<-1
- x<-matrix(x)
- colnames(x)<-clnames[-1]
- dimnames(x)[[1]]<-rwnames
- dimnames(x)[[2]]<-clnames[-1]
+ if (attr(attributes(model.frame(formula,data=data))$terms, "intercept") == 1) {
+ x <- as.matrix(x[,-1])
+ colnames(x)<-clnames[-1]
+ dimnames(x)[[1]]<-rwnames
+ clnames <- clnames[-1]
}
- else k<-dim(x)[[2]]
+ x <- as.matrix(x)
+ k <- dim(x)[2]
+
+
-
-
-
## det. number of groups and df
N<-length(unique(ind))
n<-N
- #x<-matrix(x,length(ind),k)
## det. max. group numerosity
- T<-max(tapply(x[,1],ind,length))
+ T<-max(tapply(tind,ind,length))
+
## det. total number of obs. (robust vs. unbalanced panels)
NT<-length(ind)
@@ -263,22 +259,22 @@
if (model == "error"){
dm<-function(A) trash<-unlist(tapply(A,inde,function(TT) lag.listw(listw,TT), simplify=TRUE))
wxt<-apply(xt,2,dm)
- colnames(wxt)<-paste('Lag.',colnames(x), sep="")
+ # colnames(wxt)<-paste('Lag.',colnames(x), sep="")
wx<-apply(x,2,dm)
- colnames(wx)<-paste('lag.',colnames(x), sep="")
+ # colnames(wx)<-paste('lag.',colnames(x), sep="")
}
if (model == "sarar"){
dm<-function(A) trash<-unlist(tapply(A,inde,function(TT) lag.listw(listw2,TT), simplify=TRUE))
wxt<-apply(xt,2,dm)
- colnames(wxt)<-paste('Lag.',colnames(x), sep="")
+ # colnames(wxt)<-paste('Lag.',colnames(x), sep="")
wx<-apply(x,2,dm)
- colnames(wx)<-paste('lag.',colnames(x), sep="")
+ # colnames(wx)<-paste('lag.',colnames(x), sep="")
}
+# print(clnames)
+colnames(xt)<- clnames
-colnames(xt)<-dimnames(x)[[2]]
-
assign("yt",yt, envir=env)
@@ -307,6 +303,8 @@
}
+
+
assign("verbose", !quiet, envir = env)
# assign("first_time", TRUE, envir = env)
assign("LAPACK", con$LAPACK, envir = env)
@@ -317,6 +315,7 @@
assign("con", con, envir=env)
+
if (!quiet)
cat(paste("\nSpatial fixed effects model\n", "Jacobian calculated using "))
@@ -358,7 +357,7 @@
# timings[[nm]] <- proc.time() - .ptime_start
# .ptime_start <- proc.time()
- RES<- sperrorlm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess, LeeYu = LeeYu, effects = effects)
+ RES <- sperrorlm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess, LeeYu = LeeYu, effects = effects)
res.eff<-feerror(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method =method, rho=RES$rho, legacy = legacy)
}
@@ -403,15 +402,17 @@
type <- paste("fixed effects", model)
-var<-RES$asyvar1
+if (!Hess) var <- RES$asyv
-if(model == "lag"){
+else{
+
+if(model == "lag" ){
var<-matrix(0,(ncol(RES$asyvar1)+1),(ncol(RES$asyvar1)+1))
var[1,1]<- RES$lambda.se
var[(2:ncol(var)),(2:ncol(var))]<-RES$asyvar1
}
-if(model == "error"){
+if(model == "error" ){
var<-matrix(0,(ncol(RES$asyvar1)+1),(ncol(RES$asyvar1)+1))
var[1,1]<- RES$rho.se
var[(2:ncol(var)),(2:ncol(var))]<-RES$asyvar1
@@ -424,6 +425,7 @@
var[(3:ncol(var)),(3:ncol(var))]<-RES$asyvar1
}
+}
spmod <- list(coefficients=Coeff, errcomp=NULL,
vcov = var ,spat.coef=spat.coef,
More information about the Splm-commits
mailing list