[Splm-commits] r102 - / pkg pkg/R pkg/chm pkg/data pkg/man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Apr 20 23:39:05 CEST 2011


Author: gpiras
Date: 2011-04-20 23:39:04 +0200 (Wed, 20 Apr 2011)
New Revision: 102

Added:
   pkg/
   pkg/ChangeLog
   pkg/DESCRIPTION
   pkg/NAMESPACE
   pkg/R/
   pkg/R/.Rapp.history
   pkg/R/REmod.R
   pkg/R/bsktest.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/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/spfeml.R
   pkg/R/spgm.R
   pkg/R/sphtest.R
   pkg/R/spreml.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/man/
   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/sphettest.Rd
   pkg/man/splm-package.Rd
   pkg/man/spreml.Rd
   pkg/man/summary.effects.splm.Rd.old
   pkg/man/summary.splm.Rd
   pkg/man/usaww.Rd
   pkg/man/write.effects.splm.Rd
Log:
Deleted simultaneous equations and test bsjk, added hausman test

Added: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog	                        (rev 0)
+++ pkg/ChangeLog	2011-04-20 21:39:04 UTC (rev 102)
@@ -0,0 +1,25 @@
+Changes in Version 0.8-01
+  o added the spatial hausman test
+  o deleted bsjk test and corrected the bsk tests
+  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. 
+ 
+  


Property changes on: pkg/ChangeLog
___________________________________________________________________
Added: svn:executable
   + 

