[Splm-commits] r244 - in pkg: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 26 13:20:17 CEST 2022
Author: the_sculler
Date: 2022-07-26 13:20:17 +0200 (Tue, 26 Jul 2022)
New Revision: 244
Modified:
pkg/DESCRIPTION
pkg/R/sarem2srREmod.R
pkg/R/saremgREmod.R
pkg/R/semgREmod.R
Log:
Added flexible choice of optimizers to GSRE estimators
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2022-07-26 09:39:14 UTC (rev 243)
+++ pkg/DESCRIPTION 2022-07-26 11:20:17 UTC (rev 244)
@@ -1,6 +1,6 @@
Package: splm
Title: Econometric Models for Spatial Panel Data
-Version: 1.6-0
+Version: 1.6-1
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"),
Modified: pkg/R/sarem2srREmod.R
===================================================================
--- pkg/R/sarem2srREmod.R 2022-07-26 09:39:14 UTC (rev 243)
+++ pkg/R/sarem2srREmod.R 2022-07-26 11:20:17 UTC (rev 244)
@@ -103,7 +103,7 @@
upper.bounds <- c(1e+09, 0.999, 0.999, 0.999) # lag-specific line (idem)
## constraints as cA %*% theta + cB >= 0
- ## equivalent to: phi>=0, -1<=(rho, lambda, psi)<=1
+ ## equivalent to: phi>=0, -1<=(psi, rho, lambda)<=1
## NB in maxLik() optimization cannot start at the boundary of the
## parameter space !
cA <- cbind(c(1, rep(0,6)),
Modified: pkg/R/saremgREmod.R
===================================================================
--- pkg/R/saremgREmod.R 2022-07-26 09:39:14 UTC (rev 243)
+++ pkg/R/saremgREmod.R 2022-07-26 11:20:17 UTC (rev 244)
@@ -1,5 +1,5 @@
saremgREmod <-
-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", ...)
{
@@ -30,18 +30,29 @@
## 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)
+ ## constraints as cA %*% theta + cB >= 0
+ ## equivalent to: phi>=0, -1<=(rho1, rho2, lambda)<=1
+ ## NB in maxLik() optimization cannot start at the boundary of the
+ ## parameter space !
+ cA <- cbind(c(1, rep(0,6)),
+ c(0,1,-1,rep(0,4)),
+ c(rep(0,3), 1, -1, rep(0,2)),
+ c(rep(0,5), 1, -1))
+ cB <- c(0, rep(1,6))
+
## 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]
@@ -51,7 +62,7 @@
## because lambda is last one
## calc inverse
invSigma <- kronecker(Jbart,
- solve(t*phi*solve(crossprod(B(rho1,w)))+
+ 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...
@@ -58,12 +69,12 @@
invSigma
}
- detSigma <- function(philambda, t, w) {
+ detSigma <- function(philambda, t., w) {
## used in some formulations
- Jt <- matrix(1, ncol = t, nrow = t)
+ 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
phi <- philambda[1]
@@ -71,36 +82,14 @@
rho2 <- philambda[3]
sigma <- kronecker(Jbart,
- t*phi*solve(crossprod(B(rho1,w))) + solve(crossprod(B(rho2,w)))) +
+ 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) {
+ ll.c <- function(philambda, y, X, n, t., w, w2, wy) {
## retrieve parms
phi <- philambda[1]
rho1 <- philambda[2]
@@ -107,7 +96,7 @@
rho2 <- philambda[3]
lambda <- philambda[4] # lag-specific line
## calc inverse sigma
- sigma.1 <- invSigma(philambda, n, t, w2)
+ sigma.1 <- invSigma(philambda, n, t., w2)
## lag y
Ay <- y - lambda * wy # lag-specific line
## do GLS step to get e, s2e
@@ -115,22 +104,11 @@
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 + ...
-
+ zero <- t.*log(detB(lambda, w)) # lag-specific line (else zero <- 0)
## 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))
+ 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
@@ -155,7 +133,7 @@
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))
}
@@ -173,7 +151,7 @@
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,
+ 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)
@@ -182,15 +160,39 @@
myll <- -optimum$objective
## retrieve optimal parms
myparms <- optimum$par
-
+ myHessian <- fdHess(myparms,
+ function(x) -ll.c(x, y, X, n, t., w, w2, wy))$Hessian
+ # lag-specific line: wy
} else {
- stop("Optim. methods other than 'nlminb' not implemented yet")
+ #stop("Optim. methods other than 'nlminb' not implemented yet")
+ ## 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, w2)
+ 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
@@ -199,15 +201,19 @@
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]
+ #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 # lag-specific line: -1
+ ## error handler here for singular Hessian cases
+ covTheta <- try(solve(-myHessian), silent=TRUE)
+ if(inherits(covTheta, "try-error")) {
+ covTheta <- matrix(NA, ncol=nvcovpms+1, # lag-specific line: +1
+ nrow=nvcovpms+1) # lag-specific line: +1
+ warning("Hessian matrix is not invertible")
+ }
+ covAR <- covTheta[nvcovpms+1, nvcovpms+1, drop=FALSE] # lag-specific line
covPRL <- covTheta[1:nvcovpms, 1:nvcovpms, drop=FALSE]
-
- #covAR <- NULL
- #covPRL <- covTheta
## final parms
betas <- as.vector(beta[[1]])
Modified: pkg/R/semgREmod.R
===================================================================
--- pkg/R/semgREmod.R 2022-07-26 09:39:14 UTC (rev 243)
+++ pkg/R/semgREmod.R 2022-07-26 11:20:17 UTC (rev 244)
@@ -1,5 +1,5 @@
semgREmod <-
-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", ...)
{
@@ -7,6 +7,9 @@
## Trieste, 31/8/2021. From the materials to the semregmod paper
## (Millo, J spat Econometrics, 2022).
+ ## general errors a la
+ ## Baltagi, Egger and Pfaffermayr
+
## 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.
@@ -24,9 +27,6 @@
## - 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")
@@ -33,18 +33,28 @@
## 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)
+
+ ## constraints as cA %*% theta + cB >= 0
+ ## equivalent to: phi>=0, -1<=(rho1, rho2)<=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))
## 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]
@@ -54,7 +64,7 @@
## because psi is last one
## calc inverse
invSigma <- kronecker(Jbart,
- solve(t*phi*solve(crossprod(B(rho1,w)))+
+ 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...
@@ -61,12 +71,12 @@
invSigma
}
- detSigma <- function(philambda, t, w) {
+ detSigma <- function(philambda, t., w) {
## used in some formulations
- Jt <- matrix(1, ncol = t, nrow = t)
+ 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
phi <- philambda[1]
@@ -74,7 +84,7 @@
rho2 <- philambda[3]
sigma <- kronecker(Jbart,
- t*phi*solve(crossprod(B(rho1,w))) + solve(crossprod(B(rho2,w)))) +
+ t.*phi*solve(crossprod(B(rho1,w))) + solve(crossprod(B(rho2,w)))) +
kronecker(Et, solve(crossprod(B(rho2,w))))
detSigma <- det(sigma)
detSigma
@@ -81,34 +91,23 @@
}
## 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]
rho1 <- philambda[2]
rho2 <- philambda[3]
## calc inverse sigma
- sigma.1 <- invSigma(philambda, n, t, w2)
+ 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))
+ 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
@@ -121,7 +120,7 @@
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))
}
@@ -139,7 +138,7 @@
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,
+ 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)
@@ -148,16 +147,42 @@
myll <- -optimum$objective
## retrieve optimal parms
myparms <- optimum$par
+ myHessian <- fdHess(myparms,
+ function(x) -ll.c(x, y, X, n, t., w, w2, wy))$Hessian
+ # lag-specific line: wy
} else {
- stop("Optim. methods other than 'nlminb' not implemented yet")
+ #stop("Optim. methods other than 'nlminb' not implemented yet")
+
+ ## 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)
beta <- GLSstep(X, y, sigma.1)
## final vcov(beta)
@@ -165,8 +190,17 @@
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
+# 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)
+ ## error handler here for singular Hessian cases
+ covTheta <- try(solve(-myHessian), silent=TRUE)
+ if(inherits(covTheta, "try-error")) {
+ covTheta <- matrix(NA, ncol=nvcovpms,
+ nrow=nvcovpms)
+ warning("Hessian matrix is not invertible")
+ }
covAR <- NULL
covPRL <- covTheta
More information about the Splm-commits
mailing list