[Splm-commits] r78 - in pkg: . R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Jul 13 18:52:02 CEST 2010


Author: gpiras
Date: 2010-07-13 18:52:02 +0200 (Tue, 13 Jul 2010)
New Revision: 78

Added:
   pkg/R/
   pkg/R/Ggsararsp.R
   pkg/R/LMHtest.R
   pkg/R/LMHtest.model.R
   pkg/R/REmod.R
   pkg/R/UTILITIES_SPSEML.R
   pkg/R/arg.R
   pkg/R/arg1.R
   pkg/R/arg2.R
   pkg/R/bsjktest.R
   pkg/R/bsjktest.formula.R
   pkg/R/bsktest.R
   pkg/R/bsktest.formula.R
   pkg/R/bsktest.lm.R
   pkg/R/bsktest.splm.R
   pkg/R/clmltest.R
   pkg/R/clmltest.model.R
   pkg/R/clmmtest.R
   pkg/R/clmmtest.model.R
   pkg/R/conclikpan.M.R
   pkg/R/conclikpan.R
   pkg/R/conclikpan.sp.R
   pkg/R/dosparsepanel.R
   pkg/R/effects.R
   pkg/R/fixed_effects.R
   pkg/R/fs.R
   pkg/R/likelihoodserror.R
   pkg/R/likelihoodslag.R
   pkg/R/listw2U_Matrix.R
   pkg/R/listw2U_spam.R
   pkg/R/listw2dgCMatrix.R
   pkg/R/modtslssp.R
   pkg/R/pbsjkARtest.R
   pkg/R/pbsjkJtest.R
   pkg/R/pbsjkREtest.R
   pkg/R/pbsjkSDtest.R
   pkg/R/print.effects.splm.R
   pkg/R/print.splm.R
   pkg/R/print.summary.splm.R
   pkg/R/pw.R
   pkg/R/sarREmod.R
   pkg/R/sarpanelerror.M.R
   pkg/R/sarpanelerror.R
   pkg/R/sarpanelerror.sp.R
   pkg/R/searg.R
   pkg/R/semREmod.R
   pkg/R/semarREmod.R
   pkg/R/semarmod.R
   pkg/R/semmod.R
   pkg/R/similar.listw_Matrix.R
   pkg/R/similar.listw_spam.R
   pkg/R/slm1test.R
   pkg/R/slm1test.model.R
   pkg/R/slm2test.R
   pkg/R/slm2test.model.R
   pkg/R/sperrorlm.R
   pkg/R/spfegm.R
   pkg/R/spfeml.R
   pkg/R/splaglm.R
   pkg/R/spregm.R
   pkg/R/spreml.R
   pkg/R/spsegm.R
   pkg/R/spseml.R
   pkg/R/ssrREmod.R
   pkg/R/ssrmod.R
   pkg/R/summary.effects.splm.R.old
   pkg/R/summary.splm.R
   pkg/R/sumres.R
   pkg/R/tslssp.R
   pkg/R/tss.R
   pkg/R/tw.R
   pkg/R/utilities.R
   pkg/R/vcov.splm.R
   pkg/R/write.effects.splm.R
Log:
import new functions

Added: pkg/R/Ggsararsp.R
===================================================================
--- pkg/R/Ggsararsp.R	                        (rev 0)
+++ pkg/R/Ggsararsp.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,22 @@
+`Ggsararsp` <-
+function (W, u, zero.policy = FALSE) 
+{
+      n <- length(u)
+      tt<-matrix(0,n,1)
+      tr<-sum(unlist(W$weights)^2)
+      wu<-lag.listw(W,u)
+      wwu<-lag.listw(W,wu)
+    	uu <- crossprod(u, u)
+    	uwu <- crossprod(u, wu)
+ 	uwpuw <- crossprod(wu, wu)
+    	uwwu <- crossprod(u, wwu)
+    	wwupwu <- crossprod(wwu, wu)
+    	wwupwwu <- crossprod(wwu, wwu)
+    	bigG <- matrix(0, 3, 3)
+    	bigG[, 1] <- c(2 * uwu, 2 * wwupwu, (uwwu + uwpuw))/n
+    	bigG[, 2] <-  -c(uwpuw, wwupwwu, wwupwu)/n
+    	bigG[, 3] <- c(1, tr/n, 0)
+    	litg <- c(uu, uwpuw, uwu)/n
+    	list(bigG = bigG, litg = litg)
+}
+


