Nothing
Infin <- 1e16
funcLogL <- function(parameters, quiet = F){
st <- Sys.time()
alg <- GenzBretz(abseps = 1e-200)
#logging("Algorithm: GenzBretz(maxpts = 1e5, abseps = 1e-12)")
y <- envirGet("y")
X <- envirGet("X")
k <- ncol(X)
W_y <- envirGet("W_y")
W_v <- envirGet("W_v")
W_u <- envirGet("W_u")
costFrontier <- envirGet("costFrontier")
inefficiency <- envirGet("inefficiency")
isSpY <- !is.null(W_y)
isSpV <- !is.null(W_v)
isSpU <- !is.null(W_u)
isTN <- (inefficiency == "truncated")
costf <- ifelse(costFrontier, -1, 1)
p <- paramsFromVector(parameters, k, isSpY, isSpV, isSpU, isTN)
if(!quiet){
pbeta <- olsenReparamBack(p)
paramsToVector(pbeta, olsen = F)
counter = envirCounter("funcLogL")
logging(paste("Evaluating 'funcLogL' for parameters","[run=",counter,"]:"),paramsToVector(pbeta, olsen = F))
}
n <- length(y)
I <- diag(n)
ret <- -Inf
if ((p$nu>0) && (p$lambda>0) && (is.null(p$rhoY) || abs(p$rhoY)<1)
&& (is.null(p$rhoV) || abs(p$rhoV)<1)
&& (is.null(p$rhoU) || abs(p$rhoU)<1)){
e <- y - X %*% p$gamma/p$nu
if (isSpY){
e <- e - p$rhoY * W_y %*% y
}
tryCatch({
sigmaV <- 1/(p$nu*sqrt(1+p$lambda^2))
sigmaU <- p$lambda * sigmaV
mu <- 0
if(isTN){
mu <- p$mu
}
vMu = rep(mu, n)
r <- 0
if (isSpV || isSpU){
SpV <- I
SpU <- I
if (isSpV)
SpV <- solve(I-p$rhoV*W_v)
if (isSpU)
SpU <- solve(I-p$rhoU*W_u)
mSigmaV = sigmaV^2*SpV%*%t(SpV)
mSigmaU = sigmaU^2*SpU%*%t(SpU)
mSigma = mSigmaV + mSigmaU
imSigma <- solve(mSigma)
mDelta <- mSigmaU%*%imSigma%*%mSigmaV
mGamma <- -costf*mSigmaU%*%imSigma
rownames(mDelta) <- colnames(mDelta)
f1 <- pmvnorm(lower=-Inf,upper = 0, mean=-vMu, sigma=mSigmaU, algorithm=alg)
r <- -log(f1)
r <- r + dmvnorm(x=as.vector(e),mean=-costf*vMu, sigma=mSigma,log=TRUE)
f2 <- pmvnorm(lower=-Inf,upper = as.vector(mGamma%*%(e+costf*vMu)), mean=-vMu, sigma=mDelta, algorithm=alg)
r <- r + log(f2)
}else{
r <- r - n*log(2*pi)/2
r <- r - n*log(pnorm(mu/sigmaU))
r <- r + n * log(p$nu)
r <- r - 0.5 * p$nu^2*sum((e + costf*vMu)^2)
r <- r + sum(log(pnorm(mu*p$nu/p$lambda - costf*p$lambda*e*p$nu)))
}
if (isSpY)
r <- r + determinant(I-p$rhoY*W_y,logarithm = TRUE)$modulus
ret <- r
}, error = function(e){
logging(e$message, level="debug")
})
}else{
if(!quiet) logging("Parameters are out of space")
}
if(!quiet) logging(paste("LogL =",ret))
if (ret == -Inf || is.nan(ret)) ret<- -Infin #For fake finite differences
logging(paste("Elapsed time: ",Sys.time()-st,"\n"))
return(as.numeric(ret))
}
funcLogL.direct <- function(parameters){
y <- envirGet("y")
X <- envirGet("X")
k <- ncol(X)
W_y <- envirGet("W_y")
W_v <- envirGet("W_v")
W_u <- envirGet("W_u")
inefficiency <- envirGet("inefficiency")
isSpY <- !is.null(W_y)
isSpV <- !is.null(W_v)
isSpU <- !is.null(W_u)
isTN <- (inefficiency == "truncated")
p <- paramsFromVector(parameters, k, isSpY, isSpV, isSpU, isTN, olsen = F)
polsen <- olsenReparam(p)
return(funcLogL(paramsToVector(polsen), quiet = F))
}
funcGradient<-function(parameters){
logging("Evaluating 'funcGradient' for parameters:",parameters)
y = envirGet("y")
X = envirGet("X")
k <- ncol(X)
W_y <- envirGet("W_y")
W_v <- envirGet("W_v")
W_u <- envirGet("W_u")
inefficiency <- envirGet("inefficiency")
costFrontier <- envirGet("costFrontier")
isSpY <- !is.null(W_y)
isSpV <- !is.null(W_v)
isSpU <- !is.null(W_u)
isTN <- (inefficiency == "truncated")
costf <- ifelse(costFrontier, -1, 1)
p <- paramsFromVector(parameters, k, isSpY, isSpV, isSpU, isTN)
yst <- y
if(isSpY){
yst <- y-p$rhoY*W_y%*%y
}
mu <- 0
if (isTN){
mu <- p$mu
}
omega <- p$nu * yst - X%*%p$gamma
a <- -costf*omega * p$lambda + mu*p$nu/p$lambda
delta <- exp(dnorm(a, log=T)-pnorm(a, log.p=T))
#delta <- dnorm(a)/pnorm(a)
N <- length(y)
sigmaU <- p$lambda/(p$nu*sqrt(1+p$lambda^2))
a2 <- mu/sigmaU
delta2 <- (dnorm(a2)/pnorm(a2))
dLdGamma <- t(omega+costf*mu*p$nu)%*%X + costf*t(delta)%*%X * p$lambda
dLdNu <- -t(omega+costf*mu*p$nu)%*%(y+mu) - t(delta)%*%(costf*y*p$lambda-mu/p$lambda) + N * 1/p$nu - N*delta2*a2/p$nu
dLdLambda <- -t(delta)%*%(costf*omega+mu*p$nu/p$lambda^2)+N*delta2*mu*p$nu/(p$lambda^2*sqrt(1+p$lambda^2))
grad = c(dLdGamma)
names(grad) = paste("dGamma", seq(length(dLdGamma)), sep = "")
if(isSpY){
dLdRhoY <- NaN
try({
mat = solve(diag(N)-p$rhoY*W_y) %*% (-W_y)
dLdRhoY = p$nu*(t(omega)+costf*p$lambda*t(delta))%*%W_y%*%y+sum(diag(mat))
})
ns <- names(grad)
grad <- c(grad, ifelse(is.nan(dLdRhoY),0,dLdRhoY))
names(grad) <- c(ns, "dRhoY")
}
ns <- names(grad)
grad = c(grad, dLdNu,dLdLambda)
names(grad) = c(ns, "dNu", "dLambda")
if(isTN){
dLdMu <- -N*delta2/sigmaU - p$nu*sum(omega+costf*mu*p$nu)+sum(delta)*p$nu/p$lambda
#dLdMu = -(mu - 1)
ns <- names(grad)
grad <- c(grad, dLdMu)
names(grad) <- c(ns, "dMu")
}
logging("Gradient",grad)
if (any(is.na(grad))) {
grad <- NULL
}
return(grad)
}
#' @title Calculation of the log likelihood function for the spatial stochastic frontier model
#'
#' @description
#' \code{logLikelihood} returns a value of the log likelihood function
#' for the spatial stochastic frontier model
#'
#' @details
#' This function is exported from the package for testing and presentation purposes
#' A list of arguments of the function exactly matches the corresponding list of the \code{\link{spfrontier}} function
#'
#'
#' @param formula an object of class "\code{\link{formula}}"
#' @param data data frame, containing the variables in the model
#' @param W_y a spatial weight matrix for spatial lag of the dependent variable
#' @param W_v a spatial weight matrix for spatial lag of the symmetric error term
#' @param W_u a spatial weight matrix for spatial lag of the inefficiency error term
#' @param values a vector of log likelihood function parameters
#' @param logging an optional level of logging. Possible values are 'quiet','warn','info','debug'.
#' By default set to quiet.
#' @param inefficiency sets the distribution for inefficiency error component. Possible values are 'half-normal' (for half-normal distribution) and 'truncated' (for truncated normal distribution).
#' By default set to 'half-normal'.
#' @param costFrontier is designed for selection of cost or production frontier
#'
#' @export
logLikelihood <- function(formula, data,
W_y = NULL, W_v = NULL,W_u = NULL,
inefficiency = "half-normal",
values,
logging = c("quiet", "info", "debug"), costFrontier=F){
logging <- match.arg(logging)
con <- list(grid.beta0 = 1, grid.sigmaV = 1, grid.sigmaU = 1, grid.rhoY = 1, grid.rhoU = 10, grid.rhoV = 10, grid.mu = 1)
initEnvir(W_y=W_y, W_v=W_v,W_u=W_u,inefficiency=inefficiency, logging=logging)
logging("Estimator started", level="info")
mf <- model.frame(formula, data)
y <- as.matrix(model.response(mf))
X <- as.matrix(mf[-1])
tm <- attr(mf, "terms")
intercept <- attr(tm, "intercept") == 1
if (intercept){
X <- cbind(Intercept=1L,X)
}
k <- ncol(X)
n <- length(y)
envirAssign("X", X)
envirAssign("y", y)
envirAssign("costFrontier", costFrontier)
isSpY <- !is.null(W_y)
isSpV <- !is.null(W_v)
isSpU <- !is.null(W_u)
isTN <- (inefficiency == "truncated")
res <- funcLogL(paramsToVector(olsenReparam(paramsFromVector(values, k, isSpY, isSpV, isSpU, isTN, olsen=F))))
names(res)<-"Log-likelihood function"
return(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.