[Splm-commits] r148 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 21 15:49:23 CET 2013
Author: gpiras
Date: 2013-03-21 15:49:23 +0100 (Thu, 21 Mar 2013)
New Revision: 148
Modified:
pkg/R/.Rapp.history
pkg/R/likelihoodsFE.R
pkg/R/spfeml.R
Log:
updated spfeml
Modified: pkg/R/.Rapp.history
===================================================================
--- pkg/R/.Rapp.history 2013-02-27 17:13:06 UTC (rev 147)
+++ pkg/R/.Rapp.history 2013-03-21 14:49:23 UTC (rev 148)
@@ -1,371 +1,260 @@
-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")
+install.packages("spatialprobit")
+library(spatialprobit)
+help(spatialprobit)
+install.packages("tmvtnorm")
+library(tmvtnorm)
+install.packages("mvtnorm")
+library(mvtnorm)
+library(Matrix)
+I_n_dense <- diag(10000)
+print(object.size(I_n_dense), units = "Mb")
+rm(I_n_dense)
+I_n_dense <- Diagonal(10000)
+print(object.size(I_n_dense), units = "Mb")
+I_n_dense
+rm(I_n_dense)
+I_n_dense <- sparseMatrix(i=1:10000, j =1:10000, x=1)
+I_n_dense
+print(object.size(I_n_dense), units = "Mb")
+rm(I_n_dense)
+I_n_dense <- Diagonal(10000)
+print(object.size(I_n_dense), units = "b")
+I_n_dense <- sparseMatrix(i=1:10000, j =1:10000, x=1)
+print(object.size(I_n_dense), units = "b")
+library(Matrix)
+I_n_dense <- diag(10000)
+print(object.size(I_n_dense), units = "b")
+print(object.size(I_n_dense), units = "Mb")
+rm(I_n_dense)
+I_n_dense <- sparseMatrix(i=1:10000, j =1:10000, x=1)
+print(object.size(I_n_dense), units = "Mb")
+print(object.size(I_n_dense), units = "b")
+rm(I_n_dense)
+I_n_dense <- Diagonal(10000)
+print(object.size(I_n_dense), units = "Mb")
+print(object.size(I_n_dense), units = "b")
+library(spatialprobit)
+n<- 400
+beta <- c(0, 1, -1)
+rho <- 0.75
+X <- cbind(intercept = 1, x = rnorm(n), y = rnorm(n))
+I_n <- sparseMatrix(i = 1:n, j = 1:n, x = 1)
+kNearestNeighbors
+W <- kNearestNeighbors(x = rnorm(n), y = rnorm(n), k = 6)
+W
+kNearestNeighbors
+knearneigh
+knearneighbours
+W <- kNearestNeighbors(x = 1:n, y =1:n, k = 6)
+help(kNearestNeighbors)
+library(spdep)
+eps <- rnorm(n = n, mean = 0, sd = 1)
+eps
+z <- solve(qr(I_n - rho * W), X %*% beta + eps)
+qr(I_n - rho * W)
+y <- as.vector(z >= 0)
+y
+sarprobit.fit1 <- sar_probit_mcmc(y, X, W, ndraw = 1000, burn.in = 200, thinning = 1, m = 10)
+summary(sarprobit.fit1)
+plot(sarprobit.fit1)
+library(igraph)
+install.packages("igraph")
+library(igraph)
+set.seed(1.2345)
+n <- 200
+branch <- 3
+probability <- branch/n
+probability
+grandom <- igraph::erdos.renyi.game(n = n, p.or.m = probability, type = "gnp", directed = F, loops = F)
+grandom
+V(grandom)$name <- 1:n
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+A <- igraph::get.adjacency(grandom, type = "both", sparse = T)
+A
+W <- A/rowSums(A)
+plot(grandom, vertex.label.family = "sans", vertex.size = 2, vertex.label = "", layout = layout.fruchterman.reingold)
+x <- rnorm(n)
+X <- cbind(intercept = rep(1, n), x = x)
+p <- 0.3
+beta <- c(-1, 2)
+I_n <- sparseMatrix(i = 1:n, j = 1:n, x = 1)
+z <- solve(qr(I_n - p * W), X %*% beta + rnorm(n))
+y <- as.real(z >= 0)
+sarprobit.fit <- sar_probit_mcmc(y, X, W, ndraw = 3000, burn.in = 200, thinning = 1)
+y
+z <- solve(qr(I_n - p * W), X %*% beta + rnorm(n))
+z
+I_n <- sparseMatrix(i = 1:n, j = 1:n, x = 1)
+I_n
+x <- rnorm(n)#
+X <- cbind(intercept = rep(1, n), x = x)
+X
+p <- 0.3#
+ beta <- c(-1, 2)#
+ I_n <- sparseMatrix(i = 1:n, j = 1:n, x = 1)#
+z <- solve(qr(I_n - p * W), X %*% beta + rnorm(n))
+z
+qr(I_n - p * W)
+beta
+qr(I_n - p * W), X %*% beta
+solve(qr(I_n - p * W), X %*% beta + rnorm(n))
+x <- rnorm(n)
+x
+X <- cbind(intercept = rep(1, n), x = x)
+X
+p <- 0.3
+beta <- c(-1, 2)
+beta
+p <- 0.3
+I_n <- sparseMatrix(i = 1:n, j = 1:n, x = 1)
+solve(qr(I_n - p * W), X %*% beta + rnorm(n))
+W <- A/rowSums(A)
+W
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+grandom
+V(grandom)$name <- 1:n
+igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+help(igraph::get.adjacency)
+W <- A/rowSums(A)
+W
+grandom <- igraph::erdos.renyi.game(n = n, p.or.m = probability, type = "gnp", directed = F, loops = F)
+grandom
+V(grandom)$name <- 1:n
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+load("CMK.RData")
+library(igraph)#
+set.seed(1.2345)#
+n <- 200#
+branch <- 3#
+probability <- branch/n#
+grandom <- igraph::erdos.renyi.game(n = n, p.or.m = probability, type = "gnp", directed = F, loops = F)#
+V(grandom)$name <- 1:n
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+W <- A/rowSums(A)
+A <- igraph::get.adjacency(grandom, type = "both", binary = T, sparse = T)
+W <- A/rowSums(A)
+t = 2
+T=2
+IT<- diag(T)
+IT
+N=3
+IN <- diag(N)
+IN
+JT<- matrix(1,2,2)
+JT
+Jbar<- 1/T*JT
+Jbar
+Q=IT-Jbar;
+Q
+Q=IT-Jbar
+Q
+eigen(Q)
+vec<-eigen(Q)
+vec
+vec$vectors[vec$vectors==1]
+vec$vectors[vec$vectors==1,]
+vec$vectors[vec$values==1]
+vec$vectors[,vec$values==1]
+kroneker(vec$vectors[,vec$values==1] , IN)
+kronecker(vec$vectors[,vec$values==1] , IN)
+kronecker(t(vec$vectors[,vec$values==1] ), IN)
+vec$vectors[vec$values==1]
+vec$vectors[vec$values==1]
+vec$vectors[vec$vectors==1,]
+vec$vectors[vec$values==1,]
+kronecker(t(vec$vectors[,vec$values==1] ), IN)
+Fmat <- vec$vectors[,vec$values==1]
+Fmat
+kronecker(t(Fmat), IN)
+Ftm <- kronecker(t(Fmat), IN)
+Ftm
+IT <- Diagonal(T)
+IT <- Diagonal(T)
+IT <- Diag(T)
+library(Matrix)
+IT <- Diag(T)
+IT <- Diag(T)
+IT <- Diagonal(T)
+IT
+IN <- diag(n)
+IN <- Diagonal(n)
+n=N
+IN <- Diagonal(n)
+IN
+JT <- matrix(1,T,T)
+JT
+Jbar <- 1/T * JT
+Jbar
+Qmat <-IT - Jbar
+Qmat
+vec <- eigen(Qmat)
+vec
+Fmat <- vec$vectors[,vec$values==1]
+Fmat
+Ftm <- kronecker(t(Fmat), IN)
+Ftm
+Fmat <- vec$vectors[,vec$values==L1]
+Fmat <- vec$vectors[,vec$values==1L]
+Fmat
+Ftm <- kronecker(t(Fmat), IN)
+Ftm
+iotan <- matrix(1,1,n)
+iotan
+iotan <- matrix(1,n,1)
+iotan
+Jnbar <-1/n * iotan %*% t(iotan)
+Jnbar
+Qmat1 <- IN - Jnbar
+Qmat1
+vec1 <- eigen(Qmat1)
+vec1
+Fmat1 <- vec1$vectors[,vec1$values==1L]
+Fmat1
+FFmat<- kronecker(t(Fmat), t(Fmat1))
+FFmat
+FFmat<- kronecker(Fmat, Fmat1)
+FFmat
+Fmat
+Fmat1
+Fmat1
+Fmat
+Fmat <- matrix(vec$vectors[,vec$values==1L], T, T-1)
+Fmat
+Ftm <- kronecker(t(Fmat), IN)
+Ftm
+Fmat1 <- matrix(vec1$vectors[,vec1$values==1L], n, n-1)
+Fmat1
+FFmat<- kronecker(Fmat, Fmat1)
+FFmat
+FFmat<- kronecker(t(Fmat), Fmat1)
+FFmat
+Ftm
+FFmat<- kronecker(t(Fmat), Fmat1)
+FFmat
+FFmat<- kronecker(Fmat, Fmat1)
+FFmat
+FFmat<- kronecker(Fmat, t(Fmat1))
+FFmat
+FFmat<- kronecker(t(Fmat), t(Fmat1))
+FFmat
+help(spfeml)
+help(spml)
+help(spfeml)
+help(spml)
+help(splm)
+help(spml)
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
+help(spml)
+help(spgm)
+library(spdep)
+library(spdep)
+help(sacsarlm)
+sacsarlm
+help(sarlm)
+lagsarlm
+help(lagsarlm)
+errorsarlm
Modified: pkg/R/likelihoodsFE.R
===================================================================
--- pkg/R/likelihoodsFE.R 2013-02-27 17:13:06 UTC (rev 147)
+++ pkg/R/likelihoodsFE.R 2013-03-21 14:49:23 UTC (rev 148)
@@ -1,3 +1,4 @@
+#### SAR PANEL
`conclikpan` <- function(lambda, env){
e0e0 <- get("e0e0", envir = env)
@@ -10,19 +11,159 @@
sigma2 <- Nsig/NT
ldet <- do_ldet(lambda, env)
- #ret <- T * ldet - ( (NT)/2 )*log(2*pi*sigma2) - (1/(2*sigma2))*Nsig
+
ret <- - (NT/2)*log(Nsig) + T * ldet
-# ret <- (NT/2)*log(Nsig) - T * ldet
-#print(T*log(det))
+
+
if (get("verbose", envir=env))
- cat("rho:\t", lambda, "\tfunction value:\t", ret,
+ cat("lambda:\t", lambda, "\tfunction value:\t", ret,
"\n")
ret
}
+
+splaglm<-function(env, zero.policy = zero.policy, interval = interval){
+
+xt <- get("xt", envir = env)
+yt <- get("yt", envir = env)
+wyt <- get("wyt", envir = env)
+con<-get("con", envir = env)
+NT<-get("NT", envir = env)
+T<-get("T", envir = env)
+listw<-get("listw", envir = env)
+inde<-get("inde", envir = env)
+interval1 <- get("interval1", envir = env)
+
+ XpX<-crossprod(xt)
+ b0<-solve(XpX,crossprod(xt,yt)) ####y on X
+ b1<-solve(XpX,crossprod(xt,wyt)) ####Wy on x
+ e0<-yt - xt%*% b0
+ e1<-wyt - xt%*% b1
+ e0e0<-crossprod(e0)
+ e1e1<-crossprod(e1)
+ e0e1<-t(e1)%*%e0
+
+assign("e0e0", e0e0, envir = env)
+assign("e1e1", e1e1, envir = env)
+assign("e0e1", e0e1, envir = env)
+
+
+opt <- optimize(conclikpan, interval = interval1, maximum = TRUE, env = env, tol = con$tol.opt)
+
+
+#opt <- nlminb(0.02138744, conclikpan, lower = interval[1], upper= interval[2], env = env)
+
+ rho <- opt$maximum
+
+ if (isTRUE(all.equal(rho, interval[1])) || isTRUE(all.equal(rho,
+ interval[2])))
+ warning("rho on interval bound - results should not be used")
+
+ names(rho) <- "lambda"
+ LL <- opt$objective
+ optres <- opt
+
+ nm <- paste(method, "opt", sep = "_")
+ timings[[nm]] <- proc.time() - .ptime_start
+ .ptime_start <- proc.time()
+
+ lm.lag <- lm((yt - rho * wyt) ~ xt - 1)
+ p <- lm.lag$rank
+ r <- residuals(lm.lag)
+ fit <- yt - r
+ names(r) <- names(fit)
+ betas <- coefficients(lm.lag)
+ names(betas) <- colnames(xt)
+ coefs <- c(rho, betas)
+
+ SSE <- deviance(lm.lag)
+ s2 <- SSE/NT
+ coefsl <- c(s2, rho, betas)
+
+ timings[["coefs"]] <- proc.time() - .ptime_start
+ .ptime_start <- proc.time()
+ assign("first_time", TRUE, envir = env)
+
+
+### add numerical hessian
+if(!Hess){
+
+ fd <- fdHess(coefsl, f_sarpanel_hess, env)
+ mat <- fd$Hessian
+ fdHess<- solve(-(mat), tol.solve = tol.solve)
+ rownames(fdHess) <- colnames(fdHess) <- c("s2", "lambda", colnames(xt))
+
+ rho.se <- fdHess[2, 2]
+ sig.se <- fdHess[1, 1]
+ rest.se <- vcov(lm.target)
+
+
+}
+
+else{
+
+ tr <- function(A) sum(diag(A))
+ W <-listw2dgCMatrix(listw, zero.policy = zero.policy)
+ A <- solve(diag(NT/T) - rho * W)
+ WA <- W %*% A
+ one <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
+
+ lag<-function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(WA %*% TT), simplify=TRUE))
+ lag2<-function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(t(WA)%*%TT), simplify=TRUE))
+ WAxt<-apply(as.matrix(xt),2,lag)
+# print(WAxt)
+ WAWAxt<-apply(WAxt,2,lag2)
+ # print(WAWAxt)
+ xtWAWAxt <- crossprod(xt,WAWAxt)
+ xtWAxt <- crossprod(xt,WAxt)
+ xtxt <- crossprod(xt)
+ two <- 1/as.numeric(s2) * t(betas) %*% xtWAWAxt %*% betas
+ V <- one + two
+ zero <- rbind(rep(0, length(betas)))
+ col1 <- rbind(NT/(2 * (s2^2)), T*tr(WA)/s2, t(zero))
+ three <- (1/as.numeric(s2)) * xtWAxt %*% betas
+ col2 <- rbind(T*tr(WA)/s2, V, three )
+ col3 <- rbind(zero, t(three), 1/as.numeric(s2)* xtxt)
+ asyvar <- cbind(col1, col2, col3)
+ asyv <- solve(asyvar, tol = con$tol.solve)
+ rownames(asyv) <- colnames(asyv) <- c("sigma","lambda", colnames(xt))
+
+ rho.se <- sqrt(asyv[2, 2])
+ rest.se <- sqrt(diag(asyv))[-c(1:2)]
+ sig.se <- sqrt(asyv[1, 1])
+ asyvar1 <- asyv[-1,-1]
+
+ rownames(asyvar1) <- colnames(asyvar1) <- c("lambda", colnames(xt))
+
+}
+
+ return<-list(coeff=betas,rho=rho,s2=s2, rest.se=rest.se, rho.se=rho.se, sig.se = sig.se, asyvar1=asyvar1)
+}
+
+
+
+f_sarpanel_hess <- function (coefs, env)
+{
+
+ T<-get("T", envir = env)
+ NT<-get("NT", envir = env)
+ s2<- coefs[1]
+ lambda <- coefs[2]
+ beta <- coefs[-c(1,2)]
+ SSE <- s2*n
+ n <- NT/T
+ ldet <- do_ldet(lambda, env, which = 1)
+ ret <- (T * ldet - ((n*T/2) * log(2 * pi)) - (n*T/2) * log(s2) -
+ (1/(2 * s2)) * SSE)
+ if (get("verbose", envir = env))
+ cat("lambda:", lambda, " function:", ret,
+ " Jacobian:", ldet," SSE:", SSE, "\n")
+ ret
+}
+
+####### ERROR MODEL
sarpanelerror<-function (lambda, env=env)
{
-
yt<- get("yt", envir = env)
xt<- get("xt", envir = env)
wyt<- get("wyt", envir = env)
@@ -46,12 +187,127 @@
ret <- T*ldet - (NT/2) * log(SSE)
if (get("verbose", envir = env))
- cat("lambda:", lambda, " function:", ret, " Jacobian:", ldet, " SSE:", SSE, "\n")
+ cat("rho:", lambda, " function:", ret, " Jacobian:", ldet, " SSE:", SSE, "\n")
ret
}
+
+
+sperrorlm<-function(env, zero.policy = zero.policy, interval = interval){
+
+xt <- get("xt", envir = env)
+yt <- get("yt", envir = env)
+wyt <- get("wyt", envir = env)
+wxt<-get("wxt", envir = env)
+
+wy <- get("wy", envir = env)
+wx<-get("wx", envir = env)
+
+con<-get("con", envir = env)
+NT<-get("NT", envir = env)
+T<-get("T", envir = env)
+listw<-get("listw", envir = env)
+inde<-get("inde", envir = env)
+interval <- get("interval1", envir = env)
+
+opt <- optimize(sarpanelerror, interval = interval, maximum = TRUE, env = env, tol = con$tol.opt)
+
+
+#opt <- nlminb(0.5,sarpanelerror,lower = interval[1], upper= interval[2], env = env)
+#print(opt)
+
+ lambda <- opt$maximum
+ names(lambda) <- "rho"
+ LL <- opt$objective
+
+
+ lm.target <- lm(I(yt - lambda * wyt) ~ I(xt - lambda * wxt) -
+ 1)
+ r <- as.vector(residuals(lm.target))
+ p <- lm.target$rank
+ s2 <- crossprod(r)/NT
+ rest.se <- (summary(lm.target)$coefficients[, 2]) * sqrt((NT - p)/NT)
+ betas <- coefficients(lm.target)
+ names(betas) <- colnames(xt)
+ coefsl <- c(s2, lambda, betas)
+
+
+if(!Hess){
+
+ fd <- fdHess(coefsl, sarpanelerror_hess, env)
+ mat <- fd$Hessian
+ fdHess<- solve(-(mat), tol.solve = tol.solve)
+ rownames(fdHess) <- colnames(fdHess) <- c("s2", "rho",colnames(xt))
+
+ rho.se <- fdHess[2, 2]
+ s2.se <- fdHess[1, 1]
+
+}
+else{
+
+ tr <- function(A) sum(diag(A))
+ W <- listw2dgCMatrix(listw, zero.policy = zero.policy)
+ A <- solve(Diagonal(NT/T) - lambda * W)
+ WA <- W %*% A
+ asyvar <- matrix(0, nrow = 2 + p, ncol = 2 + p)
+ asyvar[1, 1] <- NT/(2 * (s2^2))
+ asyvar[2, 1] <- asyvar[1, 2] <- T*tr(WA)/s2
+ asyvar[2, 2] <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
+ asyvar[3:(p + 2), 3:(p + 2)] <- 1/as.numeric(s2) * (t(xt - lambda *wxt) %*% (xt - lambda * wxt))
+ asyv <- solve(asyvar, tol = con$tol.solve)
+ rownames(asyv) <- colnames(asyv) <- c("sigma","rho", colnames(xt))
+ s2.se <- sqrt(asyv[1, 1])
+ lambda.se <- sqrt(asyv[2, 2])
+ asyvar1 <- asyv[-1,-1]
+ rownames(asyvar1) <- colnames(asyvar1) <- c("rho", colnames(xt))
+
+}
+
+ return<-list(coeff=betas,lambda=lambda,s2=s2, rest.se=rest.se, lambda.se=lambda.se, s2.se = s2.se, asyvar1=asyvar1)
+}
+
+
+
+sarpanelerror_hess<-function (coef, env=env)
+{
+ yt<- get("yt", envir = env)
+ xt<- get("xt", envir = env)
+ wyt<- get("wyt", envir = env)
+ wxt<- get("wxt", envir = env)
+ wy<- get("wy", envir = env)
+ wx<- get("wx", envir = env)
+
+ listw<- get("listw", envir = env)
+ NT<- get("NT", envir = env)
+ inde<- get("inde", envir = env)
+ T<- get("T", envir = env)
+
+# coefsl <- c(s2, lambda, betas)
+
+ s2 <- coef[1]
+ lambda <- coef[2]
+ bb <- coef[-c(1,2)]
+
+ # yco <- yt - lambda * wyt
+ # xco <- xt - lambda * wxt
+ # bb<- solve(crossprod(xco),crossprod(xco, yco) )
+
+ # ehat<- yco - xco %*% bb
+ SSE <- s2 * NT
+ ldet <- do_ldet(lambda, env)
+
+ ret <- T*ldet - (NT/2) * log(SSE)
+
+if (get("verbose", envir = env))
+ cat("rho:", lambda, " function:", ret, " Jacobian:", ldet, " SSE:", SSE, "\n")
+ ret
+}
+
+
+###SARAR MODEL
+
sacsarpanel<-function (coefs, env)
{
T<-get("T", envir = env)
@@ -61,8 +317,8 @@
ldet1 <- do_ldet(coefs[1], env, which = 1)
ldet2 <- do_ldet(coefs[2], env, which = 2)
- ret <- (T * ldet1 + T * ldet2 - ((n*T/2) * log(2 * pi)) - (n*T/2) * log(s2) -
- (1/(2 * (s2))) * SSE)
+ ret <-(T * (ldet1 + ldet2)) - ((n*T/2) * log(2 * pi)) - (n*T/2) * log(s2)
+ - (1/(2 * (s2))) * SSE)
if (get("verbose", envir = env)) cat("rho:", coefs[1], " lambda:", coefs[2], " function:",
ret, " Jacobian1:", ldet1, " Jacobian2:", ldet2,
@@ -102,20 +358,27 @@
T<-get("T", envir = env)
listw<-get("listw", envir = env)
inde<-get("inde", envir = env)
+interval1 <- get("interval1", envir = env)
+interval2 <- get("interval2", envir = env)
pars <- con$pars
- lower <- con$lower
- upper <- con$upper
-if (!is.null(llprof)) {
+ lower <- c(interval1[1], interval2[1])
+ upper <- c(interval1[2], interval2[2])
+
+ if (!is.null(llprof)) {
llrho <- NULL
lllambda <- NULL
- if (length(llprof) == 1) {
+ if (length(llprof) == 1L) {
llrho <- seq(lower[1], upper[1], length.out = llprof)
lllambda <- seq(lower[2], upper[2], length.out = llprof)
llprof <- as.matrix(expand.grid(llrho, lllambda))
}
ll_prof <- numeric(nrow(llprof))
-for (i in 1:nrow(llprof)) ll_prof[i] <- sacsarpanel(llprof[i,], env = env)
+ for (i in 1:nrow(llprof)) ll_prof[i] <- sacsarpanel(llprof[i,
+ ], env = env)
+ nm <- paste(method, "profile", sep = "_")
+ timings[[nm]] <- proc.time() - .ptime_start
+ .ptime_start <- proc.time()
}
if (is.null(pars)) {
if (con$npars == 4L) {
@@ -124,25 +387,59 @@
mpars <- cbind(xseq, yseq)
}
else {
- xseq <- seq(lower[1], upper[1], (upper[1] - lower[1])/2) * 0.8
- yseq <- seq(lower[2], upper[2], (upper[2] - lower[2])/2) * 0.8
+ xseq <- seq(lower[1], upper[1], (upper[1] - lower[1])/2) *
+ 0.8
+ yseq <- seq(lower[2], upper[2], (upper[2] - lower[2])/2) *
+ 0.8
mpars <- as.matrix(expand.grid(xseq, yseq))
}
}
-
else {
mxs <- NULL
}
-
+ if (con$opt_method == "nlminb") {
if (is.null(pars)) {
- mxs <- apply(mpars, 1, function(pp) -nlminb(pp, sacsarpanel, env = env, control = con$opt_control, lower = lower, upper = upper)$objective)
+ mxs <- apply(mpars, 1, function(pp) -nlminb(pp, sacsarpanel,
+ env = env, control = con$opt_control, lower = lower,
+ upper = upper)$objective)
pars <- mpars[which.max(mxs), ]
- optres <- nlminb(pars, sacsarpanel, env = env, control = con$opt_control, lower = lower, upper = upper)
+ optres <- nlminb(pars, sacsarpanel, env = env, control = con$opt_control,
+ lower = lower, upper = upper)
}
else {
- optres <- nlminb(pars, sacsarpanel, env = env, control = con$opt_control, lower = lower, upper = upper)
+ optres <- nlminb(pars, sacsarpanel, env = env, control = con$opt_control,
+ lower = lower, upper = upper)
}
-
+ }
+ else if (con$opt_method == "L-BFGS-B") {
+ if (is.null(pars)) {
+ mxs <- apply(mpars, 1, function(pp) -optim(pars,
+ sacsarpanel, env = env, method = "L-BFGS-B", control = con$opt_control,
+ lower = lower, upper = upper)$objective)
+ pars <- mpars[which.max(mxs), ]
+ optres <- optim(pars, sacsarpanel, env = env, method = "L-BFGS-B",
+ control = con$opt_control, lower = lower, upper = upper)
+ }
+ else {
+ optres <- optim(pars, sacsarpanel, env = env, method = "L-BFGS-B",
+ control = con$opt_control, lower = lower, upper = upper)
+ }
+ }
+ else {
+ if (is.null(pars)) {
+ mxs <- apply(mpars, 1, function(pp) -optim(pars,
+ sacsarpanel, env = env, method = con$opt_method,
+ control = con$opt_control)$objective)
+ pars <- mpars[which.max(mxs), ]
+ optres <- optim(pars, sacsarpanel, env = env, method = con$opt_method,
+ control = con$opt_control)
+ }
+ else {
+ optres <- optim(pars, sacsarpanel, env = env, method = con$opt_method,
+ control = con$opt_control)
+ }
+ }
+
LL <- -optres$objective
if (optres$convergence != 0)
warning(paste("convergence failure:", optres$message))
@@ -151,8 +448,11 @@
names(rho) <- "rho"
lambda <- optres$par[1]
names(lambda) <- "lambda"
+ nm <- paste(method, "opt", sep = "_")
+ timings[[nm]] <- proc.time() - .ptime_start
+ .ptime_start <- proc.time()
+
-
lm.target <- lm(I(yt - lambda * wyt - rho * w2yt + rho * lambda *
w2wyt) ~ I(xt - rho * wxt) - 1)
@@ -163,16 +463,15 @@
s2 <- SSE/NT
betas <- coefficients(lm.target)
names(betas) <- colnames(xt)
+ coefs <- c(rho, lambda, betas)
- coefs <- c(rho, lambda, betas)
+###Add the vc matrix exact
fd <- fdHess(coefs, f_sacpanel_hess, env)
mat <- fd$Hessian
fdHess<- solve(-(mat), tol.solve = tol.solve)
rownames(fdHess) <- colnames(fdHess) <- c("lambda", "lambda",colnames(xt))
-
-
rho.se <- fdHess[2, 2]
lambda.se <- fdHess[1, 1]
rest.se <- vcov(lm.target)
@@ -202,155 +501,17 @@
ret
}
-sar_sac_hess_sse_panel <- function (rho, lambda, beta, env)
-{
- yl <- get("yt", envir = env) - lambda * get("wyt", envir = env) -
- rho * get("w2yt", envir = env) + rho * lambda * get("w2wyt",
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/splm -r 148
More information about the Splm-commits
mailing list