Property changes on: pkg/R/Ggsararsp.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/LMHtest.R
===================================================================
--- pkg/R/LMHtest.R	                        (rev 0)
+++ pkg/R/LMHtest.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,97 @@
+`LMHtest` <-
+function(formula, data, index=NULL, listw){
+    ## depends on listw2dgCMatrix.R
+  if(!is.null(index)) { ####can be deleted when using the wrapper
+    require(plm)
+    data <- plm.data(data, index)
+    }
+
+  index <- data[,1]
+  tindex <- data[,2]
+
+  x<-model.matrix(formula,data=data)
+  y<-model.response(model.frame(formula,data=data))
+   cl<-match.call()
+	names(index)<-row.names(data)
+  ind<-index[which(names(index)%in%row.names(x))]
+  tind<-tindex[which(names(index)%in%row.names(x))]
+  ## reorder data by cross-sections, then time
+  oo<-order(tind,ind)
+  x<-x[oo,]
+  y<-y[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  ## det. number of groups and df
+  N<-length(unique(ind))
+  k<-dim(x)[[2]]
+  ## det. max. group numerosity
+  T<-max(tapply(x[,1],ind,length))
+  ## det. total number of obs. (robust vs. unbalanced panels)
+  NT<-length(ind)
+	ols<-lm(y~x)
+	XpXi<-solve(crossprod(x))
+   n<-dim(ols$model)[1]
+
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) ####indicator to get the cross-sectional observations
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) ####indicator to get the time periods observations
+	bOLS<-coefficients(ols)
+	e<-as.matrix(residuals(ols))
+	ee<-crossprod(e)
+####calculate the elements of LMj, LM1, SLM1
+
+		JIe<-tapply(e,inde1,sum)
+		JIe<-rep(JIe,T) ####calculates (J_T kronecker I_N)*u
+		G<-(crossprod(e,JIe)/crossprod(e))-1 ###calculate G in LMj (same notation as in the paper)
+tr<-function(R) sum(diag(R))
+		LM1<-sqrt((NT/(2*(T-1))))*as.numeric(G) ###same notation as in Baltagi et al.
+
+
+####calculate the elements of LMj, LM1, SLM1
+		Wst<-listw2dgCMatrix(listw) ###transform the listw object in a sparse matrix
+		Ws<-t(Wst)  ### this is the real W since listw2dgCMatrix generate W'
+		WWp<-(Ws+Wst)/2 ##generate (W+W')/2
+yy<-function(q){ #### for very big dimension of the data this can be changed looping over the rows and columns of W or either the listw object
+	wq<-WWp%*%q
+	wq<-as.matrix(wq)
+	}
+		IWWpe<-unlist(tapply(e,inde,yy)) ####calculates (I_T kronecker (W+W')/2)*u
+		H<-crossprod(e,IWWpe)/crossprod(e) #calculate H (same notation as in the paper)
+		W2<-Ws%*%Ws ####generate W^2
+		WW<-crossprod(Ws) ####generate W'*W
+		b<-tr(W2+WW) ###generates b (same notation as the paper)
+#		LMj<-(NT/(2*(T-1)))*as.numeric(G)^2 + ((N^2*T)/b)*as.numeric(H)^2   ###LMj as in the paper
+		LM2<-sqrt((N^2*T)/b)*as.numeric(H)^2 ###same notation as in Baltagi et al.
+if (LM1<=0){
+		if (LM2<=0) JOINT<-0
+		else JOINT<-LM2^2
+		}		####this is chi-square_m in teh notation of the paper.
+	else{
+		if (LM2<=0) JOINT<-LM1^2
+		else JOINT<-LM1^2 + LM2^2
+		}
+STAT<- qchisq(0.05,1,lower.tail=FALSE)
+STAT1<- qchisq(0.05,2,lower.tail=FALSE)
+if (JOINT>=2.952) {
+		if (JOINT<7.289 & JOINT>=4.321) pval<-0.05
+		if (JOINT >= 7.289) pval<-0.01
+		if (JOINT<= 4.321)	pval<-0.1
+	}
+else pval<-1
+
+	statistics<-JOINT
+
+  names(statistics)="LM-H"
+	method<- "Baltagi, Song and Koh LM-H one-sided joint test"
+	#alt<-"serial corr. in error terms, sub RE and spatial dependence"
+  ##(insert usual htest features)
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(formula), alternative="Random Regional Effects and Spatial autocorrelation")
+  class(RVAL) <- "htest"
+  return(RVAL)
+}
+