Added: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	                        (rev 0)
+++ pkg/DESCRIPTION	2011-04-20 21:39:04 UTC (rev 102)
@@ -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-20 21:39:04 UTC (rev 102)
@@ -0,0 +1,23 @@
+importFrom(stats, model.matrix, model.response, aggregate, effects)
+import(nlme)
+import(spdep)
+import(Matrix)
+importFrom(bdsmatrix,bdsmatrix)
+importFrom(MASS,ginv)
+
+export(bsktest, sphtest,
+effects.splm, print.effects.splm, write.effects.splm, 
+print.splm, spfeml, spgm, spreml, summary.splm, lrtest.splm, sphtest, listw2dgCMatrix)
+
+
+
+S3method(print, splm)
+S3method(print, summary.splm)
+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-20 21:39:04 UTC (rev 102)
@@ -0,0 +1,371 @@
+elag_ML_eigen <- lagsarlm(eform, data = elect80, listw = elw, method = "eigen")
+summary(elag_ML_eigen)
+system.time(elag_ML_eigen <- lagsarlm(eform, data = elect80, listw = elw, method = "eigen"))
+system.time(elag_ML_Matrix <- lagsarlm(eform, data = elect80, listw = elw,method = "Matrix"))
+Matrix()
+system.time(elag_ML_Matrix <- lagsarlm(eform, data = elect80, listw = elw,method = "Matrix"))
+summary(elag_ML_Matrix)
+eW <- as(as_dgRMatrix_listw(elw), "CsparseMatrix")
+eW
+set.seed(100831)
+etr_MC <- trW(eW, m = 24, type = "MC")
+etr_mult <- trW(eW, m = 24, type = "mult")
+eWs <- spdep:::listw2U_Matrix(spdep:::similar.listw_Matrix(elw))
+etr_mom <- trW(eWs, m = 24, type = "moments")
+cl <- makeSOCKcluster(4)
+library(snow)
+cl <- makeSOCKcluster(4)
+set.ClusterOption(cl)
+etr_mom1 <- trW(eWs, m = 24, type = "moments")
+stopCluster(cl)
+all.equal(etr_mom, etr_mom1, check.attributes = FALSE)
+system.time(etr_MC <- trW(eW, m = 24, type = "MC"))
+system.time(etr_mult <- trW(eW, m = 24, type = "mult"))
+system.time(etr_mom <- trW(eWs, m = 24, type = "moments"))
+eWs <- spdep:::listw2U_Matrix(spdep:::similar.listw_Matrix(elw))
+system.time(etr_mom <- trW(eWs, m = 24, type = "moments"))
+eW <- as(as_dgRMatrix_listw(elw), "CsparseMatrix")#
+set.seed(100831)#
+#
+system.time(etr_MC <- trW(eW, m = 24, type = "MC"))#
+system.time(etr_mult <- trW(eW, m = 24, type = "mult"))#
+#
+eWs <- spdep:::listw2U_Matrix(spdep:::similar.listw_Matrix(elw))
+system.time(etr_mom <- trW(eWs, m = 24, type = "moments"))
+etr_mom
+tr.W
+trW
+system.time(etr_mom <- trW(eWs, m = 24, type = "moments")
+)
+etr_mom <- trW(eWs, m = 24, type = "moments")
+summary.connection
+eW <- as(as_dgRMatrix_listw(elw), "CsparseMatrix")
+set.seed(100831)
+system.time(etr_MC <- trW(eW, m = 24, type = "MC"))#
+system.time(etr_mult <- trW(eW, m = 24, type = "mult"))
+eWs <- spdep:::listw2U_Matrix(spdep:::similar.listw_Matrix(elw))
+eWs
+system.time(etr_mom <- trW(eWs, m = 24, type = "moments"))
+etr_mom <- trW(eWs, m = 24, type = "moments")
+2007 - 1963
+map('state', "west virginia", fill = TRUE, col="#FFCC00", bg="#003366")
+par(mfrow=c(2,2))
+map('state', "west virginia", fill = TRUE, col="#FFCC00", bg="#003366")
+library(maps)
+install.packages("maps")
+library(maps)
+map('state', "west virginia", fill = TRUE, col="#FFCC00", bg="#003366")
+map('state', "west virginia", fill = TRUE, col="#003366", bg="#FFCC00")
+map('state', "west virginia", fill = TRUE, col="blue", bg="gold")
+map('state', "west virginia", fill = TRUE, col="gold", bg="blue")
+map('state', "west virginia")
+par(mfrow=c(2,2))#
+map('state', "west virginia", fill = TRUE, col="#FFCC00", bg="#003366")#
+map('state', "west virginia", fill = TRUE, col="#003366", bg="#FFCC00")#
+map('state', "west virginia", fill = TRUE, col="blue", bg="gold")#
+map('state', "west virginia", fill = TRUE, col="gold", bg="blue")
+col="#FFCC00"
+map('state', "west virginia", fill = TRUE, col="#FFCC00", bg="#003366")
+library(splm)
+help(spfeml)
+spfeml
+help(sphtest)
+help(spreml)
+data(Produc, package = "Ecdat")
+data(usaww)
+Produc <- Produc[Produc$year<1975, ]
+fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
+sarargm<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "fixed", control = list())
+summary(sarargm)
+sarargm<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "random", control = list())
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "fixed", control = list())#
+summary(sarargmfe)#
+#
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "random", control = list())#
+summary(sarargmre)
+coef(sarargmfe)
+coef(sarargmre)
+vcov(sarargmfe)
+sarargmfe$vcov
+sarargmre$vcov
+spgm
+library(splm)
+data(Produc, package = "Ecdat")#
+data(usaww)
+Produc <- Produc[Produc$year<1975, ]#
+fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "fixed", control = list())#
+summary(sarargmfe)
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "random", control = list())#
+summary(sarargmre)
+coef(sarargmfe)
+sarargmfe$vcov
+coef(sarargmre)
+sarargmre$vcov
+sarargmre$legacy
+sarargmfe$legacy
+all.equal(sarargmfe$legacy, sarargmre$legacy)
+sarargmre$effects
+sarargmre$ef.sph
+x$ef.sph == x2$ef.sph
+sarargmre$ef.sph == sarargmfe$ef.sph
+names(coef(sarargmre))
+match(names(coef(sarargmre)), names(coef(sarargmre)) )
+names(coef(sarargmre))
+ names(coef(sarargmre))
+match(names(coef(sarargmfe)), names(coef(sarargmre)) )
+x$ef.sph
+names(coef(sarargmre$ef.sph
+)
+sarargmre$ef.sph
+match("random", c(sarargmfe$ef.sph, sarargmre$ef.sph))
+    tc <- match(names(coef(xwith)), names(coef(xbetw)) )
+x <- sarargmfe
+x2 <- sarargmre
+if (!all.equal(x$legacy, x2$legacy)) stop("The model are different")
+if(x$ef.sph == x2$ef.sph) stop("Effects should be different")
+    ran <- match("random", c(x$ef.sph, x2$ef.sph))
+    ran
+	xwith <- x#
+	xbetw <- x2
+    tc <- match(names(coef(xwith)), names(coef(xbetw)) )
+    coef.wi <- coef(xwith)
+    coef.re <- coef(xbetw)[tc]
+    vcov.wi <- xwith$vcov
+    vcov.re <- xbetw$vcov[tc,tc]
+    coef.wi
+    coef.re
+    vcov.wi
+    vcov.re
+    dbeta <- coef.wi - coef.re
+    dbeta
+    df <- length(dbeta)
+    dvcov <- vcov.re - vcov.wi
+    stat <- abs(t(dbeta) %*% solve(dvcov) %*% dbeta)
+    pval <- pchisq(stat, df = df, lower.tail = FALSE)
+    names(stat) <- "chisq"
+#
+    parameter <- df
+#
+    names(parameter) <- "df"
+#
+    data.name <- paste(deparse(x$call$formula))
+    alternative <- "one model is inconsistent"
+#
+    res <- list(statistic = stat, p.value = pval, parameter = parameter,
+#
+        method = "Hausman test for spatial models", data.name = data.name, alternative = alternative)
+#
+    class(res) <- "htest"
+    res
+library(splm)#
+data(Produc, package = "Ecdat")#
+data(usaww)#
+Produc <- Produc[Produc$year<1975, ]#
+fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "fixed", control = list())#
+summary(sarargmfe)
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "random", control = list())#
+summary(sarargmre)
+sphtest(sarargmfe, sarargmre)
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= FALSE, spatial.error=TRUE, effects = "fixed", control = list())#
+summary(sarargmfe)#
+#
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= FALSE, spatial.error=TRUE, effects = "random", control = list())#
+summary(sarargmre)
+sphtest(sarargmfe, sarargmre)
+    tc <- match(names(coef(xwith)), names(coef(xbetw)) )
+if (!all.equal(x$legacy, x2$legacy)) stop("The model are different")#
+if(x$ef.sph == x2$ef.sph) stop("Effects should be different")
+x<-sarargmfe
+x2<-sarargmre
+if (!all.equal(x$legacy, x2$legacy)) stop("The model are different")#
+if(x$ef.sph == x2$ef.sph) stop("Effects should be different")
+    ran <- match("random", c(x$ef.sph, x2$ef.sph))
+    ran
+	xwith <- x2#
+	xbetw <- x
+	xwith <- x#
+	xbetw <- x2
+	xwith
+	xbetw
+    tc <- match(names(coef(xwith)), names(coef(xbetw)) )
+    tc
+names(coef(xwith)
+)
+names(coef(xbetw))
+xbetw <- x2
+xbetw
+names(coef(xbetw))
+class(coef(xbetw))
+names(coef(as.numeric(xbetw)))
+as.numeric(coef(xbetw))
+names(as.numeric(coef(xbetw)))
+coef(xbetw)
+attr(coef(xbetw), names)
+attr(coef(xbetw), names, colnames(coef(xbetween)))
+attr(coef(xbetw), names, colnames(coef(xbetw)))
+help(attr)
+attr(coef(xbetw), "names")
+attr(coef(xbetw), "names")<-colnames(coef(xbetw))
+names(as.numeric(coef(xbetw))) <- colnames(coef(xbetw))
+as.numeric(coef(xbetw))
+attr(as.numeric(coef(xbetw)), "names")<-colnames(coef(xbetw))
+colnames(coef(xbetw))
+rownames(coef(xbetw))
+attr(as.numeric(coef(xbetw)), "names")<-rownames(coef(xbetw))
+coef <- as.numeric(coef(xbetw))
+coef
+names(coef)<-rownames(coef(xbetw))
+library(splm)#
+data(Produc, package = "Ecdat")#
+data(usaww)#
+Produc <- Produc[Produc$year<1975, ]#
+fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= FALSE, spatial.error=TRUE, effects = "fixed", control = list())#
+summary(sarargmfe)
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= FALSE, spatial.error=TRUE, effects = "random", control = list())#
+summary(sarargmre)
+sphtest(sarargmfe, sarargmre)
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "fixed", control = list())#
+summary(sarargmfe)#
+#
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=TRUE, effects = "random", control = list())#
+summary(sarargmre)
+sarargmfe<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=FALSE, effects = "fixed", control = list())#
+summary(sarargmfe)#
+#
+sarargmre<-spgm(fm, data=Produc,  listw = mat2listw(usaww),  lag= TRUE, spatial.error=FALSE, effects = "random", control = list())#
+summary(sarargmre)
+sphtest(sarargmfe, sarargmre)
+help(tapply)
+spgm
+library(splm)
+spgm
+Ws<-Matrix(,7,7)
+Ws
+Ws<-Matrix(1,7,7)
+Ws
+Diagonal(Ws)
+help(Diagonal)
+Ws<-Matrix(0,7,7)
+Diagonal(Ws)<-1
+diagonal(Ws)<-1
+diag(Ws)<-1
+Ws
+diag(Ws)
+Diagonal(30)
+library(splm)#
+data(Produc, package = "Ecdat")#
+data(usaww)#
+Produc <- Produc[Produc$year<1975, ]#
+fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
+	ml <- spfeml(formula, data=data, index=index, listw, model="error", effects="pooled")
+	ml <- spfeml(fm, data=data, index = NULL, listw= mat2listw(usaww), model="error", effects="pooled")
+	ml <- spfeml(fm, data=Produc, index = NULL, listw= mat2listw(usaww), model="error", effects="pooled")
+summary(ml)
+data = Produc
+  index <- data[,1]#
+  tindex <- data[,2]#
+  X<-model.matrix(formula,data=data)
+formula = fm
+  index <- data[,1]#
+  tindex <- data[,2]#
+  X<-model.matrix(formula,data=data)
+  y<-model.response(model.frame(formula,data=data))#
+  names(index)<-row.names(data)#
+  ind<-index[which(names(index)%in%row.names(X))]#
+  tind<-tindex[which(names(index)%in%row.names(X))]#
+  oo<-order(tind,ind)#
+  X<-X[oo,]#
+  y<-y[oo]#
+  ind<-ind[oo]#
+  tind<-tind[oo]#
+  N<-length(unique(ind))#
+  k<-dim(X)[[2]]#
+  T<-max(tapply(X[,1],ind,length))#
+  NT<-length(ind)
+	indic<-seq(1,T)#
+	inde<-as.numeric(rep(indic,each=N)) #
+	ind1<-seq(1,N)#
+	inde1<-as.numeric(rep(ind1,T))#
+#
+	lambda<-ml$spat.coef#
+	eML<-residuals(ml)#
+#
+ 	Ws<-listw2dgCMatrix(listw) #
+	Wst<-t(Ws)  #
+	B<- Diagonal(N)-lambda*Ws
+listw<- mat2listw(usaww)
+	indic<-seq(1,T)#
+	inde<-as.numeric(rep(indic,each=N)) #
+	ind1<-seq(1,N)#
+	inde1<-as.numeric(rep(ind1,T))#
+#
+	lambda<-ml$spat.coef#
+	eML<-residuals(ml)#
+#
+ 	Ws<-listw2dgCMatrix(listw) #
+	Wst<-t(Ws)  #
+	B<- Diagonal(N)-lambda*Ws
+	BpB<-crossprod(B)#
+	BpB2 <- BpB %*% BpB#
+	BpBi<- solve(BpB)#
+tr<-function(R) sum(diag(R))#
+	trBpB<-tr(BpB)
+	eme<-tapply(eML,inde1,mean)#
+	emme<-eML - rep(eme,T)#
+	sigmav2<-crossprod(eML,emme)/(N*(T-1))	sigmav4<-sigmav2^2
+	sigmav2<-crossprod(eML,emme)/(N*(T-1)#
+	sigmav4<-sigmav2^2
+	sigmav2<-crossprod(eML,emme)/(N*(T-1))
+	sigmav4<-sigmav2^2
+	sigmav2
+	sigmav4
+yybis<-function(q){ #
+	wq<-rep(q,T)#
+	tmp<-wq%*%eML#
+					}
+	BBu<-apply(BpB2,1,yybis)#
+	BBu<-rep(BBu,T)
+	upBBu<-crossprod(eML,BBu)
+	BBu
+	upBBu<-crossprod(eML,BBu)
+	upBBu
+	Dmu<- -(T/(2*sigmav2))*trBpB + (1/(2*sigmav4))*upBBu
+	Dmu
+	WpB<-Wst%*%B#
+	BpW<-crossprod(B, Ws)#
+	WpBplBpW <-WpB + BpW#
+	G<-WpBplBpW %*% BpBi
+	g<-tr(G)
+	g
+	WpB<-Wst%*%B#
+	BpW<-crossprod(B, Ws)#
+	WpBplBpW <-WpB + BpW#
+	bigG<-WpBplBpW %*% BpBi
+	smalle<-tr(BpB2)#
+	smalld<-tr(WpBplBpW)#
+	smallh<-trBpB#
+	smallg<-tr(bigG)#
+	smallc<-tr(bigG%*%bigG)
+	NUM<- ((2*sigmav4)/T)*((N*sigmav4*smallc)-(sigmav4*smallg^2))   ###equation 2.30 in the paper
+	NUM
+	NUM<- ((2 * sigmav4)/T) * ((N*sigmav4*smallc)-(sigmav4*smallg^2))   ###equation 2.30 in the paper
+	NUM
+	DENft<- NT*sigmav4* smalle * smallc
+	DENft
+	DENst<- N*sigmav4* smalld^2
+	DENtt<- T*sigmav4* smallg^2 * smalle
+DENtt
+DENst
+DENfot<- 2* sigmav4 *smallg * smallh* smalld
+DENfot
+	DENfit<- sigmav4 * smallh^2* smallc
+DENfit
+	DEN<- DENft - DENst - DENtt + DENfot - DENfit
+DEN
+	LMmu <- Dmu^2*NUM / DEN
+LMmu
+	LMmustar<- sqrt(LMmu)
+LMmustar


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

Added: pkg/R/REmod.R
===================================================================
--- pkg/R/REmod.R	                        (rev 0)
+++ pkg/R/REmod.R	2011-04-20 21:39:04 UTC (rev 102)
@@ -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/bsktest.R
===================================================================
--- pkg/R/bsktest.R	                        (rev 0)
+++ pkg/R/bsktest.R	2011-04-20 21:39:04 UTC (rev 102)
@@ -0,0 +1,1044 @@
+`bsktest` <-
+function(x,...){
+  UseMethod("bsktest")
+}
+
+
+`bsktest.splm` <-
+function(x, listw, index=NULL, test=c("CLMlambda","CLMmu"), ...){
+	
+	switch(match.arg(test), CLMlambda = {
+
+    bsk = clmltest.model(x,listw, index, ...)
+
+  }, CLMmu = {
+
+    bsk = clmmtest.model(x,listw, index, ... )
+
+  })
+
+  return(bsk)
+}
+
+
+`bsktest.lm` <-
+function(x, listw, index=NULL, test=c("SLM1","SLM2","LMJOINT"), ...){
+	
+	switch(match.arg(test), SLM1 = {
+
+    bsk = slm1test.model(x,listw, index, ...)
+
+  }, SLM2 = {
+
+    bsk = slm2test.model(x,listw, index, ... )
+
+  }, LMJOINT = {
+
+    bsk = LMHtest.model(x,listw, index, ...)
+
+  })
+
+  return(bsk)
+}
+
+
+`bsktest.formula` <-
+function(x, data, listw, test=c("SLM1","SLM2","LMJOINT","CLMlambda","CLMmu"), index=NULL, ...){
+  
+
+switch(match.arg(test), SLM1 = {
+
+    bsk = slm1test(x, data, index,  listw, ...)
+
+  }, SLM2 = {
+
+    bsk = slm2test(x, data, index,  listw, ...)
+
+  }, LMJOINT = {
+
+    bsk = LMHtest(x, data, index,  listw, ...)
+
+  }, CLMlambda = {
+
+    bsk = clmltest(x, data, index,  listw, ...)
+
+  }, CLMmu = {
+
+    bsk = clmmtest(x, data, index,  listw, ...)
+
+  })
+
+  return(bsk)
+
+}
+
+
+
+
+
+`slm1test.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 listw 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)
+
+	XpXi<-solve(crossprod(x))
+
+   cl<-match.call()
+  oo<-order(tind,ind)
+  x<-x[oo,]
+  y<-y[oo]
+  e<-e[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  N<-length(unique(ind))
+  k<-dim(x)[[2]]
+  T<-max(tapply(x[,1],ind,length))
+  NT<-length(ind)
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N))
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) 
+
+
+
+		JIe<-tapply(e,inde1,sum)
+		JIe<-rep(JIe,T) 
+		G<-(crossprod(e,JIe)/ee)-1 
+tr<-function(R) sum(diag(R))
+
+		LM1<-sqrt((NT/(2*(T-1))))*as.numeric(G) 
+		
+		s<-NT-k 
+		B<-XpXi%*%t(x)   
+		
+fun<-function(Q) tapply(Q,inde1,sum) 
+		JIx<-apply(x,2,fun)
+		JIX<-matrix(,NT,k)
+for (i in 1:k) JIX[,i]<-rep(JIx[,i],T) ## "NOTE ON THE TRACE.R"
+		di<-numeric(NT)
+		XpJIX<-crossprod(x,JIX)
+		d1<-NT-tr(XpJIX%*%XpXi) 
+		Ed1<-d1/s 
+
+		di2<-numeric(NT)
+		JIJIx<-apply(JIX,2,fun)
+		JIJIX<-matrix(,NT,k)
+for (i in 1:k) JIJIX[,i]<-rep(JIJIx[,i],T)
+		JIJIxxpx<-JIJIX%*%XpXi
+		di1<- crossprod(x, JIJIxxpx)
+		tr1<-tr(di1)
+		XpIJX<-crossprod(x,JIX)
+		fp<-XpIJX%*%B
+		sp<-JIX%*%XpXi
+		tr3<-tr(fp%*%sp)
+		fintr<-NT*T-2*tr1+tr3 
+		Vd1<-2*(s*fintr - (d1^2))/s^2*(s+2) 
+
+SLM1<-((G+1)- Ed1)/sqrt(Vd1) 
+
+STAT2<- qnorm(0.95,lower.tail=TRUE)
+	statistics<-SLM1
+  pval <- pnorm(SLM1, lower.tail=FALSE)
+	
+  names(statistics)="SLM1"
+	method<- "Baltagi, Song and Koh SLM1 marginal test"
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(formula), alternative="Random effects")
+  class(RVAL) <- "htest"
+  return(RVAL)
+}
+
+
+
+`slm2test.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]
+
+  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)
+  XpXi<-solve(crossprod(x))
+
+   cl<-match.call()
+  oo<-order(tind,ind)
+  x<-x[oo,]
+  y<-y[oo]
+  e<-e[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  N<-length(unique(ind))
+  k<-dim(x)[[2]]
+  T<-max(tapply(x[,1],ind,length))
+  NT<-length(ind)
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) 
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) 
+
+
+
+		Ws<-listw2dgCMatrix(listw) 
+		Wst<-t(Wst)  
+		WWp<-(Ws+Wst)/2 
+
+yy<-function(q){ 
+	wq<-WWp%*%q
+	wq<-as.matrix(wq)
+	}
+
+		IWWpe<-unlist(tapply(e,inde,yy)) 
+		H<-crossprod(e,IWWpe)/crossprod(e) 
+		W2<-Ws%*%Ws 
+		WW<-crossprod(Ws) 
+    tr<-function(R) sum(diag(R))
+	b<-tr(W2+WW) 
+		LM2<-sqrt((N^2*T)/b)*as.numeric(H)
+		s<-NT-k
+lag<-function(QQ)lag.listw(listw,QQ)
+fun2<-function(Q) unlist(tapply(Q,inde,lag))
+	Wx<-apply(x,2,fun2)
+	WX<-matrix(Wx,NT,k)
+	XpWx<-crossprod(x,WX)
+	D2M<-XpWx%*%XpXi 
+	Ed2<- (T*sum(diag(Ws)) - tr(D2M))/s 
+
+	WWx<-apply(WX,2,fun2)
+	WWX<-matrix(WWx,NT,k)
+	XpWWX<-crossprod(x,WWX)				
+	spb<-XpWWX%*%XpXi
+	spbb<-tr(spb)
+	tpb<-XpWx%*%XpXi%*%XpWx%*%XpXi
+	fintr2<-T*tr(W2) - 2* spbb + tr(tpb)
+	Vd2<-2*(s*fintr2 - (sum(diag(D2M))^2))/s^2*(s+2) 
+	We<-unlist(tapply(e,inde,function(W) lag.listw(listw,W)))
+	d2<-crossprod(e,We)/ee
+	
+	SLM2<- (d2-Ed2)/sqrt(Vd2) 
+
+STAT2<- qnorm(0.95,lower.tail=TRUE)
+	statistics<-SLM2
+  pval <- pnorm(SLM2, lower.tail=FALSE)
+
+  names(statistics)="SLM2"
+	method<- "Baltagi, Song and Koh SLM2 marginal test"
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(formula), alternative="Spatial autocorrelation")
+  class(RVAL) <- "htest"
+  return(RVAL)
+
+
+}
+
+
+
+`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]
+
+  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)
+	XpXi<-solve(crossprod(x))
+
+   cl<-match.call()
+  oo<-order(tind,ind)
+  x<-x[oo,]
+  y<-y[oo]
+  e<-e[oo]
+  ind<-ind[oo]
+  tind<-tind[oo]
+
+  N<-length(unique(ind))
+  k<-dim(x)[[2]]
+  T<-max(tapply(x[,1],ind,length))
+  NT<-length(ind)
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) 
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) 
+
+
+
+		JIe<-tapply(e,inde1,sum)
+		JIe<-rep(JIe,T) 
+		G<-(crossprod(e,JIe)/ee)-1 
+      tr<-function(R) sum(diag(R))
+
+		LM1<-sqrt((NT/(2*(T-1))))*as.numeric(G) 
+
+
+####calculate the elements of LMj, LM1, SLM1
+		Ws<-listw2dgCMatrix(listw) 
+		Wst<-t(Wst)  
+		WWp<-(Ws+Wst)/2 
+
+yy<-function(q){
+	wq<-WWp%*%q
+	wq<-as.matrix(wq)
+	}
+	
+		IWWpe<-unlist(tapply(e,inde,yy)) 
+		H<-crossprod(e,IWWpe)/crossprod(e) 
+		W2<-Ws%*%Ws 
+		WW<-crossprod(Ws) 
+      tr<-function(R) sum(diag(R))
+	   b<-tr(W2+WW) 
+		LM2<-sqrt((N^2*T)/b)*as.numeric(H)
+
+
+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"
+
+  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)
+
+
+}
+
+
+
+
+`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
+
+
+
+  oo<-order(tind,ind)
+  X<-X[oo,]
+  y<-y[oo]
+  N<-length(unique(ind))
+  k<-dim(X)[[2]]
+  T<-max(tapply(X[,1],ind,length))
+  NT<-length(ind)
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) 
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) 
+
+
+
+	eme<-tapply(eML,inde1,mean)
+	emme<-eML - rep(eme,T)
+	sigmav<-crossprod(eML,emme)/(N*(T-1)) 
+	sigma1<-crossprod(eML,rep(eme,T))/N 
+	
+	c1<-sigmav/sigma1^2
+	c2<-1/sigmav
+	c1e<-as.numeric(c1)*eML
+	Ws<-listw2dgCMatrix(listw) 
+	Wst<-t(Ws) 
+	WpsW<-Wst+Ws
+
+yybis<-function(q){
+	wq<-(WpsW)%*%q
+	wq<-as.matrix(wq)
+	}
+
+	Wc1e<-unlist(tapply(eML,inde,yybis)) 
+	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) 
+
+	W2<-Ws%*%Ws 
+	WW<-crossprod(Ws) 
+tr<-function(R) sum(diag(R))
+	b<-tr(W2+WW) 
+	LMl1<-D^2 / (((T-1)+as.numeric(sigmav)^2/as.numeric(sigma1)^2)*b) 
+
+	LMlstar<-sqrt(LMl1) 
+	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)"
+
+  dname <- deparse(formula)
+  RVAL <- list(statistic = statistics,
+               method = method,
+               p.value = pval, data.name=deparse(formula), alternative="Spatial autocorrelation")
+  class(RVAL) <- "htest"
+  return(RVAL)
+
+}
+
+
+
+
+`clmmtest.model` <-
+function(x, listw, index, ...){
+
+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
+
+  oo<-order(tind,ind)
+  X<-X[oo,]
+  y<-y[oo]
+
+  N<-length(unique(ind))
+  k<-dim(X)[[2]]
+  T<-max(tapply(X[,1],ind,length))
+  NT<-length(ind)
+	indic<-seq(1,T)
+	inde<-as.numeric(rep(indic,each=N)) 
+	ind1<-seq(1,N)
+	inde1<-as.numeric(rep(ind1,T)) 
+
+	lambda<-x$spat.coef
+
+ 	Ws<-listw2dgCMatrix(listw) 
+	Wst<-t(Ws)  
+	B<- Diagonal(N)-lambda*Ws
+
+	
+	BpB<-crossprod(B)
+	BpB2 <- BpB %*% BpB
+	BpBi<- solve(BpB)
+tr<-function(R) sum(diag(R))
+	trBpB<-tr(BpB)
+
+	eme<-tapply(eML,inde1,mean)
+	emme<-eML - rep(eme,T)
+	sigmav2<-crossprod(eML,emme)/(N*(T-1))
+	sigmav4<-sigmav2^2
+
+
+yybis<-function(q){ 
+	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<-crossprod(B, Ws)
+	WpBplBpW <-WpB + BpW
+	bigG<-WpBplBpW %*% BpBi
+	
+	smalle<-tr(BpB2)
+	smalld<-tr(WpBplBpW)
+	smallh<-trBpB
+	smallg<-tr(bigG)
+	smallc<-tr(bigG%*%bigG)
+	
+	NUM<- ((2 * sigmav4)/T) * ((N*sigmav4*smallc)-(sigmav4*smallg^2))   ###equation 2.30 in the paper
+	DENft<- NT*sigmav4* smalle * smallc
+	DENst<- N*sigmav4* smalld^2
+	DENtt<- T*sigmav4* smallg^2 * smalle
+	DENfot<- 2* sigmav4 *smallg * smallh* smalld
+	DENfit<- sigmav4 * smallh^2* smallc
+	DEN<- DENft - DENst - DENtt + DENfot - DENfit
+	LMmu <- Dmu^2*NUM / DEN
+	
+	LMmustar<- sqrt(LMmu)
+  
+  statistics<-LMmustar
[TRUNCATED]

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


More information about the Splm-commits mailing list