[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