Property changes on: pkg/R/LMHtest.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/LMHtest.model.R
===================================================================
--- pkg/R/LMHtest.model.R	                        (rev 0)
+++ pkg/R/LMHtest.model.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,102 @@
+`LMHtest.model` <-
+function(x, listw, index){
+## depends on listw2dgCMatrix.R
+
+if(!inherits(x,"lm")) stop("argument should be an object of class lm")
+
+  if(is.null(index))  stop("index should be specified to retrieve information on time and cross-sectional dimentions")
+
+  if(!inherits(listw,"listw")) stop("object w should be of class listw")
+
+  ind <- index[,1]
+  tind <- index[,2]
+
+###extract objects from x
+  y<-model.response(x$model)
+  e<-as.matrix(residuals(x))
+  	ee<-crossprod(e)
+   n<-dim(x$model)[1]
+  	bOLS<-coefficients(x)
+  	  	  form<-x$call
+  x<-model.matrix(eval(x$call),x$model)
+  #print(x)
+	XpXi<-solve(crossprod(x))
+
+   cl<-match.call()
+  ## reorder data by cross-sections, then time
+  oo<-order(tind,ind)
+  x<-x[oo,]
+  y<-y[oo]
+  e<-e[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  ## det. number of groups and df
+  N<-length(unique(ind))
+  k<-dim(x)[[2]]
+  ## det. max. group numerosity
+  T<-max(tapply(x[,1],ind,length))
+  ## det. total number of obs. (robust vs. unbalanced panels)
+  NT<-length(ind)
+#print(c(N,k,T,NT))
+#   print(ols$model)
+#	k<-dim(ols$model)[2]-1
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) ####indicator to get the cross-sectional observations
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) ####indicator to get the time periods observations
+####calculate the elements of LMj, LM1, SLM1
+
+		JIe<-tapply(e,inde1,sum)
+		JIe<-rep(JIe,T) ####calculates (J_T kronecker I_N)*u
+		G<-(crossprod(e,JIe)/crossprod(e))-1 ###calculate G in LMj (same notation as in the paper)
+tr<-function(R) sum(diag(R))
+		LM1<-sqrt((NT/(2*(T-1))))*as.numeric(G) ###same notation as in Baltagi et al.
+
+
+####calculate the elements of LMj, LM1, SLM1
+		Wst<-listw2dgCMatrix(listw) ###transform the listw object in a sparse matrix
+		Ws<-t(Wst)  ### this is the real W since listw2dgCMatrix generate W'
+		WWp<-(Ws+Wst)/2 ##generate (W+W')/2
+yy<-function(q){ #### for very big dimension of the data this can be changed looping over the rows and columns of W or either the listw object
+	wq<-WWp%*%q
+	wq<-as.matrix(wq)
+	}
+		IWWpe<-unlist(tapply(e,inde,yy)) ####calculates (I_T kronecker (W+W')/2)*u
+		H<-crossprod(e,IWWpe)/crossprod(e) #calculate H (same notation as in the paper)
+		W2<-Ws%*%Ws ####generate W^2
+		WW<-crossprod(Ws) ####generate W'*W
+		b<-tr(W2+WW) ###generates b (same notation as the paper)
+#		LMj<-(NT/(2*(T-1)))*as.numeric(G)^2 + ((N^2*T)/b)*as.numeric(H)^2   ###LMj as in the paper
+		LM2<-sqrt((N^2*T)/b)*as.numeric(H)^2 ###same notation as in Baltagi et al.
+if (LM1<=0){
+		if (LM2<=0) JOINT<-0
+		else JOINT<-LM2^2
+		}		####this is chi-square_m in teh notation of the paper.
+	else{
+		if (LM2<=0) JOINT<-LM1^2
+		else JOINT<-LM1^2 + LM2^2
+		}
+STAT<- qchisq(0.05,1,lower.tail=FALSE)
+STAT1<- qchisq(0.05,2,lower.tail=FALSE)
+if (JOINT>=2.952) {
+		if (JOINT<7.289 & JOINT>=4.321) pval<-0.05
+		if (JOINT >= 7.289) pval<-0.01
+		if (JOINT<= 4.321)	pval<-0.1
+	}
+else pval<-1
+
+	statistics<-JOINT
+
+  names(statistics)="LM-H"
+	method<- "Baltagi, Song and Koh LM-H one-sided joint test"
+	#alt<-"serial corr. in error terms, sub RE and spatial dependence"
+  ##(insert usual htest features)
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(form), alternative="Random Regional Effects and Spatial autocorrelation")
+  class(RVAL) <- "htest"
+  return(RVAL)
+}
+


Property changes on: pkg/R/LMHtest.model.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/REmod.R
===================================================================
--- pkg/R/REmod.R	                        (rev 0)
+++ pkg/R/REmod.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,101 @@
+`REmod` <-
+function(X, y, ind, tind, n, k, t, nT, w=NULL, coef0=0,
+         hess=FALSE, trace=trace, x.tol=1.5e-18, rel.tol=1e-15, ...) {
+  ## random effects panel model 'vanilla' ML estimation
+  ## based on general framework, RE structure on errors
+
+  ## MUCH room for optimization! Take framework from ssrRE,
+  ## using kinship etc. (notice inversion of order(ind,tind)!)
+
+  ## Giovanni Millo, Trieste; this version: 22/10/2008.
+
+  ## some useful pieces:
+    Jt<-matrix(1,ncol=t,nrow=t)
+    In<-diag(1,n)
+    It<-diag(1,t)
+    Jbart<-Jt/t
+    Et<-It-Jbart
+
+  ## inverse of Sigma
+    invSigma <- function(phi, n, t) {
+			 invSigma <- 1/(t*phi+1) * kronecker(Jbart, In) + kronecker(Et, In)
+                   invSigma
+                   }
+
+  ## outsource determinant of sigma with an algebraically efficient form!
+
+  ## concentrated likelihood
+    ll.c<-function(phi, y, X, n, t) {
+
+        ## perform GLS
+
+            sigma.1<-invSigma(phi,n,t)
+
+            b.hat<-solve( crossprod(X,sigma.1)%*%X, crossprod(X,sigma.1)%*%y )
+            ehat<-y-X%*%b.hat
+            sigma2ehat<-crossprod(ehat,sigma.1)%*%ehat/(n*t)
+            bhat<-list(betahat=b.hat,e=ehat,sigma2=sigma2ehat)
+
+            e <- bhat[[2]]
+            s2e <- bhat[[3]]
+
+            due <- -n*t/2*log(s2e)
+            tre <- -n/2*log(t*phi+1)
+            quattro <- -1/(2*s2e)*crossprod(e,sigma.1)%*%e
+
+            const <- -(n*t)/2*log(2*pi)
+            ll.c <- const+due+tre+quattro
+            llc <- - ll.c
+            }
+
+  ## iterate (=traballa) until convergence:
+
+  myphi0 <- coef0
+
+  optimum<-nlminb(myphi0, ll.c,
+                  lower=1e-08, upper=1e08,
+                  control=list(x.tol=x.tol, rel.tol=rel.tol, trace=trace),
+                  y=y, X=X, n=n, t=t, ...)
+
+
+  myphi<-optimum$par
+  myll <- optimum$objective
+
+
+  ## optimal values of parms:
+  phi<-myphi
+
+  ## perform GLS
+
+            sigma.1<-invSigma(phi,n,t)
+
+            b.hat<-solve( crossprod(X,sigma.1)%*%X, crossprod(X,sigma.1)%*%y )
+            ehat<-y-X%*%b.hat
+            sigma2ehat<-crossprod(ehat,sigma.1)%*%ehat/(n*t)
+            beta<-list(betahat=b.hat,e=ehat,sigma2=sigma2ehat)
+
+  ## names for coefs and error comp.s
+  nam.beta <- dimnames(X)[[2]]
+  nam.errcomp <- "phi"
+
+  ## calc. cov(b) by GLS
+  covB<-as.numeric(beta[[3]])*solve(crossprod(X,sigma.1)%*%X)
+  dimnames(covB) <- list(nam.beta, nam.beta)
+
+  ## calc. cov(phi) by numerical Hessian
+  covPRL <- solve(-fdHess(myphi, function(x) -ll.c(x,y,X,n,t))$Hessian)
+  dimnames(covPRL) <- list(nam.errcomp, nam.errcomp)
+
+
+  ## make (separate) coefficients' vectors
+  betas <- as.vector(beta[[1]])
+  errcomp <- phi
+  names(betas) <- nam.beta
+  names(errcomp) <- nam.errcomp
+
+  RES <- list(betas=betas, errcomp=errcomp,
+              covB=covB, covPRL=covPRL, ll=myll)
+
+  return(RES)
+  }
+


