[Splm-commits] r157 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 30 23:22:06 CET 2013
Author: the_sculler
Date: 2013-03-30 23:22:05 +0100 (Sat, 30 Mar 2013)
New Revision: 157
Modified:
pkg/DESCRIPTION
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/spreml.R
pkg/man/sphtest.Rd
pkg/man/splm-package.Rd
pkg/man/spml.Rd
pkg/man/usaww.Rd
pkg/man/vcov.splm.Rd
Log:
Major revision ex CSDA paper for all functions depending on spreml.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-03-28 21:12:49 UTC (rev 156)
+++ pkg/DESCRIPTION 2013-03-30 22:22:05 UTC (rev 157)
@@ -1,10 +1,10 @@
Package: splm
Title: Econometric Models for Spatial Panel Data
-Version: 1.0-05
-Date: 2012-10-24
+Version: 1.1-00
+Date: 2013-03-30
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.15.0), MASS, nlme, spdep, plm, Matrix, bdsmatrix, spam, ibdreg, car, lmtest, Ecdat
+Depends: R (>= 2.15.0), MASS, nlme, spdep, plm, Matrix, bdsmatrix, spam, ibdreg, car, lmtest, Ecdat, maxLik
License: GPL-2
LazyLoad: yes
Modified: pkg/R/sarREmod.R
===================================================================
--- pkg/R/sarREmod.R 2013-03-28 21:12:49 UTC (rev 156)
+++ pkg/R/sarREmod.R 2013-03-30 22:22:05 UTC (rev 157)
@@ -38,17 +38,26 @@
## the fact that in this case the vcov matrix is block-diagonal
## calc. Wy (spatial lag of y)
+ ## (flexible fun accepting either listws or matrices for w)
Wy <- function(y, w, tind) { # lag-specific line
+ wyt <- function(y, w) { # lag-specific line
+ if("listw" %in% class(w)) { # lag-specific line
+ wyt <- lag.listw(w, y) # lag-specific line
+ } else { # lag-specific line
+ wyt <- w %*% y # lag-specific line
+ } # lag-specific line
+ return(wyt) # lag-specific line
+ } # lag-specific line
wy<-list() # lag-specific line
for (j in 1:length(unique(tind))) { # lag-specific line
yT<-y[tind==unique(tind)[j]] # lag-specific line
- wy[[j]] <- w %*% yT # lag-specific line
+ wy[[j]] <- wyt(yT, w) # lag-specific line
} # lag-specific line
return(unlist(wy)) # lag-specific line
} # lag-specific line
## lag y once for all
- wy <- Wy(y, w2, tind) # lag-specific line
+ wy <- Wy(y, w2, tind) # lag-specific line
## the sigma matrix is inverted during the GLS step and not before as
## in the other cases, to take advantage of specialized methods in the
@@ -116,7 +125,7 @@
e <- glsres[["ehat"]]
s2e <- glsres[["sigma2"]]
## calc ll
- zero <- t*log(detB(psi, w2)) # lag-specific line (else zero <- 0)
+ zero <- t*ldetB(psi, w2) # lag-specific line (else zero <- 0)
due <- detSigma(phi, n, t)
tre <- -(n * t)/2 * log(s2e)
cinque <- -1/(2 * s2e) * crossprod(e, solve(sigma, e))
@@ -158,7 +167,7 @@
## final parms
betas <- as.vector(beta[[1]])
- sigma2 <- as.numeric(beta[["sigma2"]])
+ sigma2 <- as.numeric(beta[["sigma2"]])
arcoef <- myparms[which(nam.errcomp=="lambda")] # lag-specific line
errcomp <- myparms[which(nam.errcomp!="lambda")]
names(betas) <- nam.beta
Modified: pkg/R/sarem2REmod.R
===================================================================
--- pkg/R/sarem2REmod.R 2013-03-28 21:12:49 UTC (rev 156)
+++ pkg/R/sarem2REmod.R 2013-03-30 22:22:05 UTC (rev 157)
@@ -1,10 +1,10 @@
sarem2REmod <-
-function (X, y, ind, tind, n, k, t, nT, w, w2, coef0 = rep(0, 3),
+function (X, y, ind, tind, n, k, t., nT, w, w2, coef0 = rep(0, 3),
hess = FALSE, trace = trace, x.tol = 1.5e-18, rel.tol = 1e-15,
- ...)
+ method="nlminb", ...)
{
- ## extensive function rewriting, Giovanni Millo 29/09/2010
+ ## extensive function rewriting, Giovanni Millo 27/03/2013
## structure:
## a) specific part
## - set names, bounds and initial values for parms
@@ -17,8 +17,7 @@
## - calc final covariances
## - make list of results
- # mark
- #print("uso versione 0") # done from saremsrREmod4.R
+ ## now using flex optimization and sparse matrix methods
## set names for final parms vectors
nam.beta <- dimnames(X)[[2]]
@@ -26,18 +25,13 @@
## initialize values for optimizer
myparms0 <- coef0
- ## set bounds for optimizer
- lower.bounds <- c(1e-08, -0.999, -0.999) # lag-specific line (3rd parm)
- upper.bounds <- c(1e+09, 0.999, 0.999) # lag-specific line (idem)
## modules for likelihood
- B <- function(lambda, w) diag(1, ncol(w)) - lambda * w
- detB <- function(lambda, w) det(B(lambda, w))
- invSigma <- function(philambda, n, t, w) {
- Jt <- matrix(1, ncol = t, nrow = t)
+ invSigma <- function(philambda, n, t., w) {
+ Jt <- matrix(1, ncol = t., nrow = t.)
#In <- diag(1, n)
- It <- diag(1, t)
- Jbart <- Jt/t
+ It <- diag(1, t.)
+ Jbart <- Jt/t.
Et <- It - Jbart
## retrieve parms
phi <- philambda[1]
@@ -45,29 +39,29 @@
## psi not used: here passing 4 parms, but works anyway
## because psi is last one
## calc inverse
- BB <- crossprod(B(lambda, w))
- invSigma <- kronecker( (1/(t*phi+1)*Jbart + Et), BB )
+ BB <- xprodB(lambda, w)
+ invSigma <- kronecker( (1/(t.*phi+1)*Jbart + Et), BB )
invSigma
}
- detSigma <- function(phi, lambda, n, t, w) {
- Jt <- matrix(1, ncol = t, nrow = t)
+ detSigma <- function(phi, lambda, n, t., w) {
+ Jt <- matrix(1, ncol = t., nrow = t.)
#In <- diag(1, n)
- It <- diag(1, t)
- Jbart <- Jt/t
+ It <- diag(1, t.)
+ Jbart <- Jt/t.
Et <- It - Jbart
- detSigma <- -n/2*log( det( (t*phi+1) * Jbart + Et) ) +
- t*log(detB(lambda, w))
+ detSigma <- -n/2*log( det( (t.*phi+1) * Jbart + Et) ) +
+ t.*ldetB(lambda, w)
detSigma
}
## likelihood function, both steps included
- ll.c <- function(philambda, y, X, n, t, w, w2, wy) {
+ ll.c <- function(philambda, y, X, n, t., w, w2, wy) {
## retrieve parms
phi <- philambda[1]
lambda <- philambda[2]
psi <- philambda[3] # lag-specific line
## calc inverse sigma
- sigma.1 <- invSigma(philambda, n, t, w)
+ sigma.1 <- invSigma(philambda, n, t., w)
## lag y
Ay <- y - psi * wy # lag-specific line
## do GLS step to get e, s2e
@@ -75,66 +69,130 @@
e <- glsres[["ehat"]]
s2e <- glsres[["sigma2"]]
## calc ll
- zero <- t*log(detB(psi, w2)) # lag-specific line (else zero <- 0)
- due <- detSigma(phi, lambda, n, t, w)
- tre <- -n * t/2 * log(s2e)
- quattro <- -1/(2 * s2e) * crossprod(e, sigma.1) %*% e
- const <- -(n * t)/2 * log(2 * pi)
+ zero <- t.*ldetB(psi, w2) # lag-specific line (else zero <- 0)
+ due <- detSigma(phi, lambda, n, t., w)
+ tre <- -n * t./2 * log(s2e)
+ quattro <- -1/(2 * s2e) * t(e) %*% sigma.1 %*% e
+ const <- -(n * t.)/2 * log(2 * pi)
ll.c <- const + zero + due + tre + quattro
## invert sign for minimization
llc <- -ll.c
}
+ ## set bounds for optimizer
+ lower.bounds <- c(1e-08, -0.999, -0.999) # lag-specific line (4th parm)
+ upper.bounds <- c(1e+09, 0.999, 0.999) # lag-specific line (idem)
+
+ ## constraints as cA %*% theta + cB >= 0
+ ## equivalent to: phi>=0, -1<=(rho, lambda, psi)<=1
+ ## NB in maxLik() optimization cannot start at the boundary of the
+ ## parameter space !
+ cA <- cbind(c(1, rep(0,4)),
+ c(0,1,-1,rep(0,2)),
+ c(rep(0,3), 1, -1))
+ cB <- c(0, rep(1,4))
## generic from here
## calc. Wy (spatial lag of y)
+ ## (flexible fun accepting either listws or matrices for w)
Wy <- function(y, w, tind) { # lag-specific line
+ wyt <- function(y, w) { # lag-specific line
+ if("listw" %in% class(w)) { # lag-specific line
+ wyt <- lag.listw(w, y) # lag-specific line
+ } else { # lag-specific line
+ wyt <- w %*% y # lag-specific line
+ } # lag-specific line
+ return(wyt) # lag-specific line
+ } # lag-specific line
wy<-list() # lag-specific line
for (j in 1:length(unique(tind))) { # lag-specific line
yT<-y[tind==unique(tind)[j]] # lag-specific line
- wy[[j]] <- w %*% yT # lag-specific line
+ wy[[j]] <- wyt(yT, w) # lag-specific line
} # lag-specific line
return(unlist(wy)) # lag-specific line
} # lag-specific line
## GLS step function
GLSstep <- function(X, y, sigma.1) {
- b.hat <- solve(crossprod(X, sigma.1) %*% X,
- crossprod(X, sigma.1) %*% y)
+ b.hat <- solve(t(X) %*% sigma.1 %*% X,
+ t(X) %*% sigma.1 %*% y)
ehat <- y - X %*% b.hat
- sigma2ehat <- (crossprod(ehat, sigma.1) %*% ehat)/(n * t)
+ sigma2ehat <- (t(ehat) %*% sigma.1 %*% ehat)/(n * t.)
return(list(betahat=b.hat, ehat=ehat, sigma2=sigma2ehat))
}
## lag y once for all
wy <- Wy(y, w2, tind) # lag-specific line
- ## 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,
+ ## optimization
+
+ ## adaptive scaling
+ parscale <- 1/max(myparms0, 0.1)
+
+ if(method=="nlminb") {
+
+ 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 = parscale,
+ control = list(x.tol = x.tol,
rel.tol = rel.tol, trace = trace),
- lower = lower.bounds, upper = upper.bounds)
+ lower = lower.bounds, upper = upper.bounds)
- ## log likelihood at optimum (notice inverted sign)
- myll <- -optimum$objective
- ## retrieve optimal parms
- myparms <- optimum$par
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms and H
+ myparms <- optimum$par
+ myHessian <- fdHess(myparms, function(x) -ll.c(x,
+ y, X, n, t., w, w2, wy))$Hessian # lag-specific line: wy
+ } else {
+
+ require(maxLik)
+
+ ## initial values are not allowed to be zero
+ maxout<-function(x,a) ifelse(x>a, x, a)
+ myparms0 <- maxout(myparms0, 0.01)
+
+ ## invert sign for MAXimization
+ ll.c2 <- function(phirholambda, y, X, n, t., w, w2, wy) {
+ -ll.c(phirholambda, y, X, n, t., w, w2, wy)
+ }
+
+ ## max likelihood
+ optimum <- maxLik(logLik = ll.c2,
+ grad = NULL, hess = NULL, start=myparms0,
+ method = method,
+ parscale = parscale,
+ constraints=list(ineqA=cA, ineqB=cB),
+ y = y, X = X, n = n, t. = t., w = w, w2 = w2, wy = wy)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- optimum$maximum # this one MAXimizes
+ ## retrieve optimal parms and H
+ myparms <- optimum$estimate
+ myHessian <- optimum$hessian
+ }
+
+
## one last GLS step at optimal vcov parms
- sigma.1 <- invSigma(myparms, n, t, w)
+ sigma.1 <- invSigma(myparms, n, t., w)
Ay <- y - myparms[length(myparms)] * wy # lag-specific line
beta <- GLSstep(X, Ay, sigma.1)
## final vcov(beta)
covB <- as.numeric(beta[[3]]) *
- solve(crossprod(X, sigma.1) %*% X)
+ solve(t(X) %*% sigma.1 %*% 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) - 1
+ ## error handler here for singular Hessian cases
+ covTheta <- try(solve(-myHessian), silent=TRUE)
+ if(class(covTheta) == "try-error") {
+ covTheta <- matrix(NA, ncol=nvcovpms+1,
+ nrow=nvcovpms+1)
+ warning("Hessian matrix is not invertible")
+ }
covAR <- covTheta[nvcovpms+1, nvcovpms+1, drop=FALSE]
covPRL <- covTheta[1:nvcovpms, 1:nvcovpms, drop=FALSE]
Modified: pkg/R/saremREmod.R
===================================================================
--- pkg/R/saremREmod.R 2013-03-28 21:12:49 UTC (rev 156)
+++ pkg/R/saremREmod.R 2013-03-30 22:22:05 UTC (rev 157)
@@ -1,7 +1,7 @@
saremREmod <-
-function (X, y, ind, tind, n, k, t, nT, w, w2, coef0 = rep(0, 3),
+function (X, y, ind, tind, n, k, t., nT, w, w2, coef0 = rep(0, 3),
hess = FALSE, trace = trace, x.tol = 1.5e-18, rel.tol = 1e-15,
- ...)
+ method="nlminb", ...)
{
## extensive function rewriting, Giovanni Millo 29/09/2010
@@ -26,18 +26,16 @@
## initialize values for optimizer
myparms0 <- coef0
- ## set bounds for optimizer
- lower.bounds <- c(1e-08, -0.999, -0.999) # lag-specific line (3rd parm)
- upper.bounds <- c(1e+09, 0.999, 0.999) # lag-specific line (idem)
## modules for likelihood
- B <- function(lambda, w) diag(1, ncol(w)) - lambda * w
- detB <- function(lambda, w) det(B(lambda, w))
- invSigma <- function(philambda, n, t, w) {
- Jt <- matrix(1, ncol = t, nrow = t)
+ BB.1 <- function(lambda, w) {
+ solve(xprodB(lambda, listw=w))
+ }
+ invSigma <- function(philambda, n, t., w) {
+ Jt <- matrix(1, ncol = t., nrow = t.)
In <- diag(1, n)
- It <- diag(1, t)
- Jbart <- Jt/t
+ It <- diag(1, t.)
+ Jbart <- Jt/t.
Et <- It - Jbart
## retrieve parms
phi <- philambda[1]
@@ -45,27 +43,28 @@
## psi not used: here passing 4 parms, but works anyway
## because psi is last one
## calc inverse
- BB <- crossprod(B(lambda, w))
- invSigma <- kronecker(Jbart, solve(t * phi * In + solve(BB))) +
- kronecker(Et, BB)
+ BB <- xprodB(lambda, w)
+ invSigma <- kronecker(Jbart,
+ solve(t. * phi * In + BB.1(lambda, w))) +
+ kronecker(Et, BB)
invSigma
}
- detSigma <- function(phi, lambda, n, t, w) {
+ detSigma <- function(phi, lambda, n, t., w) {
In <- diag(1, n)
- detSigma <- -1/2 * log(det(t * phi * In +
- solve(crossprod(B(lambda, w))))) +
- (t - 1) * log(detB(lambda, w))
+ detSigma <- -1/2 * log(det(t. * phi * In +
+ BB.1(lambda, w))) +
+ (t. - 1) * ldetB(lambda, w)
detSigma
}
## likelihood function, both steps included
- ll.c <- function(philambda, y, X, n, t, w, w2, wy) {
+ ll.c <- function(philambda, y, X, n, t., w, w2, wy) {
## retrieve parms
phi <- philambda[1]
lambda <- philambda[2]
psi <- philambda[3] # lag-specific line
## calc inverse sigma
- sigma.1 <- invSigma(philambda, n, t, w)
+ sigma.1 <- invSigma(philambda, n, t., w2)
## lag y
Ay <- y - psi * wy # lag-specific line
## do GLS step to get e, s2e
@@ -73,24 +72,47 @@
e <- glsres[["ehat"]]
s2e <- glsres[["sigma2"]]
## calc ll
- zero <- t*log(detB(psi, w2)) # lag-specific line (else zero <- 0)
- due <- detSigma(phi, lambda, n, t, w)
- tre <- -n * t/2 * log(s2e)
+ zero <- t.*ldetB(psi, w) # lag-specific line (else zero <- 0)
+ due <- detSigma(phi, lambda, n, t., w2)
+ tre <- -n * t./2 * log(s2e)
quattro <- -1/(2 * s2e) * crossprod(e, sigma.1) %*% e
- const <- -(n * t)/2 * log(2 * pi)
+ const <- -(n * t.)/2 * log(2 * pi)
ll.c <- const + zero + due + tre + quattro
## invert sign for minimization
llc <- -ll.c
}
+ ## set bounds for optimizer
+ lower.bounds <- c(1e-08, -0.999, -0.999) # lag-specific line (4th parm)
+ upper.bounds <- c(1e+09, 0.999, 0.999) # lag-specific line (idem)
+
+ ## constraints as cA %*% theta + cB >= 0
+ ## equivalent to: phi>=0, -1<=(rho, lambda, psi)<=1
+ ## NB in maxLik() optimization cannot start at the boundary of the
+ ## parameter space !
+ cA <- cbind(c(1, rep(0,4)),
+ c(0,1,-1,rep(0,2)),
+ c(rep(0,3), 1, -1))
+ cB <- c(0, rep(1,4))
+
+
## generic from here
## calc. Wy (spatial lag of y)
+ ## (flexible fun accepting either listws or matrices for w)
Wy <- function(y, w, tind) { # lag-specific line
+ wyt <- function(y, w) { # lag-specific line
+ if("listw" %in% class(w)) { # lag-specific line
+ wyt <- lag.listw(w, y) # lag-specific line
+ } else { # lag-specific line
+ wyt <- w %*% y # lag-specific line
+ } # lag-specific line
+ return(wyt) # lag-specific line
+ } # lag-specific line
wy<-list() # lag-specific line
for (j in 1:length(unique(tind))) { # lag-specific line
yT<-y[tind==unique(tind)[j]] # lag-specific line
- wy[[j]] <- w %*% yT # lag-specific line
+ wy[[j]] <- wyt(yT, w) # lag-specific line
} # lag-specific line
return(unlist(wy)) # lag-specific line
} # lag-specific line
@@ -100,28 +122,65 @@
b.hat <- solve(crossprod(X, sigma.1) %*% X,
crossprod(X, sigma.1) %*% y)
ehat <- y - X %*% b.hat
- sigma2ehat <- (crossprod(ehat, sigma.1) %*% ehat)/(n * t)
+ sigma2ehat <- (crossprod(ehat, sigma.1) %*% ehat)/(n * t.)
return(list(betahat=b.hat, ehat=ehat, sigma2=sigma2ehat))
}
## lag y once for all
- wy <- Wy(y, w2, tind) # lag-specific line
+ wy <- Wy(y, w, tind) # lag-specific line
- ## 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,
+ ## optimization
+
+ ## adaptive scaling
+ parscale <- 1/max(myparms0, 0.1)
+
+ if(method=="nlminb") {
+
+ 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 = parscale,
+ control = list(x.tol = x.tol,
rel.tol = rel.tol, trace = trace),
- lower = lower.bounds, upper = upper.bounds)
+ lower = lower.bounds, upper = upper.bounds)
- ## log likelihood at optimum (notice inverted sign)
- myll <- -optimum$objective
- ## retrieve optimal parms
- myparms <- optimum$par
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms and H
+ myparms <- optimum$par
+ myHessian <- fdHess(myparms, function(x) -ll.c(x,
+ y, X, n, t., w, w2, wy))$Hessian # lag-specific line: wy
+ } else {
+
+ require(maxLik)
+
+ ## initial values are not allowed to be zero
+ maxout<-function(x,a) ifelse(x>a, x, a)
+ myparms0 <- maxout(myparms0, 0.01)
+
+ ## invert sign for MAXimization
+ ll.c2 <- function(phirholambda, y, X, n, t., w, w2, wy) {
+ -ll.c(phirholambda, y, X, n, t., w, w2, wy)
+ }
+
+ ## max likelihood
+ optimum <- maxLik(logLik = ll.c2,
+ grad = NULL, hess = NULL, start=myparms0,
+ method = method,
+ parscale = parscale,
+ constraints=list(ineqA=cA, ineqB=cB),
+ y = y, X = X, n = n, t. = t., w = w, w2 = w2, wy = wy)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- optimum$maximum # this one MAXimizes
+ ## retrieve optimal parms and H
+ myparms <- optimum$estimate
+ myHessian <- optimum$hessian
+ }
+
## one last GLS step at optimal vcov parms
- sigma.1 <- invSigma(myparms, n, t, w)
+ sigma.1 <- invSigma(myparms, n, t., w)
Ay <- y - myparms[length(myparms)] * wy # lag-specific line
beta <- GLSstep(X, Ay, sigma.1)
@@ -130,9 +189,15 @@
solve(crossprod(X, sigma.1) %*% 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) - 1
+ ## error handler here for singular Hessian cases
+ covTheta <- try(solve(-myHessian), silent=TRUE)
+ if(class(covTheta) == "try-error") {
+ covTheta <- matrix(NA, ncol=nvcovpms+1,
+ nrow=nvcovpms+1)
+ warning("Hessian matrix is not invertible")
+ }
covAR <- covTheta[nvcovpms+1, nvcovpms+1, drop=FALSE]
covPRL <- covTheta[1:nvcovpms, 1:nvcovpms, drop=FALSE]
Modified: pkg/R/saremmod.R
===================================================================
--- pkg/R/saremmod.R 2013-03-28 21:12:49 UTC (rev 156)
+++ pkg/R/saremmod.R 2013-03-30 22:22:05 UTC (rev 157)
@@ -1,10 +1,10 @@
saremmod <-
-function (X, y, ind, tind, n, k, t, nT, w, w2, coef0 = rep(0, 2),
+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,
- ...)
+ method="nlminb", ...)
{
- ## extensive function rewriting, Giovanni Millo 29/09/2010
+ ## extensive function rewriting, Giovanni Millo 27/03/2013
## structure:
## a) specific part
## - set names, bounds and initial values for parms
@@ -26,36 +26,31 @@
## initialize values for optimizer
myparms0 <- coef0
- ## set bounds for optimizer
- lower.bounds <- c(-0.999, -0.999) # lag-specific line (2nd parm)
- upper.bounds <- c(0.999, 0.999) # lag-specific line (idem)
## modules for likelihood
- B <- function(lambda, w) diag(1, ncol(w)) - lambda * w
- detB <- function(lambda, w) det(B(lambda, w))
- invSigma <- function(lambdapsi, n, t, w) {
- It <- diag(1, t)
+ invSigma <- function(lambdapsi, n, t., w) {
+ It <- diag(1, t.)
## retrieve parms
lambda <- lambdapsi[1]
## psi not used: here passing 4 parms, but works anyway
## because psi is last one
## calc inverse
- BB <- crossprod(B(lambda, w))
+ BB <- xprodB(lambda, w)
invSigma <- kronecker(It, BB)
invSigma
}
- detSigma <- function(lambda, t, w) {
- detSigma <- t * log(detB(lambda, w))
+ detSigma <- function(lambda, t., w) {
+ detSigma <- t. * ldetB(lambda, w)
detSigma
}
## likelihood function, both steps included
- ll.c <- function(lambdapsi, y, X, n, t, w, w2, wy) {
+ ll.c <- function(lambdapsi, y, X, n, t., w, w2, wy) {
## retrieve parms
lambda <- lambdapsi[1]
psi <- lambdapsi[2] # lag-specific line
## calc inverse sigma
- sigma.1 <- invSigma(lambdapsi, n, t, w)
+ sigma.1 <- invSigma(lambdapsi, n, t., w)
## lag y
Ay <- y - psi * wy # lag-specific line
## do GLS step to get e, s2e
@@ -63,72 +58,134 @@
e <- glsres[["ehat"]]
s2e <- glsres[["sigma2"]]
## calc ll
- zero <- t*log(detB(psi, w2)) # lag-specific line (else zero <- 0)
- due <- detSigma(lambda, t, w)
- tre <- -n * t/2 * log(s2e)
- quattro <- -1/(2 * s2e) * crossprod(e, sigma.1) %*% e
- const <- -(n * t)/2 * log(2 * pi)
+ zero <- t.*ldetB(psi, w2) # lag-specific line (else zero <- 0)
+ due <- detSigma(lambda, t., w)
+ tre <- -n * t./2 * log(s2e)
+ quattro <- -1/(2 * s2e) * t(e) %*% sigma.1 %*% e
+ const <- -(n * t.)/2 * log(2 * pi)
ll.c <- const + zero + due + tre + quattro
## invert sign for minimization
llc <- -ll.c
}
+ ## set bounds for optimizer
+ lower.bounds <- c(-0.999, -0.999) # lag-specific line (4th parm)
+ upper.bounds <- c(0.999, 0.999) # lag-specific line (idem)
+
+ ## constraints as cA %*% theta + cB >= 0
+ ## equivalent to: phi>=0, -1<=(rho, lambda, psi)<=1
+ ## NB in maxLik() optimization cannot start at the boundary of the
+ ## parameter space !
+ cA <- cbind(c(1,-1,rep(0,2)),
+ c(rep(0,2), 1, -1))
+ cB <- c(rep(1,4))
## generic from here
## calc. Wy (spatial lag of y)
+ ## (flexible fun accepting either listws or matrices for w)
Wy <- function(y, w, tind) { # lag-specific line
+ wyt <- function(y, w) { # lag-specific line
+ if("listw" %in% class(w)) { # lag-specific line
+ wyt <- lag.listw(w, y) # lag-specific line
+ } else { # lag-specific line
+ wyt <- w %*% y # lag-specific line
+ } # lag-specific line
+ return(wyt) # lag-specific line
+ } # lag-specific line
wy<-list() # lag-specific line
for (j in 1:length(unique(tind))) { # lag-specific line
yT<-y[tind==unique(tind)[j]] # lag-specific line
- wy[[j]] <- w %*% yT # lag-specific line
+ wy[[j]] <- wyt(yT, w) # lag-specific line
} # lag-specific line
return(unlist(wy)) # lag-specific line
} # lag-specific line
## GLS step function
GLSstep <- function(X, y, sigma.1) {
- b.hat <- solve(crossprod(X, sigma.1) %*% X,
- crossprod(X, sigma.1) %*% y)
+ b.hat <- solve(t(X) %*% sigma.1 %*% X,
+ t(X) %*% sigma.1 %*% y)
ehat <- y - X %*% b.hat
- sigma2ehat <- (crossprod(ehat, sigma.1) %*% ehat)/(n * t)
+ sigma2ehat <- (t(ehat) %*% sigma.1 %*% ehat)/(n * t.)
return(list(betahat=b.hat, ehat=ehat, sigma2=sigma2ehat))
}
## lag y once for all
wy <- Wy(y, w2, tind) # lag-specific line
- ## 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,
+ ## optimization
+
+ ## adaptive scaling
+ parscale <- 1/max(myparms0, 0.1)
+
+ if(method=="nlminb") {
+
+ 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 = parscale,
+ control = list(x.tol = x.tol,
rel.tol = rel.tol, trace = trace),
- lower = lower.bounds, upper = upper.bounds)
+ lower = lower.bounds, upper = upper.bounds)
- ## log likelihood at optimum (notice inverted sign)
- myll <- -optimum$objective
- ## retrieve optimal parms
- myparms <- optimum$par
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms and H
+ myparms <- optimum$par
+ myHessian <- fdHess(myparms, function(x) -ll.c(x,
+ y, X, n, t., w, w2, wy))$Hessian # lag-specific line: wy
+ } else {
+
+ require(maxLik)
+
+ ## initial values are not allowed to be zero
+ maxout<-function(x,a) ifelse(x>a, x, a)
+ myparms0 <- maxout(myparms0, 0.01)
+
+ ## invert sign for MAXimization
+ ll.c2 <- function(phirholambda, y, X, n, t., w, w2, wy) {
+ -ll.c(phirholambda, y, X, n, t., w, w2, wy)
+ }
+
+ ## max likelihood
+ optimum <- maxLik(logLik = ll.c2,
+ grad = NULL, hess = NULL, start=myparms0,
+ method = method,
+ parscale = parscale,
+ constraints=list(ineqA=cA, ineqB=cB),
+ y = y, X = X, n = n, t. = t., w = w, w2 = w2, wy = wy)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- optimum$maximum # this one MAXimizes
+ ## retrieve optimal parms and H
+ myparms <- optimum$estimate
+ myHessian <- optimum$hessian
+ }
+
## one last GLS step at optimal vcov parms
- sigma.1 <- invSigma(myparms, n, t, w)
+ sigma.1 <- invSigma(myparms, n, t., w)
Ay <- y - myparms[length(myparms)] * wy # lag-specific line
beta <- GLSstep(X, Ay, sigma.1)
## final vcov(beta)
covB <- as.numeric(beta[[3]]) *
- solve(crossprod(X, sigma.1) %*% X)
+ solve(t(X) %*% sigma.1 %*% 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) - 1
+ ## error handler here for singular Hessian cases
+ covTheta <- try(solve(-myHessian), silent=TRUE)
+ if(class(covTheta) == "try-error") {
+ covTheta <- matrix(NA, ncol=nvcovpms+1,
+ nrow=nvcovpms+1)
+ warning("Hessian matrix is not invertible")
+ }
covAR <- covTheta[nvcovpms+1, nvcovpms+1, drop=FALSE]
covPRL <- covTheta[1:nvcovpms, 1:nvcovpms, drop=FALSE]
## final parms
- betas <- as.vector(beta[[1]])
- sigma2 <- as.numeric(beta[["sigma2"]])
+ betas <- as.vector(beta[[1]])
+ sigma2 <- as.numeric(beta[["sigma2"]])
arcoef <- myparms[which(nam.errcomp=="lambda")] # lag-specific line
errcomp <- myparms[which(nam.errcomp!="lambda")]
names(betas) <- nam.beta
Modified: pkg/R/saremsrREmod.R
===================================================================
--- pkg/R/saremsrREmod.R 2013-03-28 21:12:49 UTC (rev 156)
+++ pkg/R/saremsrREmod.R 2013-03-30 22:22:05 UTC (rev 157)
@@ -1,7 +1,8 @@
saremsrREmod <-
-function (X, y, ind, tind, n, k, t, nT, w, w2, coef0 = rep(0, 4),
+function (X, y, ind, tind, n, k, t., nT, w, w2, coef0 = rep(0, 4),
hess = FALSE, trace = trace, x.tol = 1.5e-18, rel.tol = 1e-15,
- ...)
+ method="nlminb",
+ ...)
{
## extensive function rewriting, Giovanni Millo 29/09/2010
@@ -17,36 +18,55 @@
## - calc final covariances
## - make list of results
- # mark
- #print("uso versione 5") # fixed 'rho^2' for 'rho' in uno <- n/2 * log(1 - rho^2)
+ ## from version 2: both nlminb (fastest, some negative covariances) and
+ ## optimizers from maxLik ("BFGS", "SANN", "NM") are supported.
+ ## from version 3: analytical inverse for Vmat (see notes; eliminates
+ ## singular matrix pbs., from ca. 45'' to ca. 30'' on 281x3 'datiNY' example)
+ ## and exploit the fact that solve(crossprod(B))=tcrossprod(solve(B))
+ ## (this last not giving any benefit; but check sparse methods on B)
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/splm -r 157
More information about the Splm-commits
mailing list