Nothing
ml.truncreg <- function(param, X, y, gradient = FALSE, hessian = FALSE, fit = FALSE, point, direction, scaled){
beta <- param[1:ncol(X)]
sigma <- param[length(param)]
bX <- as.numeric(crossprod(beta, t(X)))
sgn <- ifelse(direction == "left", +1, -1)
if (!scaled){
resid <- (y - bX)
trunc <- (bX - point)
# Update of the mills function, use logs to avoid Inf
#mills <- dnorm(sgn * trunc / sigma) / pnorm(sgn * trunc / sigma)
mills <- exp(dnorm(sgn * trunc / sigma, log = TRUE) - pnorm(sgn * trunc / sigma, log.p = TRUE))
lnL <- - log(sigma) + dnorm(resid / sigma, log = TRUE) - pnorm(sgn * trunc / sigma, log.p = TRUE)
if (gradient){
gbX <- resid / sigma ^ 2 - sgn * mills / sigma
gsigma <- - 1 / sigma + resid ^ 2 / sigma ^ 3 + sgn * mills * trunc / sigma ^ 2
gradi <- cbind(gbX * X, as.numeric(gsigma))
attr(lnL, "gradient") <- gradi
}
if (fit) attr(lnL, "fit") <- bX
if (hessian){
bb <- -1 / sigma ^ 2 + mills * (sgn * trunc / sigma + mills) / sigma ^ 2
ss <- 1 / sigma ^ 2 - 3 * resid ^ 2 / sigma ^ 4 - 2 * mills * sgn * trunc / sigma ^ 3 +
mills * (sgn * trunc / sigma + mills) * trunc / sigma ^ 3
ss <- 1 / sigma ^ 2 - 3 * resid ^ 2 / sigma ^ 4 - 2 * mills * sgn * trunc / sigma ^ 3 +
mills * (sgn * trunc / sigma + mills) * trunc^2 / sigma ^ 4
bs <- - 2 * resid / sigma ^ 3 + sgn * mills / sigma ^ 2 -
mills * (mills + sgn * trunc / sigma) * trunc / sigma ^ 3
bb <- crossprod(bb * X, X)
bs <- apply(bs * X, 2, sum)
ss <- sum(ss)
h <- rbind(cbind(bb, bs), c(bs, ss))
attr(lnL,"hessian") <- h
}
}
else{
lnL <- - log(sigma) + dnorm(y / sigma - bX, log = TRUE) - pnorm(sgn * (bX - point / sigma), log.p = TRUE)
# Update of the mills function, use logs to avoid Inf (YC 2015/12/11)
# mills <- dnorm(sgn * (bX - point / sigma)) / pnorm(sgn * (bX - point / sigma))
mills <- exp(dnorm(sgn * (bX - point / sigma), log = TRUE) - pnorm(sgn * (bX - point / sigma), log.p = TRUE))
if (gradient){
gbX <- (y / sigma - bX) - mills * sgn
gsigma <- - 1 / sigma + (y / sigma - bX) * y / sigma ^ 2 - sgn * mills * point / sigma ^ 2
gradi <- cbind(gbX * X, as.numeric(gsigma))
attr(lnL, "gradient") <- gradi
}
if (fit) attr(lnL, "fit") <- bX * sigma
if(hessian){
bb <- -1 + mills * (mills + sgn * (bX - point / sigma))
bs <- - y / sigma ^ 2 + (mills + sgn * (bX - point / sigma)) * mills * point / sigma ^ 2
ss <- 1 / sigma ^ 2 - 3 * y ^ 2 / sigma ^ 4 + 2 * bX * y / sigma ^ 3 +
mills * (mills + sgn * (bX - point / sigma)) * point ^ 2 / sigma ^ 4 +
2 * sgn * mills * point / sigma ^ 3
bb <- crossprod(bb * X, X)
bs <- apply(bs * X, 2, sum)
ss <- sum(ss)
h <- rbind(cbind(bb, bs), c(bs, ss))
attr(lnL,"hessian") <- h
}
}
lnL
}
truncreg <- function(formula, data, subset, weights, na.action, point = 0, direction = "left",
model = TRUE, y = FALSE, x = FALSE, scaled = FALSE, ...){
formula.type <- TRUE
if (class(formula[[3]]) == "name"){
X <- try(eval(formula[[3]],sys.frame(which=-3)),silent = TRUE)
if (class(X) == "matrix") formula.type <- FALSE else formula.type <- TRUE
}
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
if (formula.type){
m <- match(c("formula", "data", "subset", "na.action","weights"),
names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
X <- model.matrix(formula, data = mf)
Y <- model.response(mf)
mt <- attr(mf, "terms")
}
else{
Y <- eval(formula[[2]], sys.frame(which=-3))
mt <- terms(formula)
}
## process options
direction <- match.arg(direction, c("left", "right"))
point <- rep(point, length.out = length(Y))
result <- truncreg.fit(X, Y, point, direction, scaled, ...)
result$call <- cl
result$terms <- mt
if(model) result$model <- mf
if(y) result$y <- Y
if(x) result$x <- X
result
}
truncreg.fit <- function(X, y, point, direction, scaled, ...){
## check input
if(direction == "left" & any(y < point)) stop("response not truncated, contains observations < 'point'")
if(direction == "right" & any(y > point)) stop("response not truncated, contains observations > 'point'")
dots <- list(...)
if (is.null(dots$method)) method <- "bfgs" else method <- dots$method
if (is.null(dots$iterlim)) iterlim <- 50 else iterlim <- dots$iterlim
if (is.null(dots$print.level)) print.level <- 0 else print.level <- dots$print.level
oldoptions <- options(warn = - 1)
on.exit(options(oldoptions))
start.time <- proc.time()
f <- function(param) ml.truncreg(param, X = X, y = y,
gradient = TRUE, hessian = TRUE,
fit = FALSE, point = point,
direction = direction, scaled = scaled)
linmod <- lm(y ~ X - 1)
start <- c(coef(linmod), summary(linmod)$sigma)
if (scaled) start[1:ncol(X)] <- start[1:ncol(X)] / start[ncol(X) + 1]
if (FALSE){
f0 <- ml.truncreg(start, X = X, y = y,
gradient = TRUE, hessian = FALSE,
fit = FALSE, point = point, direction = direction, scaled = scaled)
agrad <- apply(attr(f0, "gradient"), 2, sum)
ostart <- start
of <- sum(f(start))
ngrad <- c()
eps <- 1E-5
for (i in 1:length(start)){
start <- ostart
start[i] <- start[i] + eps
ngrad <- c(ngrad, (sum(f(start)) - of) / eps)
}
start <- ostart
print(cbind(agrad, ngrad))
stop()
}
maxl <- maxLik(f, start = start, method = method,
iterlim = iterlim, print.level = print.level)
coefficients <- maxl$estimate
if (scaled){
coefficients[1:ncol(X)] <- coefficients[1:ncol(X)] * coefficients[ncol(X) + 1]
f <- function(param) ml.truncreg(param, X = X, y = y,
gradient = TRUE, hessian = TRUE,
fit = FALSE, point = point,
direction = direction, scaled = FALSE)
maxl <- maxLik(f, start = coefficients, method = method,
iterlim = 0, print.level = print.level)
}
coefficients <- maxl$estimate
grad.conv <- maxl$gradient
vcov <- - solve(maxl$hessian)
fit <- attr(ml.truncreg(coefficients, X = X, y = y,
gradient = TRUE, hessian = TRUE,
fit = TRUE, point = point,
direction = direction, scaled = scaled), "fit")
names(fit) <- rownames(X)
logLik <- maxl$maximum
attr(logLik,"df") <- length(coefficients)
hessian <- maxl$hessian
convergence.OK <- maxl$code<=2
elaps.time <- proc.time() - start.time
nb.iter <- maxl$iterations
eps <- maxl$gradient%*%solve(-maxl$hessian)%*%maxl$gradient
est.stat <- list(elaps.time = elaps.time,
nb.iter = nb.iter,
eps = eps,
method = maxl$type,
message = maxl$message
)
class(est.stat) <- "est.stat"
coef.names <- c(colnames(X),"sigma")
names(coefficients) <- rownames(vcov) <- colnames(vcov) <- coef.names
result <- list(coefficients = coefficients,
vcov = vcov,
fitted.values = fit,
logLik = logLik,
gradient = grad.conv,
gradientObs = maxl$gradientObs,
nobs = length(y),
call = NULL,
terms = NULL,
model = NULL,
y = NULL,
x = NULL,
point = if(isTRUE(all.equal(rep(point[1], length(point)), point))) point[1] else point,
direction = direction,
est.stat = est.stat
)
class(result) <- c("truncreg", "maxLik")
result
}
fitted.truncreg <- function(object, ...){
object$fitted.values
}
residuals.truncreg <- function(object, ...){
model.frame(object)[[1]]-fitted(object)
}
coef.truncreg <- function(object, ...){
object$coefficients
}
vcov.truncreg <- function(object, ...){
object$vcov
}
logLik.truncreg <- function(object, ...)
structure(object$logLik, df = length(object$coefficients), nobs = object$nobs, class = "logLik")
print.truncreg <- function (x, digits = max(3, getOption("digits") - 2), width = getOption("width"), ...){
cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
if (length(coef(x))) {
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2,
quote = FALSE)
}
else cat("No coefficients\n")
cat("\n")
invisible(x)
}
summary.truncreg <- function (object,...){
b <- coef(object)
std.err <- sqrt(diag(vcov(object)))
z <- b/std.err
p <- 2*(1-pnorm(abs(z)))
CoefTable <- cbind(b,std.err,z,p)
colnames(CoefTable) <- c("Estimate","Std. Error","t-value","Pr(>|t|)")
object$coefficients <- CoefTable
class(object) <- c("summary.truncreg","truncreg")
return(object)
}
print.summary.truncreg <- function(x,digits= max(3, getOption("digits") - 2),width=getOption("width"),...){
cat("\nCall:\n")
print(x$call)
if (!is.null(x$est.stat)){
cat("\n")
print(x$est.stat)
}
cat("\n")
cat("\nCoefficients :\n")
printCoefmat(x$coefficients,digits=digits)
cat("\n")
df <- attr(x$logLik,"df")
cat(paste("Log-Likelihood: ",
signif(x$logLik,digits),
" on ",df," Df\n",sep=""))
invisible(x)
}
model.frame.truncreg <- function(formula, ...) {
if(!is.null(formula$model)) return(formula$model)
NextMethod()
}
model.matrix.truncreg <- function(object, ...)
if(!is.null(object$x)) object$x else model.matrix(object$terms, model.frame(object), ...)
predict.truncreg <- function(object, newdata = NULL, na.action = na.pass, ...)
{
if(missing(newdata)) {
rval <- object$fitted.values
} else {
mt <- delete.response(object$terms)
X <- model.matrix(mt, model.frame(mt, newdata, na.action = na.action))
rval <- drop(X %*% head(object$coefficients, -1))
}
return(rval)
}
print.est.stat <- function(x, ...){
et <- x$elaps.time[3]
i <- x$nb.iter[1]
halton <- x$halton
method <- x$method
if (!is.null(x$type) && x$type != "simple"){
R <- x$nb.draws
cat(paste("Simulated maximum likelihood with", R, "draws\n"))
}
s <- round(et,0)
h <- s%/%3600
s <- s-3600*h
m <- s%/%60
s <- s-60*m
cat(paste(method, "method\n"))
tstr <- paste(h, "h:", m, "m:", s, "s", sep="")
cat(paste(i,"iterations,",tstr,"\n"))
if (!is.null(halton)) cat("Halton's sequences used\n")
if (!is.null(x$eps)) cat(paste("g'(-H)^-1g =", sprintf("%5.3G", as.numeric(x$eps)),"\n"))
if (is.numeric(x$code)){
msg <- switch(x$code,
"1" = "gradient close to zero",
"2" = "successive function values within tolerance limits",
"3" = "last step couldn't find higher value",
"4" = "iteration limit exceeded"
)
cat(paste(msg, "\n"))
}
else cat(paste(x$code, "\n"))
}
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.