Property changes on: pkg/R/REmod.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/UTILITIES_SPSEML.R
===================================================================
--- pkg/R/UTILITIES_SPSEML.R	                        (rev 0)
+++ pkg/R/UTILITIES_SPSEML.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,151 @@
+	can.be.simmed <- function(listw) {
+	res <- is.symmetric.nb(listw$neighbours, FALSE)
+	if (res) {
+		if (attr(listw$weights, "mode") == "general")
+			res <- attr(listw$weights, "glistsym")
+	} else return(res)
+	res
+}
+
+
+similar.listw_Matrix <- function(listw) {
+	nbsym <- attr(listw$neighbours, "sym")
+	if(is.null(nbsym)) nbsym <- is.symmetric.nb(listw$neighbours, FALSE)
+	if (!nbsym) 
+		stop("Only symmetric nb can yield similar to symmetric weights")
+	if (attr(listw$weights, "mode") == "general")
+		if (!attr(listw$weights, "glistsym"))
+			stop("General weights must be symmetric")
+	n <- length(listw$neighbours)
+	if (n < 1) stop("non-positive number of entities")
+	ww <- as(as_dgRMatrix_listw(listw), "CsparseMatrix")
+	if (listw$style == "W") {
+		d <- attr(listw$weights, "comp")$d
+		d1 <- 1/(sqrt(d))
+		dd <- as(as(Diagonal(x=d), "symmetricMatrix"), "CsparseMatrix")
+		dd1 <- as(as(Diagonal(x=d1), "symmetricMatrix"),
+		    "CsparseMatrix")
+		ww1 <- dd %*% ww
+		res <- dd1 %*% ww1 %*% dd1
+	} else if (listw$style == "S") {
+		q <- attr(listw$weights, "comp")$q
+		Q <- attr(listw$weights, "comp")$Q
+		eff.n <- attr(listw$weights, "comp")$eff.n
+		q1 <- 1/(sqrt(q))
+		qq <- as(as(Diagonal(x=q), "symmetricMatrix"), "CsparseMatrix")
+		qq1 <- as(as(Diagonal(x=q1), "symmetricMatrix"),
+		    "CsparseMatrix")
+		ww0 <- (Q/eff.n) * ww
+		ww1 <- qq %*% ww0
+		sim0 <- qq1 %*% ww1 %*% qq1
+		res <- (eff.n/Q) * sim0
+	} else stop("Conversion not suitable for this weights style")
+	res
+}
+
+
+similar.listw_spam <- function(listw) {
+        if (!require(spam)) stop("spam not available")
+	nbsym <- attr(listw$neighbours, "sym")
+	if(is.null(nbsym)) nbsym <- is.symmetric.nb(listw$neighbours, FALSE)
+	if (!nbsym) 
+		stop("Only symmetric nb can yield similar to symmetric weights")
+	if (attr(listw$weights, "mode") == "general")
+		if (!attr(listw$weights, "glistsym"))
+			stop("General weights must be symmetric")
+	n <- length(listw$neighbours)
+	if (n < 1) stop("non-positive number of entities")
+	sww <- as.spam.listw(listw)
+	if (listw$style == "W") {
+		sd <- attr(listw$weights, "comp")$d
+		sd1 <- 1/(sqrt(sd))
+		sdd <- diag.spam(sd, n, n)
+		sdd1 <- diag.spam(sd1, n, n)
+		sww1 <- sdd %*% sww
+		res <- sdd1 %*% sww1 %*% sdd1
+	} else if (listw$style == "S") {
+		q <- attr(listw$weights, "comp")$q
+		Q <- attr(listw$weights, "comp")$Q
+		eff.n <- attr(listw$weights, "comp")$eff.n
+		q1 <- 1/(sqrt(q))
+		qq <- diag.spam(q, n, n)
+		qq1 <- diag.spam(q1, n, n)
+		ww0 <- (Q/eff.n) * sww
+		ww1 <- qq %*% ww0
+		sim0 <- qq1 %*% ww1 %*% qq1
+		res <- (eff.n/Q) * sim0
+	} else stop("Conversion not suitable for this weights style")
+	res
+}
+
+
+
+similar.listw <- function(listw) {
+	nbsym <- attr(listw$neighbours, "sym")
+	if(is.null(nbsym)) nbsym <- is.symmetric.nb(listw$neighbours, FALSE)
+	if (!nbsym) 
+		stop("Only symmetric nb can yield similar to symmetric weights")
+	if (attr(listw$weights, "mode") == "general")
+		if (!attr(listw$weights, "glistsym"))
+			stop("General weights must be symmetric")
+	n <- length(listw$neighbours)
+	if (n < 1) stop("non-positive number of entities")
+	cardnb <- card(listw$neighbours)
+	if (listw$style == "W") {
+		d <- attr(listw$weights, "comp")$d
+		glist <- vector(mode="list", length=n)
+		for (i in 1:n) glist[[i]] <- d[i] * listw$weights[[i]]
+		sd1 <- 1/sqrt(d)
+		for (i in 1:n) {
+			inb <- listw$neighbours[[i]]
+			icd <- cardnb[i]
+			if (icd > 0) {
+				for (j in 1:icd) {
+					glist[[i]][j] <- sd1[i] * 
+						glist[[i]][j] * sd1[inb[j]]
+				}
+			}
+		}
+		res <- listw
+		res$weights <- glist
+		attr(res$weights, "mode") <- "sim"
+		attr(res$weights, "W") <- TRUE
+		attr(res$weights, "comp") <- attr(listw$weights, "comp")
+		res$style <- "W:sim"
+	} else if (listw$style == "S") {
+		q <- attr(listw$weights, "comp")$q
+		Q <- attr(listw$weights, "comp")$Q
+		eff.n <- attr(listw$weights, "comp")$eff.n
+		glist <- vector(mode="list", length=n)
+		for (i in 1:n) {
+			glist[[i]] <- (Q/eff.n) * listw$weights[[i]]
+			glist[[i]] <- q[i] * glist[[i]]
+		}
+		sq1 <- 1/sqrt(q)
+		for (i in 1:n) {
+			inb <- listw$neighbours[[i]]
+			icd <- cardnb[i]
+			if (icd > 0) {
+				for (j in 1:icd) {
+					glist[[i]][j] <- sq1[i] * 
+						glist[[i]][j] * sq1[inb[j]]
+				}
+				glist[[i]] <- (eff.n/Q) * glist[[i]]
+			}
+		}
+		res <- listw
+		res$weights <- glist
+		attr(res$weights, "mode") <- "sim"
+		attr(res$weights, "S") <- TRUE
+		attr(res$weights, "comp") <- attr(listw$weights, "comp")
+		res$style <- "S:sim"
+	} else stop("Conversion not suitable for this weights style")
+	sym_out <- is.symmetric.glist(res$neighbours, res$weights)
+	if (!sym_out) {
+	    if (attr(sym_out, "d") < .Machine$double.eps ^ 0.5)
+		res <- listw2U(res)
+	    else stop("defective similarity")
+	}
+	res
+}
+

