[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