[Splm-commits] r98 - / pkg pkg/R pkg/chm pkg/data pkg/inst pkg/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 13 18:40:19 CEST 2011
Author: gpiras
Date: 2011-04-13 18:40:19 +0200 (Wed, 13 Apr 2011)
New Revision: 98
Added:
pkg/
pkg/ChangeLog
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/
pkg/R/.Rapp.history
pkg/R/LMHtest.R
pkg/R/LMHtest.model.R
pkg/R/REmod.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/fixed_effects.R
pkg/R/ivplm.b2sls.R
pkg/R/ivplm.ec2sls.R
pkg/R/ivplm.g2sls.R
pkg/R/ivplm.w2sls.R
pkg/R/ivsplm.R
pkg/R/likelihoodsFE.R
pkg/R/listw2dgCMatrix.R
pkg/R/lrtest.splm.R
pkg/R/olsmod.R
pkg/R/pbsjkARtest.R
pkg/R/pbsjkJtest.R
pkg/R/pbsjkREtest.R
pkg/R/pbsjkSDtest.R
pkg/R/print.splm.R
pkg/R/print.summary.splm.R
pkg/R/sarREmod.R
pkg/R/sarem2REmod.R
pkg/R/saremREmod.R
pkg/R/saremmod.R
pkg/R/saremsrREmod.R
pkg/R/saremsrmod.R
pkg/R/sarmod.R
pkg/R/sarsrREmod.R
pkg/R/sarsrmod.R
pkg/R/sem2REmod.R
pkg/R/semREmod.R
pkg/R/semmod.R
pkg/R/semsrREmod.R
pkg/R/semsrmod.R
pkg/R/slm1test.R
pkg/R/slm1test.model.R
pkg/R/slm2test.R
pkg/R/slm2test.model.R
pkg/R/spfeml.R
pkg/R/spgm.R
pkg/R/sphtest.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/tss.R
pkg/R/utilities_GM.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/listw2dgCMatrix.Rd
pkg/man/print.effects.splm.Rd
pkg/man/print.splm.Rd
pkg/man/spfeml.Rd
pkg/man/spgm.Rd
pkg/man/splm-package.Rd
pkg/man/spreml.Rd
pkg/man/spsegm.Rd
pkg/man/spseml.Rd
pkg/man/summary.effects.splm.Rd.old
pkg/man/summary.splm.Rd
pkg/man/usaww.Rd
pkg/man/write.effects.splm.Rd
Log:
upload version 0.8.2
Added: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog (rev 0)
+++ pkg/ChangeLog 2011-04-13 16:40:19 UTC (rev 98)
@@ -0,0 +1,23 @@
+Changes in Version 0.8-01
+ o added spgm: general function that deals with all the GM estimators
+ o added the methodologies in Mutl and Pfaffermeyer (2011) and Piras (2011)
+ for the estimation of the GM models sperrorgm and spsarargm
+ o includes the following estimators: ivplm.w2sls, ivplm.b2sls, ivplm.ec2sls, ivplm.g2sls
+ along with ivsplm that is the wrapper to use them.
+
+Changes in Version 0.2-04
+ o dependency changed from kinship to bdsmatrix; removed require(kinship) from all functions
+
+Changes in Version 0.2-02
+ o spfeml: Added methods for Jacobian
+
+
+Changes in Version 0.2-01
+
+
+ o spregm: modified to allow for additional endogenous variables and lag of the dependent variable
+ o Added spfegm
+ o Added spseml
+ o spsegm: improved substantially and now reads a list of formulas.
+
+
Added: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION (rev 0)
+++ pkg/DESCRIPTION 2011-04-13 16:40:19 UTC (rev 98)
@@ -0,0 +1,10 @@
+Package: splm
+Title: Econometric Models for Spatial Panel Data
+Version: 0.8-02
+Date: 2011-04-13
+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 (>= 2.11.1), MASS, nlme, spdep, plm, Matrix, bdsmatrix, spam
+License: GPL-2
+LazyLoad: yes
Property changes on: pkg/DESCRIPTION
___________________________________________________________________
Added: svn:executable
+
Added: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE (rev 0)
+++ pkg/NAMESPACE 2011-04-13 16:40:19 UTC (rev 98)
@@ -0,0 +1,25 @@
+importFrom(stats, model.matrix, model.response, aggregate, effects)
+import(nlme)
+import(spdep)
+import(Matrix)
+importFrom(bdsmatrix,bdsmatrix)
+importFrom(MASS,ginv)
+
+export(bsjktest, bsktest,
+effects.splm, print.effects.splm, write.effects.splm,
+print.splm, spfeml, spgm, spreml, summary.splm,
+spseml, spsegm, lrtest.splm, sphtest, listw2dgCMatrix)
+
+
+
+S3method(print,splm)
+S3method(print,summary.splm)
+S3method(bsjktest,formula)
+S3method(effects,splm)
+S3method(print,effects.splm)
+S3method(bsktest,formula)
+S3method(bsktest,lm)
+S3method(bsktest, splm)
+S3method(sphtest,formula)
+S3method(sphtest, splm)
+
Property changes on: pkg/NAMESPACE
___________________________________________________________________
Added: svn:executable
+
Added: pkg/R/.Rapp.history
===================================================================
--- pkg/R/.Rapp.history (rev 0)
+++ pkg/R/.Rapp.history 2011-04-13 16:40:19 UTC (rev 98)
@@ -0,0 +1,16 @@
+system.time()
+help(time)
+proc.time()
+date()
+date()[4]
+sys.time()
+Sys.time()
+p1<-Sys.time()
+p2<-Sys.time()
+p2-p1
+lagsarlm
+library(spdep)
+lagsarlm
+proc.time()
+help(Matrix)
+ p0
Added: pkg/R/LMHtest.R
===================================================================
--- pkg/R/LMHtest.R (rev 0)
+++ pkg/R/LMHtest.R 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -0,0 +1,157 @@
+REmod <-
+function (X, y, ind, tind, n, k, t, nT, w, w2, coef0 = rep(0, 2),
+ hess = FALSE, trace = trace, x.tol = 1.5e-18, rel.tol = 1e-15,
+ ...)
+{
+
+## optimizing version 1:
+ ##
+ ## exploit ordering reversal
+ ## and bdsmatrix functions as in ssrmod, sarsrmod, sarREmod...
+ ##
+ ## a) lag y etc.
+ ## b) reverse ordering and exploit bds nature of vcov(e)
+ ##
+ ## maybe exploit analytical inverse of the submatrix block (gains on
+ ## large-N problems)?? but how likely is it to have laaaarge T?
+
+ ## extensive function rewriting, Giovanni Millo 04/10/2010
+ ## structure:
+ ## a) specific part
+ ## - set names, bounds and initial values for parms
+ ## - define building blocks for likelihood and GLS as functions of parms
+ ## - define likelihood
+ ## b) generic part(independent from ll.c() and #parms)
+ ## - fetch covariance parms from max lik
+ ## - calc last GLS step
+ ## - fetch betas
+ ## - calc final covariances
+ ## - make list of results
+
+ ## change this to 'bdsmatrix'
+ #require(kinship)
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("phi")
+
+ ## initialize values for optimizer
+ myparms0 <- coef0
+ ## set bounds for optimizer
+ lower.bounds <- c(1e-08)
+ upper.bounds <- c(1e+09)
+
+ ## rearranging module
+ ## save this for eventually re-rearranging output
+ oo.0 <- order(tind, ind)
+ ## reorder as stacked time series, as in std. panels
+ oo <- order(ind, tind)
+ X <- X[oo, ]
+ y <- y[oo]
+ #wy <- wy[oo]
+ ind <- ind[oo]
+ tind <- tind[oo]
+
+ ## modules for likelihood
+ bSigma <- function(phipsi, n, t, w) {
+ ## single block of the original *scaled* covariance
+ ## maintain w for homogeneity with generic part
+ Jt <- matrix(1, ncol = t, nrow = t)
+ It <- diag(1, t)
+ ## retrieve parms
+ phi <- phipsi[1]
+ ## psi not used: here passing 2 parms, but works anyway
+ ## because psi is last one
+ ## calc inverse
+ bSigma <- phi * Jt + It
+ bSigma
+ }
+ detSigma <- function(phi, n, t) {
+ detSigma <- -n/2 * log(t * phi + 1)
+ detSigma
+ }
+ fullSigma <- function(phipsi, n, t, w) {
+ sigma.i <- bSigma(phipsi, n, t, w)
+ fullSigma <- bdsmatrix(rep(t, n), rep(as.numeric(sigma.i),
+ n))
+ fullSigma
+ }
+
+ ## likelihood function, both steps included
+ ll.c <- function(phipsi, y, X, n, t, w, w2, wy) {
+ ## retrieve parms
+ phi <- phipsi[1]
+ ## calc inverse sigma
+ sigma <- fullSigma(phipsi, n, t, w)
+ ## do GLS step to get e, s2e
+ glsres <- GLSstepBDS(X, y, sigma)
+ e <- glsres[["ehat"]]
+ s2e <- glsres[["sigma2"]]
+ ## calc ll
+ due <- detSigma(phi, n, t)
+ tre <- -n * t/2 * log(s2e)
+ quattro <- -1/(2 * s2e) * crossprod(e, solve(sigma, e))
+ const <- -(n * t)/2 * log(2 * pi)
+ ll.c <- const + due + tre + quattro
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## GLS step function for bdsmatrices
+ GLSstepBDS <- function(X, y, sigma) {
+ b.hat <- solve(crossprod(X, solve(sigma, X)), crossprod(X,
+ solve(sigma, y)))
+ ehat <- y - X %*% b.hat
+ sigma2ehat <- crossprod(ehat, solve(sigma, ehat))/(n * t)
+ return(list(betahat=b.hat, ehat=ehat, sigma2=sigma2ehat))
+ }
+
+ ## lag y unneeded here, keep parm for compatibility
+ wy <- NULL
+
+ ## max likelihood
+ optimum <- nlminb(start = myparms0, objective = ll.c,
+ gradient = NULL, hessian = NULL,
+ y = y, X = X, n = n, t = t, w = w, w2 = w2, wy = wy,
+ scale = 1, control = list(x.tol = x.tol,
+ rel.tol = rel.tol, trace = trace),
+ lower = lower.bounds, upper = upper.bounds)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms
+ myparms <- optimum$par
+
+ ## one last GLS step at optimal vcov parms
+ sigma <- fullSigma(myparms, n, t, w)
+ beta <- GLSstepBDS(X, y, sigma)
+
+ ## final vcov(beta)
+ covB <- as.numeric(beta[[3]]) *
+ solve(crossprod(X, solve(sigma, X)))
+
+ ## final vcov(errcomp)
+ covTheta <- solve(-fdHess(myparms, function(x) -ll.c(x,
+ y, X, n, t, w, w2, wy))$Hessian) # lag-specific line: wy
+ nvcovpms <- length(nam.errcomp)
+ covAR <- NULL
+ covPRL <- covTheta[1:nvcovpms, 1:nvcovpms, drop=FALSE]
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- NULL
+ errcomp <- myparms[which(nam.errcomp!="psi")]
+ names(betas) <- nam.beta
+ names(errcomp) <- nam.errcomp[which(nam.errcomp!="psi")]
+
+ dimnames(covB) <- list(nam.beta, nam.beta)
+ dimnames(covPRL) <- list(names(errcomp), names(errcomp))
+
+ ## result
+ RES <- list(betas = betas, arcoef=arcoef, errcomp = errcomp,
+ covB = covB, covAR=covAR, covPRL = covPRL, ll = myll)
+
+ return(RES)
+}
Property changes on: pkg/R/REmod.R
___________________________________________________________________
Added: svn:executable
+
Added: pkg/R/bsjktest.R
===================================================================
--- pkg/R/bsjktest.R (rev 0)
+++ pkg/R/bsjktest.R 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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 2011-04-13 16:40:19 UTC (rev 98)
@@ -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
___________________________________________________________________
Added: svn:executable
+
Added: pkg/R/clmmtest.model.R
===================================================================
--- pkg/R/clmmtest.model.R (rev 0)
+++ pkg/R/clmmtest.model.R 2011-04-13 16:40:19 UTC (rev 98)
@@ -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
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/splm -r 98
More information about the Splm-commits
mailing list