Added: pkg/R/arg.R
===================================================================
--- pkg/R/arg.R	                        (rev 0)
+++ pkg/R/arg.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,11 @@
+`arg` <-
+function (rhopar, v, verbose = FALSE) 
+{
+    vv <-  v$bigG %*% c(rhopar[1], rhopar[1]^2, rhopar[2]) - v$smallg
+    value <- sum(vv^2)
+    if (verbose) 
+        cat("function:", value, "rho:", rhopar[1], "sig2:", 
+            rhopar[2], "\n")
+    value
+}
+


Property changes on: pkg/R/arg.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/arg1.R
===================================================================
--- pkg/R/arg1.R	                        (rev 0)
+++ pkg/R/arg1.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,18 @@
+`arg1` <-
+function (rhopar, v, ss,SS,T, verbose = FALSE) 
+{
+	Ga<-cbind( (1/(T-1))*ss^2,0)
+	Gb<-cbind( 0, SS^2)
+	Gc<-rbind(Ga,Gb)
+	Gamma<-kronecker(Gc,diag(3)) 
+	Gammainv<-solve(Gamma)
+    vv <-  v$GG %*% c(rhopar[1], rhopar[1]^2, rhopar[2], rhopar[3]) - v$gg
+    value <- t(vv)%*% Gammainv %*% vv
+    if (verbose) 
+	cat("function:", value, "rho:", rhopar[1], "sig2:", 
+		rhopar[2], "\n")
+    value
+}
+
+
+


Property changes on: pkg/R/arg1.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/arg2.R
===================================================================
--- pkg/R/arg2.R	                        (rev 0)
+++ pkg/R/arg2.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,32 @@
+`arg2` <-
+function (rhopar, v, ss,SS,T,TW, verbose = FALSE) 
+{
+	Ga<-cbind( (1/(T-1))*ss^2,0)
+	Gb<-cbind( 0, SS^2)
+	Gc<-rbind(Ga,Gb)
+	Gamma<-kronecker(Gc,TW) 
+	Gammainv<-solve(Gamma)
+    vv <-  v$GG %*% c(rhopar[1], rhopar[1]^2, rhopar[2], rhopar[3]) - v$gg
+    value <- t(vv)%*% Gammainv %*% vv
+    if (verbose) 
+	cat("function:", value, "rho:", rhopar[1], "sig2:", 
+		rhopar[2], "\n")
+    value
+}
+
+
+
+`arg3` <-
+function (rhopar, v, ss,T,TW, verbose = FALSE) 
+{
+	Ga<-(1/(T-1))*ss^2
+	Gamma<-as.numeric(Ga)*TW 
+	Gammainv<-solve(Gamma)
+    vv <-  v$bigG %*% c(rhopar[1], rhopar[1]^2, rhopar[2]) - v$smallg
+    value <- t(vv)%*% Gammainv %*% vv
+    if (verbose) 
+	cat("function:", value, "rho:", rhopar[1], "sig2:", 
+		rhopar[2], "\n")
+    value
+}
+


