[Splm-commits] r13 - in pkg: . R chm data inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Sep 29 23:27:04 CEST 2009


Author: gpiras
Date: 2009-09-29 23:27:04 +0200 (Tue, 29 Sep 2009)
New Revision: 13

Added:
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/
   pkg/R/Ggsararsp.R
   pkg/R/LMHtest.R
   pkg/R/LMHtest.model.R
   pkg/R/REmod.R
   pkg/R/arg.R
   pkg/R/arg1.R
   pkg/R/arg2.R
   pkg/R/bsjktest.R
   pkg/R/bsjktest.formula.R
   pkg/R/bsjktest.splm.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.splm.R
   pkg/R/feerror.R
   pkg/R/felag.R
   pkg/R/fs.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/spfeml.R
   pkg/R/splaglm.R
   pkg/R/spregm.R
   pkg/R/spreml.R
   pkg/R/spsegm.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/vcov.splm.R
   pkg/chm/
   pkg/chm/00Index.html
   pkg/chm/Rchm.css
   pkg/chm/bsjktest.html
   pkg/chm/bsktest.html
   pkg/chm/effects.splm.html
   pkg/chm/logo.jpg
   pkg/chm/print.splm.html
   pkg/chm/spfeml.html
   pkg/chm/splm-package.html
   pkg/chm/splm.hhp
   pkg/chm/splm.toc
   pkg/chm/spregm.html
   pkg/chm/spreml.html
   pkg/chm/spsegm.html
   pkg/chm/summary.splm.html
   pkg/data/
   pkg/data/usaww.rda
   pkg/inst/
   pkg/inst/doc/
   pkg/man/
   pkg/man/bsjktest.Rd
   pkg/man/bsktest.Rd
   pkg/man/effects.splm.Rd
   pkg/man/print.effects.splm.Rd
   pkg/man/print.splm.Rd
   pkg/man/spfeml.Rd
   pkg/man/splm-package.Rd
   pkg/man/spregm.Rd
   pkg/man/spreml.Rd
   pkg/man/spsegm.Rd
   pkg/man/summary.effects.splm.Rd.old
   pkg/man/summary.splm.Rd
   pkg/man/usaww.Rd
   pkg/man/write.effects.splm.Rd
Log:
initial import

Added: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	                        (rev 0)
+++ pkg/DESCRIPTION	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,10 @@
+Package: splm
+Title: Econometric Models for Spatial Panel Data
+Version: 0.1-3
+Date: 2009-09-09
+Author: Giovanni Millo <giovanni.millo at generali.com>, Gianfranco Piras <gpiras at mac.com>
+Maintainer: Giovanni Millo <giovanni.millo at generali.com>
+Description: ML and GM estimation and diagnostic testing of econometric models for spatial panel data.
+Depends: R (>= 1.8.0), MASS, nlme, spdep, plm, Matrix, kinship
+License: GPL-2
+LazyLoad: yes


Property changes on: pkg/DESCRIPTION
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	                        (rev 0)
+++ pkg/NAMESPACE	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,21 @@
+importFrom(stats, model.matrix, model.response, aggregate, effects)
+import(nlme)
+import(spdep)
+import(Matrix)
+importFrom(kinship,bdsmatrix)
+importFrom(MASS,ginv)
+
+export(bsjktest, bsktest, effects.splm, print.effects.splm, write.effects.splm, 
+print.splm,spfeml,spregm,spreml,spsegm,summary.splm)
+
+
+
+S3method(print,splm)
+S3method(print,summary.splm)
+S3method(bsjktest,formula)
+S3method(effects,splm)
+S3method(print,effects.splm)
+S3method(write,effects.splm)
+S3method(bsktest,formula)
+S3method(bsktest,lm)
+S3method(bsktest, splm)


Property changes on: pkg/NAMESPACE
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/Ggsararsp.R
===================================================================
--- pkg/R/Ggsararsp.R	                        (rev 0)
+++ pkg/R/Ggsararsp.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/LMHtest.R
===================================================================
--- pkg/R/LMHtest.R	                        (rev 0)
+++ pkg/R/LMHtest.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/LMHtest.model.R
===================================================================
--- pkg/R/LMHtest.model.R	                        (rev 0)
+++ pkg/R/LMHtest.model.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/REmod.R
===================================================================
--- pkg/R/REmod.R	                        (rev 0)
+++ pkg/R/REmod.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/arg.R
===================================================================
--- pkg/R/arg.R	                        (rev 0)
+++ pkg/R/arg.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/arg1.R
===================================================================
--- pkg/R/arg1.R	                        (rev 0)
+++ pkg/R/arg1.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,16 @@
+`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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/arg2.R
===================================================================
--- pkg/R/arg2.R	                        (rev 0)
+++ pkg/R/arg2.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,16 @@
+`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
+}
+


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

Added: pkg/R/bsjktest.R
===================================================================
--- pkg/R/bsjktest.R	                        (rev 0)
+++ pkg/R/bsjktest.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,5 @@
+`bsjktest` <-
+function(x,...){
+  UseMethod("bsjktest")
+}
+


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

Added: pkg/R/bsjktest.formula.R
===================================================================
--- pkg/R/bsjktest.formula.R	                        (rev 0)
+++ pkg/R/bsjktest.formula.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/bsjktest.splm.R
===================================================================
--- pkg/R/bsjktest.splm.R	                        (rev 0)
+++ pkg/R/bsjktest.splm.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,5 @@
+`bsjktest.splm` <-
+function(x, x2, ...){
+## make this a function testing SD residuals on a panelmodel?? es. BSK test
+}
+


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

Added: pkg/R/bsktest.R
===================================================================
--- pkg/R/bsktest.R	                        (rev 0)
+++ pkg/R/bsktest.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,5 @@
+`bsktest` <-
+function(x,...){
+  UseMethod("bsktest")
+}
+


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

Added: pkg/R/bsktest.formula.R
===================================================================
--- pkg/R/bsktest.formula.R	                        (rev 0)
+++ pkg/R/bsktest.formula.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/bsktest.lm.R
===================================================================
--- pkg/R/bsktest.lm.R	                        (rev 0)
+++ pkg/R/bsktest.lm.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/bsktest.splm.R
===================================================================
--- pkg/R/bsktest.splm.R	                        (rev 0)
+++ pkg/R/bsktest.splm.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/clmltest.R
===================================================================
--- pkg/R/clmltest.R	                        (rev 0)
+++ pkg/R/clmltest.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/clmltest.model.R
===================================================================
--- pkg/R/clmltest.model.R	                        (rev 0)
+++ pkg/R/clmltest.model.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/clmmtest.R
===================================================================
--- pkg/R/clmmtest.R	                        (rev 0)
+++ pkg/R/clmmtest.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -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
___________________________________________________________________
Name: svn:executable
   + 

Added: pkg/R/clmmtest.model.R
===================================================================
--- pkg/R/clmmtest.model.R	                        (rev 0)
+++ pkg/R/clmmtest.model.R	2009-09-29 21:27:04 UTC (rev 13)
@@ -0,0 +1,105 @@
+`clmmtest.model` <-
+function(x, listw, index){
+    ## depends on listw2dgCMatrix.R, spfeml.R
+
+if(!inherits(x,"splm")) stop("argument should be an object of class splm")
+frm<-x$call
+if(x$type != "fixed effects error") 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
+
+	lambda<-x$spat.coef
+	#print(lambda)
+#	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)"
[TRUNCATED]

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


More information about the Splm-commits mailing list