[Splm-commits] r243 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 26 11:39:14 CEST 2022
Author: the_sculler
Date: 2022-07-26 11:39:14 +0200 (Tue, 26 Jul 2022)
New Revision: 243
Added:
pkg/R/saremgREmod.R
pkg/R/semgREmod.R
Modified:
pkg/DESCRIPTION
pkg/R/spreml.R
pkg/man/spreml.Rd
Log:
Added new generalized spatial random effects (GSRE) estimator from JSPE paper
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2022-07-26 07:12:13 UTC (rev 242)
+++ pkg/DESCRIPTION 2022-07-26 09:39:14 UTC (rev 243)
@@ -1,7 +1,7 @@
Package: splm
Title: Econometric Models for Spatial Panel Data
-Version: 1.5-5
-Date: 2022-07-24
+Version: 1.6-0
+Date: 2022-07-26
Authors at R: c(person(given = "Giovanni", family = "Millo", role = c("aut", "cre"), email = "giovanni.millo at deams.units.it"),
person(given = "Gianfranco", family = "Piras", role = c("aut"), email = "gpiras at mac.com"),
person("Roger", "Bivand", role = c("ctb"), email = "Roger.Bivand at nhh.no", comment=c(ORCID="0000-0003-2392-6140")))
Added: pkg/R/saremgREmod.R
===================================================================
--- pkg/R/saremgREmod.R (rev 0)
+++ pkg/R/saremgREmod.R 2022-07-26 09:39:14 UTC (rev 243)
@@ -0,0 +1,229 @@
+saremgREmod <-
+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", ...)
+{
+
+ ## Trieste, 22/3/2022. This is the SAR extension of semgREmod2 from the
+ ## materials to the semregmod paper (Millo, J spat Econometrics, 2022)
+
+ ## (summary from general version)
+ ## extensive function rewriting, Giovanni Millo 29/09/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
+
+ ## this version sem2gre (orig.: 14/10/2010): general errors
+ ## a la Baltagi, Egger and Pfaffermayr
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("phi", "rho1", "rho2", "lambda")
+
+ ## initialize values for optimizer
+ myparms0 <- coef0
+ ## set bounds for optimizer
+ lower.bounds <- c(1e-08, -0.999, -0.999, -0.999) # lag-specific line (4th parm)
+ upper.bounds <- c(1e+09, 0.999, 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)
+ #In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+ ## retrieve parms
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+ ## lambda not used: here passing 4 parms, but works anyway
+ ## because lambda is last one
+ ## calc inverse
+ invSigma <- kronecker(Jbart,
+ solve(t*phi*solve(crossprod(B(rho1,w)))+
+ solve(crossprod(B(rho2,w))))) +
+ kronecker(Et, crossprod(B(rho2,w))) # fixed error:
+ # was solve(crossprod(B...
+
+ invSigma
+ }
+ detSigma <- function(philambda, t, w) {
+ ## used in some formulations
+ Jt <- matrix(1, ncol = t, nrow = t)
+ #In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+
+ sigma <- kronecker(Jbart,
+ t*phi*solve(crossprod(B(rho1,w))) + solve(crossprod(B(rho2,w)))) +
+ kronecker(Et, solve(crossprod(B(rho2,w))))
+ detSigma <- det(sigma)
+ detSigma
+ }
+
+ detOmega <- function(philambda, t, w, s2e) {
+ ## used in some formulations; probably wrong!
+ Jt <- matrix(1, ncol = t, nrow = t)
+ #In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+
+ s2e <- as.numeric(s2e)
+
+ omega <- kronecker(Jbart,
+ t * (phi*s2e) * solve(crossprod(B(rho1,w))) +
+ s2e * solve(crossprod(B(rho2,w)))) +
+ s2e * kronecker(Et, solve(crossprod(B(rho2,w))))
+ detOmega <- det(omega)
+ detOmega
+ }
+
+ ## likelihood function, both steps included
+ ll.c <- function(philambda, y, X, n, t, w, w2, wy) {
+ ## retrieve parms
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+ lambda <- philambda[4] # lag-specific line
+ ## calc inverse sigma
+ sigma.1 <- invSigma(philambda, n, t, w2)
+ ## lag y
+ Ay <- y - lambda * wy # lag-specific line
+ ## do GLS step to get e, s2e
+ glsres <- GLSstep(X, Ay, sigma.1) # lag-specific line (Ay for y)
+ e <- glsres[["ehat"]]
+ s2e <- glsres[["sigma2"]]
+ ## calc ll
+ zero <- t*log(detB(lambda, w)) # lag-specific line (else zero <- 0)
+ #due <- detSigma(phi, lambda, n, t, w)
+ # this is log(detSigma)
+ due <- -1/2*log(det(t*phi*solve(crossprod(B(rho1,w2)))+
+ solve(crossprod(B(rho2,w2))))) - # was: +
+ (t-1)/2*log(detB(rho2,w2)) # was: (t-1)
+ tre <- -(n*t)/2 * log(s2e)
+ quattro <- -1/(2 * s2e) * crossprod(e, sigma.1) %*% e # tre or quattro are alternative!
+ const <- -(n * t)/2 * log(2 * pi)
+ ## log-likelihood "the usual way", from BEP paper
+ #ll.c <- const + zero + due + quattro # was: ... + tre + ...
+
+ ## ex appendix to BEP, p.6, concentrated ll
+ ## (notice that, from GLS step, s2e = 1/nt * crossprod(e, sigma.1) %*% e)
+ ll.c <- zero - (n*t)/2 * (log(2*pi) + 1) - (n*t)/2 * log(s2e) - # lag-specific
+ 1/2 * log(detSigma(philambda, t, w2))
+ # or: log(detOmega(philambda, t, w, s2e))
+
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## calc. Wy (spatial lag of y)
+ ## NB problems with NAs on y, see examples_sar.R in /Articoli/semregmod
+ Wy <- function(y, w, tind) { # 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
+ } # 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)
+ ehat <- y - X %*% b.hat
+ 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, w, tind) # lag-specific line
+
+ ## max likelihood
+
+ ## 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 = 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
+
+ } else {
+
+ stop("Optim. methods other than 'nlminb' not implemented yet")
+
+ }
+
+ ## one last GLS step at optimal vcov parms
+ sigma.1 <- invSigma(myparms, n, t, w2)
+ Ay <- y - myparms[length(myparms)] * wy # lag-specific line
+ beta <- GLSstep(X, Ay, sigma.1) # lag-specific line
+
+ ## final vcov(beta)
+ covB <- as.numeric(beta[[3]]) *
+ 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
+ covAR <- covTheta[nvcovpms+1, nvcovpms+1, drop=FALSE]
+ covPRL <- covTheta[1:nvcovpms, 1:nvcovpms, drop=FALSE]
+
+ #covAR <- NULL
+ #covPRL <- covTheta
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- myparms[which(nam.errcomp=="lambda")] # lag-specific line
+ errcomp <- myparms[which(nam.errcomp!="lambda")]
+ names(betas) <- nam.beta
+ names(arcoef) <- "lambda" # lag-specific line
+ names(errcomp) <- nam.errcomp[which(nam.errcomp!="lambda")]
+
+ dimnames(covB) <- list(nam.beta, nam.beta)
+ dimnames(covAR) <- list(names(arcoef), names(arcoef))
+ 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)
+}
Added: pkg/R/semgREmod.R
===================================================================
--- pkg/R/semgREmod.R (rev 0)
+++ pkg/R/semgREmod.R 2022-07-26 09:39:14 UTC (rev 243)
@@ -0,0 +1,188 @@
+semgREmod <-
+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", ...)
+{
+
+ ## Trieste, 31/8/2021. From the materials to the semregmod paper
+ ## (Millo, J spat Econometrics, 2022).
+
+ ## This is a fix of the version from 29/09/2010, had (B'B)^(-1) instead
+ ## of just B'B in the last element of invSigma(). Works now.
+
+ ## (summary from previous version)
+ ## extensive function rewriting, Giovanni Millo 29/09/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
+
+ ## general errors a la
+ ## Baltagi, Egger and Pfaffermayr
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("phi", "rho1", "rho2")
+
+ ## initialize values for optimizer
+ myparms0 <- coef0
+ ## set bounds for optimizer
+ lower.bounds <- c(1e-08, -0.999, -0.999)
+ upper.bounds <- c(1e+09, 0.999, 0.999)
+
+ ## 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)
+ #In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+ ## retrieve parms
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+ ## psi not used: here passing 4 parms, but works anyway
+ ## because psi is last one
+ ## calc inverse
+ invSigma <- kronecker(Jbart,
+ solve(t*phi*solve(crossprod(B(rho1,w)))+
+ solve(crossprod(B(rho2,w))))) +
+ kronecker(Et, crossprod(B(rho2,w))) # fixed error:
+ # was solve(crossprod(B...
+
+ invSigma
+ }
+ detSigma <- function(philambda, t, w) {
+ ## used in some formulations
+ Jt <- matrix(1, ncol = t, nrow = t)
+ #In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+
+ sigma <- kronecker(Jbart,
+ t*phi*solve(crossprod(B(rho1,w))) + solve(crossprod(B(rho2,w)))) +
+ kronecker(Et, solve(crossprod(B(rho2,w))))
+ detSigma <- det(sigma)
+ detSigma
+ }
+
+ ## likelihood function, both steps included
+ ll.c <- function(philambda, y, X, n, t, w, w2, wy) {
+ ## retrieve parms
+ phi <- philambda[1]
+ rho1 <- philambda[2]
+ rho2 <- philambda[3]
+ ## calc inverse sigma
+ sigma.1 <- invSigma(philambda, n, t, w2)
+ ## do GLS step to get e, s2e
+ glsres <- GLSstep(X, y, sigma.1)
+ e <- glsres[["ehat"]]
+ s2e <- glsres[["sigma2"]]
+ ## calc ll
+ #due <- detSigma(phi, lambda, n, t, w)
+
+ # this is log(detSigma)
+ # due <- -1/2*log(det(t*phi*solve(crossprod(B(rho1,w2)))+
+ # solve(crossprod(B(rho2,w2))))) - # was: +
+ # (t-1)/2*log(detB(rho2,w2)) # was: (t-1)
+ # tre <- -(n*t)/2 * log(s2e)
+ # quattro <- -1/(2 * s2e) * crossprod(e, sigma.1) %*% e # tre or quattro are alternative!
+ # const <- -(n * t)/2 * log(2 * pi)
+ ## log-likelihood "the usual way", from BEP paper
+ #ll.c <- const + due + quattro # was: ... + tre + ...
+
+ ## ex appendix to BEP, p.6, concentrated ll
+ ## (notice that, from GLS step, s2e = 1/nt * crossprod(e, sigma.1) %*% e)
+ ll.c <- - (n*t)/2 * (log(2*pi) + 1) - (n*t)/2 * log(s2e) - 1/2 *
+ log(detSigma(philambda, t, w2))
+
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## GLS step function
+ GLSstep <- function(X, y, sigma.1) {
+ 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)
+ return(list(betahat=b.hat, ehat=ehat, sigma2=sigma2ehat))
+ }
+
+ ## lag y parm kept for compatibility
+ wy <- NULL
+
+ ## max likelihood
+
+ ## 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 = 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
+
+ } else {
+
+ stop("Optim. methods other than 'nlminb' not implemented yet")
+
+ }
+
+
+ ## one last GLS step at optimal vcov parms
+ sigma.1 <- invSigma(myparms, n, t, w)
+ beta <- GLSstep(X, y, sigma.1)
+
+ ## final vcov(beta)
+ covB <- as.numeric(beta[[3]]) *
+ 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
+ covAR <- NULL
+ covPRL <- covTheta
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- NULL
+ errcomp <- myparms[which(nam.errcomp!="psi")]
+ names(betas) <- nam.beta
+ names(errcomp) <- nam.errcomp
+
+ 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)
+}
Modified: pkg/R/spreml.R
===================================================================
--- pkg/R/spreml.R 2022-07-26 07:12:13 UTC (rev 242)
+++ pkg/R/spreml.R 2022-07-26 09:39:14 UTC (rev 243)
@@ -1,7 +1,8 @@
spreml <-
function (formula, data, index = NULL, w, w2=w, lag = FALSE,
errors = c("semsrre", "semsr", "srre", "semre",
- "re", "sr", "sem","ols", "sem2srre", "sem2re"),
+ "re", "sr", "sem","ols", "sem2srre",
+ "sem2re", "semgre"),
pvar = FALSE, hess = FALSE, quiet = TRUE,
initval = c("zeros", "estimate"),
x.tol = 1.5e-18, rel.tol = 1e-15, ...)
@@ -61,7 +62,7 @@
## manage initial values
sv.length <- switch(match.arg(errors), semsrre = 3, semsr = 2,
srre = 2, semre = 2, re = 1, sr = 1, sem = 1, ols = 0,
- sem2srre = 3, sem2re = 2)
+ sem2srre = 3, sem2re = 2, semgre = 3)
errors. <- match.arg(errors)
if (is.numeric(initval)) {
if (length(initval) != sv.length) {
@@ -73,29 +74,43 @@
switch(match.arg(initval), zeros = {
coef0 <- rep(0, sv.length)
}, estimate = {
- if (nchar(errors.) < 4) {
- stop("Pre-estimation of unique vcov parm is meaningless: \n please select (default) option 'zeros' or supply a scalar")
+ if(errors. == "semgre") {
+
+ ## if errors are GSRE the pre-estimation is not implemented
+ ## as there would be little gains
+ coef0 <- rep(0, sv.length)
+ warning("Pre-estimation of error components is not implemented for 'semgre': starting values were set to zero.")
+
+ } else {
+
+ ## determine individual error parms to be pre-estimated and do
+ ## estimation of single.parameter models
+
+ if (nchar(errors.) < 4) {
+ stop("Pre-estimation of unique vcov parm is meaningless: \n please select (default) option 'zeros' or supply a scalar")
+ }
+ coef0 <- NULL
+ if (grepl("re", errors.)) {
+ REmodel <- REmod(X, y, ind, tind, n, k, t, nT,
+ w, coef0 = 0, hess = FALSE, trace = trace,
+ x.tol = 1.5e-18, rel.tol = 1e-15, ...)
+ coef0 <- c(coef0, REmodel$errcomp)
+ }
+ if (grepl("sr", errors.)) {
+ ARmodel <- ssrmod(X, y, ind, tind, n, k, t, nT,
+ w, coef0 = 0, hess = FALSE, trace = trace,
+ x.tol = 1.5e-18, rel.tol = 1e-15, ...)
+ coef0 <- c(coef0, ARmodel$errcomp)
+ }
+ if (grepl("sem", errors.)) {
+ SEMmodel <- semmod(X, y, ind, tind, n, k, t,
+ nT, w, coef0 = 0, hess = FALSE, trace = trace,
+ x.tol = 1.5e-18, rel.tol = 1e-15, ...)
+ coef0 <- c(coef0, SEMmodel$errcomp)
+ }
}
- coef0 <- NULL
- if (grepl("re", errors.)) {
- REmodel <- REmod(X, y, ind, tind, n, k, t, nT,
- w, coef0 = 0, hess = FALSE, trace = trace,
- x.tol = 1.5e-18, rel.tol = 1e-15, ...)
- coef0 <- c(coef0, REmodel$errcomp)
- }
- if (grepl("sr", errors.)) {
- ARmodel <- ssrmod(X, y, ind, tind, n, k, t, nT,
- w, coef0 = 0, hess = FALSE, trace = trace,
- x.tol = 1.5e-18, rel.tol = 1e-15, ...)
- coef0 <- c(coef0, ARmodel$errcomp)
- }
- if (grepl("sem", errors.)) {
- SEMmodel <- semmod(X, y, ind, tind, n, k, t,
- nT, w, coef0 = 0, hess = FALSE, trace = trace,
- x.tol = 1.5e-18, rel.tol = 1e-15, ...)
- coef0 <- c(coef0, SEMmodel$errcomp)
- }
- })
+ }
+ )
}
## switch actual estimator function
@@ -120,6 +135,8 @@
sarmod
}, sem2re = {
sarem2REmod
+ }, semgre = {
+ saremgREmod
})
coef0 <- c(coef0, 0)
} else {
@@ -143,6 +160,8 @@
olsmod
}, sem2re = {
sem2REmod
+ }, semgre = {
+ semgREmod
})
arcoef <- NULL
}
Modified: pkg/man/spreml.Rd
===================================================================
--- pkg/man/spreml.Rd 2022-07-26 07:12:13 UTC (rev 242)
+++ pkg/man/spreml.Rd 2022-07-26 09:39:14 UTC (rev 243)
@@ -8,7 +8,8 @@
\usage{
spreml(formula, data, index = NULL, w, w2=w, lag = FALSE,
errors = c("semsrre", "semsr", "srre", "semre",
- "re", "sr", "sem","ols", "sem2srre", "sem2re"),
+ "re", "sr", "sem","ols", "sem2srre",
+ "sem2re", "semgre"),
pvar = FALSE, hess = FALSE, quiet = TRUE,
initval = c("zeros", "estimate"),
x.tol = 1.5e-18, rel.tol = 1e-15, ...)
More information about the Splm-commits
mailing list