Property changes on: pkg/R/arg2.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/bsjktest.R
===================================================================
--- pkg/R/bsjktest.R	                        (rev 0)
+++ pkg/R/bsjktest.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,5 @@
+`bsjktest` <-
+function(x,...){
+  UseMethod("bsjktest")
+}
+


Property changes on: pkg/R/bsjktest.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/bsjktest.formula.R
===================================================================
--- pkg/R/bsjktest.formula.R	                        (rev 0)
+++ pkg/R/bsjktest.formula.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,34 @@
+`bsjktest.formula` <-
+function(x, data, w, test=c(paste("C",1:3,sep="."),"J"), index=NULL, ...){
+  ## transform data if needed
+
+  if(!is.null(index)) {
+    require(plm)
+    data <- plm.data(data, index)
+    }
+
+  gindex <- dimnames(data)[[2]][1]
+  tindex <- dimnames(data)[[2]][2]
+
+  switch(match.arg(test), C.1 = {
+
+    bsjk = pbsjkSDtest(formula=x, data=data, w=w, index=index, ...)
+
+  }, C.2 = {
+
+    bsjk = pbsjkARtest(formula=x, data=data, w=w, index=index, ...)
+
+  }, C.3 = {
+
+    bsjk = pbsjkREtest(formula=x, data=data, w=w, index=index, ...)
+
+  }, J = {
+
+    bsjk = pbsjkJtest(formula=x, data=data, w=w, index=index, ...)
+
+  })
+
+  return(bsjk)
+
+}
+


Property changes on: pkg/R/bsjktest.formula.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/bsktest.R
===================================================================
--- pkg/R/bsktest.R	                        (rev 0)
+++ pkg/R/bsktest.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,5 @@
+`bsktest` <-
+function(x,...){
+  UseMethod("bsktest")
+}
+


Property changes on: pkg/R/bsktest.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/bsktest.formula.R
===================================================================
--- pkg/R/bsktest.formula.R	                        (rev 0)
+++ pkg/R/bsktest.formula.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,30 @@
+`bsktest.formula` <-
+function(x, data, w, test=c("SLM1","SLM2","LMJOINT","CLMlambda","CLMmu"), index=NULL, ...){
+  
+
+switch(match.arg(test), SLM1 = {
+
+    bsk = slm1test(x, data, index,  w)
+
+  }, SLM2 = {
+
+    bsk = slm2test(x, data, index,  w)
+
+  }, LMJOINT = {
+
+    bsk = LMHtest(x, data, index,  w)
+
+  }, CLMlambda = {
+
+    bsk = clmltest(x, data, index,  w)
+
+  }, CLMmu = {
+
+    bsk = clmmtest(x, data, index,  w)
+
+  })
+
+  return(bsk)
+
+}
+


Property changes on: pkg/R/bsktest.formula.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/bsktest.lm.R
===================================================================
--- pkg/R/bsktest.lm.R	                        (rev 0)
+++ pkg/R/bsktest.lm.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,20 @@
+`bsktest.lm` <-
+function(x, w, index=NULL, test=c("SLM1","SLM2","LMJOINT"), ...){
+	
+	switch(match.arg(test), SLM1 = {
+
+    bsk = slm1test.model(x,w, index)
+
+  }, SLM2 = {
+
+    bsk = slm2test.model(x,w, index )
+
+  }, LMJOINT = {
+
+    bsk = LMHtest.model(x,w, index)
+
+  })
+
+  return(bsk)
+}
+


Property changes on: pkg/R/bsktest.lm.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/bsktest.splm.R
===================================================================
--- pkg/R/bsktest.splm.R	                        (rev 0)
+++ pkg/R/bsktest.splm.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,16 @@
+`bsktest.splm` <-
+function(x, w, index=NULL, test=c("CLMlambda","CLMmu"), ...){
+	
+	switch(match.arg(test), CLMlambda = {
+
+    bsk = clmltest.model(x,w, index)
+
+  }, CLMmu = {
+
+    bsk = clmmtest.model(x,w, index )
+
+  })
+
+  return(bsk)
+}
+


