[Splm-commits] r92 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 15 00:01:04 CET 2010
Author: the_sculler
Date: 2010-12-15 00:01:04 +0100 (Wed, 15 Dec 2010)
New Revision: 92
Added:
pkg/R/lrtest.splm.R
pkg/R/olsmod.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/semsrREmod.R
pkg/R/semsrmod.R
pkg/R/sphtest.R
Log:
Added more new comp. engine functions of spreml
Added: pkg/R/lrtest.splm.R
===================================================================
--- pkg/R/lrtest.splm.R (rev 0)
+++ pkg/R/lrtest.splm.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,61 @@
+lrtest.splm <- function(x, y, ...) {
+ ## correspondence with component names
+ compnames <- c("RE", "AR(1)", "SEM", "SAR")
+ compnames.long <- c("random effects", "AR(1) errors",
+ "spatial errors", "spatial lag")
+ names(compnames) <- c("phi", "rho", "lambda", "psi")
+
+ ## first: check same betas!
+ if(!identical(names(x$coef), names(y$coef))) {
+ stop("Models estimated on different regressors")
+ }
+ ## ...and Nobs
+ if(!identical(length(x$residuals), length(y$residuals))) {
+ stop("Models estimated on different number of obs.")
+ }
+
+ ## build error and ev. SAR components vector
+ ecompsx <- c(names(x$errcomp), names(x$arcoef))
+ ecompsy <- c(names(y$errcomp), names(y$arcoef))
+
+ ## check that models be different
+ if(identical(ecompsx, ecompsy)) stop("The model is the same")
+
+ ## which model is bigger?
+ if(length(ecompsx)>=length(ecompsy)) {
+ m1 <- x
+ m0 <- y
+ ecomps1 <- ecompsx
+ ecomps0 <- ecompsy
+ } else {
+ m1 <- y
+ m0 <- x
+ ecomps1 <- ecompsy
+ ecomps0 <- ecompsx
+ }
+
+ ## check if nested
+ if(!all(ecompsx %in%ecompsy) & !all(ecompsy %in% ecompsx)) {
+ stop("Models are not nested")
+ }
+
+ ll1 <- m1$logLik
+ ll0 <- m0$logLik
+
+ testedparms <- ecomps1[-which(ecomps1 %in% ecomps0)]
+
+ LRstat <- 2*(ll1-ll0)
+ df <- abs(length(ecompsx) - length(ecompsy))
+ names(df) <- "df"
+ pLR <- pchisq(LRstat, df = df, lower.tail=FALSE)
+
+ names(LRstat) <- "LR test"
+ RVAL <- list(statistic = LRstat, parameter = df,
+ method = paste("Likelihood ratio for exclusion of",
+ paste(compnames[testedparms], collapse = ", "), "\n\n",
+ "from panel model with ",
+ paste(compnames[ecomps1], collapse = ", ")),
+ p.value = pLR, data.name = "dname here (see bgtest)")
+ class(RVAL) <- "htest"
+ return(RVAL)
+}
Property changes on: pkg/R/lrtest.splm.R
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/R/olsmod.R
===================================================================
--- pkg/R/olsmod.R (rev 0)
+++ pkg/R/olsmod.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,100 @@
+olsmod <-
+function (X, y, ind, tind, n, k, t, nT, w, w2, coef0 = rep(0, dim(X)[[2]]),
+ hess = FALSE, trace = trace, x.tol = 1.5e-18, rel.tol = 1e-15,
+ ...)
+{
+
+ ## 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 function just for compatibility:
+ ## ML estimation of OLS model; produces (hopefully same)
+ ## results as lm() and logLik.lm but as a 'splm' object
+
+ ## good for any GLS estimation by ML, just supply the
+ ## right sigma.1() function
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ #nam.errcomp <- c("psi")
+
+ ## initialize values for optimizer
+ ## (NULL is currently passed as coef0 if ols and lag=F)
+ coef0 = rep(0, dim(X)[[2]])
+ myparms0 <- coef0
+ ## set bounds for optimizer
+ lower.bounds <- -Inf
+ upper.bounds <- Inf
+
+ ## modules for likelihood
+ ## (none)
+
+ ## likelihood function, both steps included
+ ll.c <- function(betas, y, X, n, t, w, w2, wy) {
+ ## get e, s2e as function of betas
+ e <- y - X %*% betas # lag-specific line (Ay for y)
+ s2e <- crossprod(e)/(n*t)
+ ## calc ll
+ tre <- -n * t/2 * log(s2e)
+ quattro <- -1/(2 * s2e) * crossprod(e)
+ const <- -(n * t)/2 * log(2 * pi)
+ ll.c <- const + tre + quattro
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+
+ ## GLS step function suppressed
+
+ ## max likelihood
+ optimum <- nlminb(start = myparms0, objective = ll.c,
+ gradient = NULL, hessian = NULL,
+ y = y, X = X, n = n, t = t, w = w, w2 = w2, wy = wy,
+ scale = 1, control = list(x.tol = x.tol,
+ rel.tol = rel.tol, trace = trace),
+ lower = lower.bounds, upper = upper.bounds)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms
+ betas <- optimum$par
+
+ ## one last GLS step at optimal vcov parms suppressed
+
+ ## final vcov(beta)
+ e <- y - X %*% betas # lag-specific line (Ay for y)
+ s2e <- crossprod(e)/(n*t)
+ covB <- as.numeric(s2e) * solve(crossprod(X))
+
+ ## final vcov(errcomp)
+ covAR <- NULL # ols.errors-specific
+ covPRL <- NULL # ols.errors-specific
+
+ ## final parms
+ #betas ok
+ arcoef <- NULL # ols.errors-specific line
+ errcomp <- NULL # ols.errors-specific
+ names(betas) <- nam.beta
+ #names(arcoef) <- "psi" # lag-specific line
+ #names(errcomp) <- nam.errcomp[which(nam.errcomp!="psi")]
+
+ 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)
+}
Property changes on: pkg/R/olsmod.R
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/R/sarem2REmod.R
===================================================================
--- pkg/R/sarem2REmod.R (rev 0)
+++ pkg/R/sarem2REmod.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,158 @@
+sarem2REmod <-
+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,
+ ...)
+{
+
+ ## 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
+
+ # mark
+ #print("uso versione 0") # done from saremsrREmod4.R
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("phi", "lambda", "psi")
+
+ ## 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)
+ #In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+ ## retrieve parms
+ phi <- philambda[1]
+ lambda <- philambda[2]
+ ## 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 )
+ invSigma
+ }
+ 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
+ Et <- It - Jbart
+ detSigma <- -n/2*log( det( (t*phi+1) * Jbart + Et) ) +
+ t*log(detB(lambda, w))
+ detSigma
+ }
+
+ ## likelihood function, both steps included
+ 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)
+ ## lag y
+ Ay <- y - psi * 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(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)
+ ll.c <- const + zero + due + tre + quattro
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## calc. Wy (spatial lag of y)
+ 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, 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,
+ rel.tol = rel.tol, trace = trace),
+ lower = lower.bounds, upper = upper.bounds)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms
+ myparms <- optimum$par
+
+ ## one last GLS step at optimal vcov parms
+ sigma.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)
+
+ ## 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]
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- myparms[which(nam.errcomp=="psi")] # lag-specific line
+ errcomp <- myparms[which(nam.errcomp!="psi")]
+ names(betas) <- nam.beta
+ names(arcoef) <- "psi" # lag-specific line
+ names(errcomp) <- nam.errcomp[which(nam.errcomp!="psi")]
+
+ 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)
+}
Property changes on: pkg/R/sarem2REmod.R
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/R/saremREmod.R
===================================================================
--- pkg/R/saremREmod.R (rev 0)
+++ pkg/R/saremREmod.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,156 @@
+saremREmod <-
+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,
+ ...)
+{
+
+ ## 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
+
+ # mark
+ #print("uso versione 0") # done from saremsrREmod4.R
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("phi", "lambda", "psi")
+
+ ## 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)
+ In <- diag(1, n)
+ It <- diag(1, t)
+ Jbart <- Jt/t
+ Et <- It - Jbart
+ ## retrieve parms
+ phi <- philambda[1]
+ lambda <- philambda[2]
+ ## 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)
+ invSigma
+ }
+ 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
+ }
+
+ ## likelihood function, both steps included
+ 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)
+ ## lag y
+ Ay <- y - psi * 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(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)
+ ll.c <- const + zero + due + tre + quattro
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## calc. Wy (spatial lag of y)
+ 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, 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,
+ rel.tol = rel.tol, trace = trace),
+ lower = lower.bounds, upper = upper.bounds)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms
+ myparms <- optimum$par
+
+ ## one last GLS step at optimal vcov parms
+ sigma.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)
+
+ ## 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]
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- myparms[which(nam.errcomp=="psi")] # lag-specific line
+ errcomp <- myparms[which(nam.errcomp!="psi")]
+ names(betas) <- nam.beta
+ names(arcoef) <- "psi" # lag-specific line
+ names(errcomp) <- nam.errcomp[which(nam.errcomp!="psi")]
+
+ 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)
+}
Property changes on: pkg/R/saremREmod.R
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/R/saremmod.R
===================================================================
--- pkg/R/saremmod.R (rev 0)
+++ pkg/R/saremmod.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,146 @@
+saremmod <-
+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,
+ ...)
+{
+
+ ## 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
+
+ # mark
+ #print("uso versione 0") # done from saremsrREmod4.R
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("lambda", "psi")
+
+ ## 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)
+ ## 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))
+ invSigma <- kronecker(It, BB)
+ invSigma
+ }
+ detSigma <- function(lambda, t, w) {
+ detSigma <- t * log(detB(lambda, w))
+ detSigma
+ }
+
+ ## likelihood function, both steps included
+ 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)
+ ## lag y
+ Ay <- y - psi * 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(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)
+ ll.c <- const + zero + due + tre + quattro
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## calc. Wy (spatial lag of y)
+ 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, 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,
+ rel.tol = rel.tol, trace = trace),
+ lower = lower.bounds, upper = upper.bounds)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms
+ myparms <- optimum$par
+
+ ## one last GLS step at optimal vcov parms
+ sigma.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)
+
+ ## 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]
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- myparms[which(nam.errcomp=="psi")] # lag-specific line
+ errcomp <- myparms[which(nam.errcomp!="psi")]
+ names(betas) <- nam.beta
+ names(arcoef) <- "psi" # lag-specific line
+ names(errcomp) <- nam.errcomp[which(nam.errcomp!="psi")]
+
+ 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)
+}
Property changes on: pkg/R/saremmod.R
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/R/saremsrREmod.R
===================================================================
--- pkg/R/saremsrREmod.R (rev 0)
+++ pkg/R/saremsrREmod.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,164 @@
+saremsrREmod <-
+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,
+ ...)
+{
+
+ ## 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
+
+ # mark
+ #print("uso versione 5") # fixed 'rho^2' for 'rho' in uno <- n/2 * log(1 - rho^2)
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("phi", "rho", "lambda", "psi")
+
+ ## 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
+ Vmat <- function(rho, t) {
+ V1 <- matrix(ncol = t, nrow = t)
+ for (i in 1:t) V1[i, ] <- rho^abs(1:t - i)
+ V <- (1/(1 - rho^2)) * V1
+ }
+ B <- function(lambda, w) diag(1, ncol(w)) - lambda * w
+ detB <- function(lambda, w) det(B(lambda, w))
+ alfa2 <- function(rho) (1 + rho)/(1 - rho)
+ d2 <- function(rho, t) alfa2(rho) + t - 1
+ Jt <- matrix(1, ncol = t, nrow = t)
+ In <- diag(1, n)
+ det2 <- function(phi, rho, lambda, t) det(d2(rho, t) * (1 -
+ rho)^2 * phi * In + solve(crossprod(B(lambda, w))))
+ Z0 <- function(phi, rho, lambda, t) solve(d2(rho, t) * (1 -
+ rho)^2 * phi * In + solve(crossprod(B(lambda, w))))
+ invSigma <- function(phirholambda, n, t, w) {
+ ## retrieve parms
+ phi <- phirholambda[1]
+ rho <- phirholambda[2]
+ lambda <- phirholambda[3]
+ ## psi not used: here passing 4 parms, but works anyway
+ ## because psi is last one
+ ## calc inverse
+ invVmat <- solve(Vmat(rho, t))
+ BB <- crossprod(B(lambda, w))
+ invSi1 <- kronecker(invVmat, BB)
+ invSi2 <- 1/(d2(rho, t) * (1 - rho)^2)
+ invSi3 <- kronecker(solve(Vmat(rho, t), Jt) %*% invVmat,
+ Z0(phi, rho, lambda, t) - BB)
+ invSigma <- invSi1 + invSi2 * invSi3
+ invSigma
+ }
+ ## likelihood function, both steps included
+ ll.c <- function(phirholambda, y, X, n, t, w, w2, wy) {
+ ## retrieve parms
+ phi <- phirholambda[1]
+ rho <- phirholambda[2]
+ lambda <- phirholambda[3]
+ psi <- phirholambda[4] # lag-specific line
+ ## calc inverse sigma
+ sigma.1 <- invSigma(phirholambda, n, t, w)
+ ## lag y
+ Ay <- y - psi * 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(psi, w2)) # lag-specific line (else zero <- 0)
+ uno <- n/2 * log(1 - rho^2)
+ due <- -1/2 * log(det2(phi, rho, lambda, t))
+ tre <- -(n * t)/2 * log(s2e)
+ quattro <- (t - 1) * log(detB(lambda, w))
+ cinque <- -1/(2 * s2e) * crossprod(e, sigma.1) %*% e
+ const <- -(n * t)/2 * log(2 * pi)
+ ll.c <- const + zero + uno + due + tre + quattro + cinque
+ ## invert sign for minimization
+ llc <- -ll.c
+ }
+
+ ## generic from here
+
+ ## calc. Wy (spatial lag of y)
+ 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, 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,
+ rel.tol = rel.tol, trace = trace),
+ lower = lower.bounds, upper = upper.bounds)
+
+ ## log likelihood at optimum (notice inverted sign)
+ myll <- -optimum$objective
+ ## retrieve optimal parms
+ myparms <- optimum$par
+
+ ## one last GLS step at optimal vcov parms
+ sigma.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)
+
+ ## 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]
+
+ ## final parms
+ betas <- as.vector(beta[[1]])
+ arcoef <- myparms[which(nam.errcomp=="psi")] # lag-specific line
+ errcomp <- myparms[which(nam.errcomp!="psi")]
+ names(betas) <- nam.beta
+ names(arcoef) <- "psi" # lag-specific line
+ names(errcomp) <- nam.errcomp[which(nam.errcomp!="psi")]
+
+ 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)
+}
Property changes on: pkg/R/saremsrREmod.R
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/R/saremsrmod.R
===================================================================
--- pkg/R/saremsrmod.R (rev 0)
+++ pkg/R/saremsrmod.R 2010-12-14 23:01:04 UTC (rev 92)
@@ -0,0 +1,152 @@
+saremsrmod <-
+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,
+ ...)
+{
+
+ ## 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
+
+ # mark
+ #print("uso versione 0")
+
+ ## set names for final parms vectors
+ nam.beta <- dimnames(X)[[2]]
+ nam.errcomp <- c("rho", "lambda", "psi")
+
+ ## initialize values for optimizer
+ myparms0 <- coef0
+ ## set bounds for optimizer
+ lower.bounds <- c(-0.999, -0.999, -0.999) # lag-specific line (4th parm)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/splm -r 92
More information about the Splm-commits
mailing list