Property changes on: pkg/R/bsktest.splm.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/clmltest.R
===================================================================
--- pkg/R/clmltest.R	                        (rev 0)
+++ pkg/R/clmltest.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,81 @@
+`clmltest` <-
+function(formula, data, index=NULL, listw){
+    ## depends on listw2dgCMatrix.R, REmod.R, spreml.R
+
+	ml <- spreml(formula, data=data, w=listw2mat(listw), errors="re")
+	  if(!is.null(index)) {
+    require(plm)
+    data <- plm.data(data, index)
+    }
+  index <- data[,1]
+  tindex <- data[,2]
+  X<-model.matrix(formula,data=data)
+  y<-model.response(model.frame(formula,data=data))
+  ## reduce index accordingly
+   names(index)<-row.names(data)
+  ind<-index[which(names(index)%in%row.names(X))]
+  tind<-tindex[which(names(index)%in%row.names(X))]
+
+  ## reorder data by cross-sections, then time
+  oo<-order(tind,ind)
+  X<-X[oo,]
+  y<-y[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  ## det. number of groups and df
+  N<-length(unique(ind))
+  k<-dim(X)[[2]]
+  ## det. max. group numerosity
+  T<-max(tapply(X[,1],ind,length))
+  ## det. total number of obs. (robust vs. unbalanced panels)
+  NT<-length(ind)
+	eML<-residuals(ml)
+###maximum likelihood estimation under the null hypothesis that lambda is equal to zero. extract the residuals.
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) ####indicator to get the cross-sectional observations
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) ####indicator to get the time periods observations
+
+	eme<-tapply(eML,inde1,mean)
+	emme<-eML - rep(eme,T)
+	sigmav<-crossprod(eML,emme)/(N*(T-1)) ####estimate of the variance component sigma_v
+	sigma1<-crossprod(eML,rep(eme,T))/N ####estimate of the variance component sigma_1
+	c1<-sigmav/sigma1^2
+	c2<-1/sigmav
+	c1e<-as.numeric(c1)*eML
+	Wst<-listw2dgCMatrix(listw) ###transform the listw object in a sparse matrix
+	Ws<-t(Wst)  ### this is the real W since listw2dgCMatrix generate W'
+	WpsW<-Wst+Ws
+yybis<-function(q){ #### for very big dimension of the data this can be changed looping over the rows and columns of W or either the listw object
+	wq<-(WpsW)%*%q
+	wq<-as.matrix(wq)
+	}
+	Wc1e<-unlist(tapply(eML,inde,yybis)) #### (W'+W)*u
+	sumWc1e<-tapply(Wc1e,inde1,sum)
+	prod1<-as.numeric(c1)*rep(sumWc1e,T)/T
+	prod2<-as.numeric(c2)* (Wc1e - rep(sumWc1e,T)/T)
+	prod<-prod1+prod2
+	D<-1/2*crossprod(eML,prod) ###calculates D in the notation of the paper.
+	W2<-Ws%*%Ws ####generate W^2
+	WW<-crossprod(Ws) ####generate W'*W
+tr<-function(R) sum(diag(R))
+	b<-tr(W2+WW) ###generates b (same notation as the paper)
+	LMl1<-D^2/(((T-1)+as.numeric(sigmav)^2/as.numeric(sigma1)^2)*b) ###conditional LM test for lambda equal to zero
+	LMlstar<-sqrt(LMl1) ###one-sided version
+	statistics<-LMlstar
+  pval <- pnorm(LMlstar, lower.tail=FALSE)
+
+  names(statistics)="LM*-lambda"
+	method<- "Baltagi, Song and Koh LM*-lambda conditional LM test (assuming sigma^2_mu >= 0)"
+	#alt<-"serial corr. in error terms, sub RE and spatial dependence"
+  ##(insert usual htest features)
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(formula), alternative="Spatial autocorrelation")
+  class(RVAL) <- "htest"
+  return(RVAL)
+
+}
+


Property changes on: pkg/R/clmltest.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/clmltest.model.R
===================================================================
--- pkg/R/clmltest.model.R	                        (rev 0)
+++ pkg/R/clmltest.model.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,86 @@
+`clmltest.model` <-
+function(x, listw, index){
+    ## depends on:
+    ## listw2dgCMatrix.R
+    ## REmod.R
+    ## spreml.R
+if(!inherits(x,"splm")) stop("argument should be an object of class splm")
+frm<-x$call
+if(x$type != "random effects ML") stop("argument should be of type random effects ML")
+
+  if(is.null(index))  stop("index should be specified to retrieve information on time and cross-sectional dimentions")
+
+  if(!inherits(listw,"listw")) stop("object w should be of class listw")
+
+
+
+  ind <- index[,1]
+  tind <- index[,2]
+
+if(names(x$coefficients)[1]=="(Intercept)")  X<-data.frame(cbind(rep(1,ncol(x$model)), x$model[,-1]))
+else X<-x$model[,-1]
+  y<-x$model[,1]
+  	eML<-x$residuals
+
+  ## reduce index accordingly
+
+  ## reorder data by cross-sections, then time
+  oo<-order(tind,ind)
+  X<-X[oo,]
+  y<-y[oo]
+
+  ## det. number of groups and df
+  N<-length(unique(ind))
+  k<-dim(X)[[2]]
+  ## det. max. group numerosity
+  T<-max(tapply(X[,1],ind,length))
+  ## det. total number of obs. (robust vs. unbalanced panels)
+  NT<-length(ind)
+###maximum likelihood estimation under the null hypothesis that lambda is equal to zero. extract the residuals.
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) ####indicator to get the cross-sectional observations
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) ####indicator to get the time periods observations
+
+	eme<-tapply(eML,inde1,mean)
+	emme<-eML - rep(eme,T)
+	sigmav<-crossprod(eML,emme)/(N*(T-1)) ####estimate of the variance component sigma_v
+	sigma1<-crossprod(eML,rep(eme,T))/N ####estimate of the variance component sigma_1
+	c1<-sigmav/sigma1^2
+	c2<-1/sigmav
+	c1e<-as.numeric(c1)*eML
+	Wst<-listw2dgCMatrix(listw) ###transform the listw object in a sparse matrix
+	Ws<-t(Wst)  ### this is the real W since listw2dgCMatrix generate W'
+	WpsW<-Wst+Ws
+yybis<-function(q){ #### for very big dimension of the data this can be changed looping over the rows and columns of W or either the listw object
+	wq<-(WpsW)%*%q
+	wq<-as.matrix(wq)
+	}
+	Wc1e<-unlist(tapply(eML,inde,yybis)) #### (W'+W)*u
+	sumWc1e<-tapply(Wc1e,inde1,sum)
+	prod1<-as.numeric(c1)*rep(sumWc1e,T)/T
+	prod2<-as.numeric(c2)* (Wc1e - rep(sumWc1e,T)/T)
+	prod<-prod1+prod2
+	D<-1/2*crossprod(eML,prod) ###calculates D in the notation of the paper.
+	W2<-Ws%*%Ws ####generate W^2
+	WW<-crossprod(Ws) ####generate W'*W
+tr<-function(R) sum(diag(R))
+	b<-tr(W2+WW) ###generates b (same notation as the paper)
+	LMl1<-D^2/(((T-1)+as.numeric(sigmav)^2/as.numeric(sigma1)^2)*b) ###conditional LM test for lambda equal to zero
+	LMlstar<-sqrt(LMl1) ###one-sided version
+	statistics<-LMlstar
+  pval <- pnorm(LMlstar, lower.tail=FALSE)
+
+  names(statistics)="LM*-lambda"
+	method<- "Baltagi, Song and Koh LM*-lambda conditional LM test (assuming sigma^2_mu >= 0)"
+	#alt<-"serial corr. in error terms, sub RE and spatial dependence"
+  ##(insert usual htest features)
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(frm), alternative="Spatial autocorrelation")
+  class(RVAL) <- "htest"
+  return(RVAL)
+
+}
+


Property changes on: pkg/R/clmltest.model.R
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/R/clmmtest.R
===================================================================
--- pkg/R/clmmtest.R	                        (rev 0)
+++ pkg/R/clmmtest.R	2010-07-13 16:52:02 UTC (rev 78)
@@ -0,0 +1,104 @@
+`clmmtest` <-
+function(formula, data, index=NULL, listw){
+    ## depends on listw2dgCMatrix.R, spfeml.R
+
+	ml <- spfeml(formula, data=data, index=index, listw, model="error", effects="pooled")
+
+	  if(!is.null(index)) {
+    require(plm)
+    data <- plm.data(data, index)
+    }
+  index <- data[,1]
+  tindex <- data[,2]
+  X<-model.matrix(formula,data=data)
+  y<-model.response(model.frame(formula,data=data))
+  ## reduce index accordingly
+   names(index)<-row.names(data)
+  ind<-index[which(names(index)%in%row.names(X))]
+  tind<-tindex[which(names(index)%in%row.names(X))]
+
+  ## reorder data by cross-sections, then time
+  oo<-order(tind,ind)
+  X<-X[oo,]
+  y<-y[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  ## det. number of groups and df
+  N<-length(unique(ind))
+  k<-dim(X)[[2]]
+  ## det. max. group numerosity
+  T<-max(tapply(X[,1],ind,length))
+  ## det. total number of obs. (robust vs. unbalanced panels)
+  NT<-length(ind)
+
+###maximum likelihood estimation under the null hypothesis that lambda is equal to zero. extract the residuals.
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) ####indicator to get the cross-sectional observations
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) ####indicator to get the time periods observations
+
+	lambda<-ml$spat.coef
+	#print(lambda)
+	eML<-residuals(ml)
+#	print(length(eML))
+ 	Wst<-listw2dgCMatrix(listw) ###transform the listw object in a sparse matrix
+	Ws<-t(Wst)  ### this is the real W since listw2dgCMatrix generate W'
+	B<- -lambda*Ws
+	diag(B)<- 1
+	BpB<-crossprod(B)
+	BpBi<- solve(BpB)
+vc<-function(R) {
+	BBu<-BpBi%*%R
+	BBu<-as.matrix(BBu)
+	}
+	eme<-unlist(tapply(eML,inde,vc))
+	sigmav2<-crossprod(eML,eme)/(N*(T-1)) ####estimate of the variance component sigma_v
+	sigmav4<-sigmav2^2
+
+tr<-function(R) sum(diag(R))
+	trBpB<-tr(BpB)
+	BpB2<-BpB%*%BpB
+yybis<-function(q){ #### for very big dimension of the data this can be changed looping over the rows and columns of W or either the listw object
+	wq<-rep(q,T)
+	tmp<-wq%*%eML
+					}
+	BBu<-apply(BpB2,1,yybis)
+	BBu<-rep(BBu,T)
+	upBBu<-crossprod(eML,BBu)
+	Dmu<--((T/(2*sigmav2))*trBpB) + ((1/(2*sigmav4))*upBBu)
+	WpB<-Wst%*%B
+	BpW<-t(B)%*%Ws
+	WpBplBpW <-WpB + BpW
+	G<-WpBplBpW %*% BpBi
+	e<-tr(BpB2)
+	d<-tr(WpBplBpW)
+	h<-trBpB
+	g<-tr(G)
+	c<-tr(G%*%G)
+	#print(c(e,d,h,g,c))
+	NUM<- ((2*sigmav4)/T)*((N*sigmav4*c)-(sigmav4*g^2))   ###equation 2.30 in the paper
+	DENft<- NT*sigmav4*e*c
+	DENst<- N*sigmav4*d^2
+	DENtt<- T*sigmav4*g^2 * e
+	DENfot<- 2* sigmav4 *g*h*d
+	DENfit<- sigmav4*h^2* c
+	DEN<- DENft - DENst - DENtt + DENfot - DENfit
+	LMmu <- Dmu^2*NUM / DEN
+	LMmustar<- sqrt(LMmu)
+	statistics<-LMmustar
+  pval <- pnorm(LMmustar, lower.tail=FALSE)
+
+  names(statistics)="LM*-mu"
+	method<- "Baltagi, Song and Koh LM*- mu conditional LM test (assuming lambda may or may not be = 0)"
+	#alt<-"serial corr. in error terms, sub RE and spatial dependence"
+  ##(insert usual htest features)
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(formula), alternative="Random regional effects")
+  class(RVAL) <- "htest"
+  return(RVAL)
+
+}
+


Property changes on: pkg/R/clmmtest.R
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/splm -r 78


More information about the Splm-commits mailing list