Nothing
###########################################################
## This file contains the source of the logitx()
## functions.
##
## Created 8 December 2020, Oslo.
##
## Contents:
##
## logitxSim()
## logit()
## logitx()
## coef.logitx #extraction functions
## fitted.logitx #(all are S3 methods)
## gets.logitx
## logLik.logitx
## plot.logitx
## WIP: predict.logitx
## print.logitx
## summary.logitx
## toLatex.logitx
## vcov.logitx
##
## dlogitxSim() #create alias for logitxSim
## dlogitx() #create alias for logitx
##
###########################################################
###########################################################
## The function logitxSim() simulates from an
## autoregressive logit model with covariates.
##
## Function arguments:
##
## n: integer, the number of observations to
## simulate from
## intercept: the value of the intercept in the logit
## specification
## ar: numeric vector w/value(s) of the AR-parameters
## xreg: numeric vector, e.g. theta1*x1+theta2*x2
## verbose: logical, whether to return vI only or all
## the info from the simulations
## as.zoo: logical, whether the returned object should
## be of class 'zoo' or not
##
## The conditional probability of yt = 1:
##
## pi1t = 1/( 1 + exp(-ht) )
##
## The logit specification:
##
## ht = intercept + ar-terms + covariates
## = intercept + theta1*y_1 + ... + thetap*y_p + covariates
##
###########################################################
logitxSim <- function(n, intercept=0, ar=NULL, xreg=NULL,
verbose=FALSE, as.zoo=TRUE)
{
##orders:
p.order <- ifelse(is.null(ar), 0, length(ar))
#q.order <- length(ma) (for the future)
maxpq <- max(1, p.order)
nplusMaxpq <- n + maxpq
##xreg:
if( is.null(xreg) ){
xreg <- rep(0, n)
xregMean <- 0
}else{
xreg <- as.vector(xreg)
xregMean <- mean(xreg)
}
##p.order = 0 (recursion not necessary):
if( p.order == 0 ){
ht <- intercept + xreg
pi1 <- 1/(1 + exp(-ht))
vItmp <- runif(n)-1
vI <- as.numeric( vItmp + pi1 > 0)
}
##if p.order > 0 (recursion necessary):
if( p.order > 0 ){
##ht:
htmean <- (intercept + xregMean)/(1 - sum(ar))
ht <- rep(htmean, n)
#pi1, vI:
vI <- rep(0, n)
pi1 <- 1/(1 + exp(-ht))
vItmp <- runif(n)-1
##recursion:
xregsum <- intercept + xreg
arsum <- 0
for(i in c(maxpq+1):n){
arsum <- sum( ar*vI[ c(i-1):c(i-p.order) ] )
ht[i] <- xregsum[i] + arsum
pi1[i] <- 1/(1+exp(-ht[i]))
vI[i] <- as.numeric( vItmp[i] + pi1[i] > 0)
}
} #end if p.order > 0
## result:
if( verbose ){
result <- cbind(vI, pi1, ht)
colnames(result) <- c("I", "pi1", "h")
}else{
result <- vI
}
if( as.zoo ){
result <- as.zoo(result)
}
return(result)
} #end logitxSim
###########################################################
## The function logit() estimates a logit model
##
## Function arguments:
##
## y: vector with the binary data
## x: matrix of covariates
## initial.values: initial values of the parameters
## lower: lower bounds of the permissible parameter values
## upper: upper bounds of the permissible parameter values
## method: estimation method, 1=only estimates,
## 2=ordinary vcov, 3=robust vcov
## lag.length: integer controlling the bandwidth when method=3
## control: options passed on to nlminb()
## eps.tol: pi0 is computed as pmax(pi0, eps.tol) to
## ensure log(pi0) is not too small
## solve.tol: tolerance for singularity when inverting
## a matrix with solve()
##
## The conditional probability of yt = 1:
##
## pi1t = 1/( 1 + exp(-ht) )
##
## The logit specification:
##
## ht = intercept + ar-terms + covariates
## = intercept + theta1*y_1 + ... + thetap*y_p + covariates
##
############################################################
logit <- function(y, x, initial.values=NULL, lower=-Inf, upper=Inf,
method=2, lag.length=NULL, control=list(), eps.tol=.Machine$double.eps,
solve.tol=.Machine$double.eps)
{
##initiate:
##=========
out <- list()
out$n <- length(y)
if( is.null(x) ){
out$k <- 0
}else {
out$k <- NCOL(x)
}
out$df <- out$n - out$k
##initial values:
if( is.null(initial.values) && out$k>0 ){
out$initial.values <- rep(0.1, out$k)
}else{
out$initial.values <- initial.values
}
##bounds, tolerance:
out$lower <- lower
out$upper <- upper
out$eps.tol <- eps.tol
##create log-likelihood function:
##===============================
logitObjective <- function(pars, aux, minimise=TRUE){
##check parameters:
parsOK <- TRUE
if( any(is.na(pars)) || any(pars<aux$lower) || any(pars>aux$upper) ){
parsOK <- FALSE
}
##compute log-likelihood:
if( parsOK ){
ht <- as.vector( aux$x %*% pars )
pi1that <- 1/(1+exp(-ht))
pi0that <- pmax(1-pi1that, aux$eps.tol)
logl <- sum( aux$y*log(pi1that) + aux$oneminy*log(pi0that) )
if(is.nan(logl) || is.na(logl) || abs(logl) == Inf){
logl <- aux$logl.penalty
}
}else{
logl <- aux$logl.penalty
}
if(minimise){ logl <- -logl }
##result:
return(logl)
}
##create aux:
##===========
aux <- out
aux$y <- y
aux$x <- x
aux$oneminy <- 1-y
aux$penalty.value <- out$n*log(0.5)
if( out$k>0 ){
aux$penalty.value <- logitObjective(out$initial.values, aux)
}
##estimate:
##=========
##w/regressors:
if( out$k>0 ){
result <- nlminb(out$initial.values, logitObjective, aux=aux,
lower=lower, upper=upper, control=control)
names(result)[1] <- "coefficients"
result$objective <- -result$objective #due to minimisation
names(result)[2] <- "logl"
ht <- as.vector(x %*% result$coefficients)
result$fit <- 1/(1+exp(-ht)) #pi1hat
}
##if there are no regressors:
if( out$k==0 ){
result <- list()
result$coefficients <- NULL
result$logl <- out$n*log(0.5)
result$convergence <- 0
result$iterations <- 0
result$evaluations <- c(1,0)
names(result$evaluations) <- c("function", "gradient")
#result$message <- "convergence (no iterations)"
result$fit <- rep(0.5, out$n)
}
##merge out with result:
result <- c(out, result)
##vcov:
##=====
##compute the hessian:
if( method>1 && result$k>0 ){
xadj <- (result$fit^2 - result$fit)*x
hessian <- crossprod(xadj, x)
hessianinv <- solve(hessian, tol=solve.tol)
}
##ordinary vcov:
if( method==2 && result$k>0 ){
result$vcov <- -hessianinv
}
##robust coefficient covariance:
if( method==3 && result$k>0 ){
##lag length (bandwidth):
if( is.null(lag.length) ){
##EViews, see Wooldridge (2009, p. 430):
iL <- as.integer(4*(result$n/100)^(2/9))
}else{
iL <- as.integer(lag.length)
}
result$lag.length <- iL
vW <- 1 - 1:iL/(iL+1) #weights
##compute the "meat":
mS <- cbind((y - result$fit) * x)
mSigmahat <- crossprod(mS)
for(l in 1:iL){
mS0 <- cbind(mS[-c(1:l),])
mSj <- cbind(mS[-c(c(result$n-l+1):result$n),])
mSigmahat <-
mSigmahat + vW[l] * (crossprod(mS0,mSj) + crossprod(mSj,mS0))
}
##compute vcov:
result$vcov <- hessianinv %*% mSigmahat %*% hessianinv
}
##return result:
##==============
return(result)
} #close logit function
############################################################
## estimate model of class "logitx":
logitx <- function(y, intercept=TRUE, ar=NULL, ewma=NULL,
xreg=NULL, vcov.type=c("ordinary", "robust"), lag.length=NULL,
initial.values=NULL, lower=-Inf, upper=Inf, control=list(),
eps.tol=.Machine$double.eps, solve.tol=.Machine$double.eps,
singular.ok=TRUE, plot=NULL)
{
##auxiliary list:
aux <- list()
aux$call <- sys.call()
aux$date <- date()
aux$version <- paste0("gets ", packageVersion("gets"), " under ",
version$version.string)
aux$control <- control
##modify ewma argument:
if( !is.null(ewma) && !is.list(ewma) && is.vector(ewma) ){
ewma <- list(length=ewma)
}
##variables:
tmp <- regressorsMean(y, mc=intercept, ar=ar, ewma=ewma, mxreg=xreg,
prefix=character(0), return.regressand=TRUE, return.as.zoo=TRUE,
na.trim=TRUE, na.omit=FALSE)
whereMconst <- which( colnames(tmp)=="const" )
if( length(whereMconst)>0 ){
colnames(tmp)[ whereMconst ] <- "intercept"
}
tmpnames <- colnames(tmp)
##regressand:
aux$y <- coredata(tmp[,1])
aux$y.name <- colnames(tmp)[1]
aux$y.index <- index(tmp)
##regressors:
if( NCOL(tmp)==1 ){
tmp <- NULL
tmpnames <- NULL
}
if( NCOL(tmp)>1 ){
tmp <- cbind(coredata(tmp[,-1]))
tmpnames <- tmpnames[-1]
}
if( singular.ok && NCOL(tmp)> 1 ){
tmp <- dropvar(tmp, tol=1e-07, LAPACK=FALSE, silent=TRUE)
droppedVars <- setdiff(tmpnames, colnames(tmp))
tmpnames <- setdiff(tmpnames, droppedVars)
if( length(droppedVars)>0 ){
warning("\n", "regressor-matrix singular, so dropping: ",
droppedVars, "\n")
}
}
aux$mX <- tmp
aux$mXnames <- tmpnames
if( !is.null(aux$mX) ){
colnames(aux$mX) <- NULL
}
##determine estimation method/vcov type:
types <- c("ordinary", "robust")
whichType <- charmatch(vcov.type[1], types)
aux$logit.method <- whichType + 1
##estimate:
result <- logit(aux$y, aux$mX, initial.values=initial.values,
lower=lower, upper=upper, method=aux$logit.method,
lag.length=lag.length, control=control, eps.tol= eps.tol,
solve.tol=solve.tol)
names(result$coefficients) <- aux$mXnames
colnames(result$vcov) <- aux$mXnames
rownames(result$vcov) <- aux$mXnames
result$fit <- zoo(result$fit, order.by=aux$y.index)
result <- c(aux, result)
##plot:
if( is.null(plot) ){
plot <- getOption("plot")
if( is.null(plot) ){ plot <- FALSE }
}
if(plot){ plot.logitx(result) }
##return result:
class(result) <- c("logitx", "dlogitx")
return(result)
} #close logitx()
############################################################
## extract coefficients
coef.logitx <- function(object, ...)
{
object$coefficients
} #close coef.logitx
############################################################
## extract fitted probabilities
fitted.logitx <- function(object, zero.prob=FALSE, ...)
{
if( zero.prob ){
result <- 1-object$fit
}else{
result <- object$fit
}
return(result)
} #close fitted.logitx
############################################################
## do gets on a logitx object
gets.logitx <- function(x, t.pval=0.05, wald.pval=t.pval,
do.pet=TRUE, user.diagnostics=NULL, keep=NULL,
include.gum=FALSE, include.1cut=TRUE, include.empty=FALSE,
max.paths=NULL, turbo=TRUE, print.searchinfo=TRUE, plot=NULL,
alarm=FALSE, ...)
{
##logit() arguments:
##------------------
callArgs <- names(x$call)
userEstArgs <- list()
userEstArgs$name <- "logit"
if( "lower" %in% callArgs ){ userEstArgs$lower <- x$lower }
if( "upper" %in% callArgs ){ userEstArgs$lower <- x$upper }
userEstArgs$method <- x$logit.method
if( "lag.length" %in% callArgs ){ userEstArgs$lag.length <- x$lag.length }
userEstArgs$control <- x$control
if( "eps.tol" %in% callArgs ){ userEstArgs$eps.tol <- eval(x$call$eps.tol) }
if( "solve.tol" %in% callArgs ){
userEstArgs$solve.tol <- eval(x$call$solve.tol)
}
##do gets:
##--------
vY <- x$y
mX <- x$mX
result1 <- getsFun(vY, mX, user.estimator=userEstArgs,
gum.result=x, t.pval=t.pval, wald.pval=wald.pval, do.pet=do.pet,
user.diagnostics=user.diagnostics, keep=keep,
include.gum=include.gum, include.1cut=include.1cut,
include.empty=include.empty, max.paths=max.paths, turbo=turbo,
print.searchinfo=print.searchinfo, alarm=alarm)
result1$call <- NULL
##create gum result:
##------------------
if( x$k>0 ){
gum.result <- matrix(0, x$k, 6)
colnames(gum.result) <- c("reg.no.", "keep", "coef", "std.error",
"t-stat", "p-value")
rownames(gum.result) <- x$mXnames
gum.result[,"reg.no."] <- 1:x$k
if( !is.null(keep) ){ gum.result[keep,"keep"] <- 1 }
gum.result[,"coef"] <- x$coefficients
gum.result[,"std.error"] <- sqrt(diag(x$vcov))
gum.result[,"t-stat"] <- gum.result[,"coef"]/gum.result[,"std.error"]
gum.result[,"p-value"] <-
pt(abs(gum.result[,"t-stat"]), df=x$df, lower.tail=FALSE)
result1$gum.result <- gum.result
}
##estimate specific:
##------------------
vY <- cbind(zoo(vY, order.by=x$y.index))
colnames(vY) <- x$y.name
xfinal <- result1$terminals[[ result1$best.terminal ]]
vcov.type <- c("none", "ordinary", "robust")[ x$logit.method ]
solve.tol <- ifelse( is.null(userEstArgs$solve.tol),
.Machine$double.eps, userEstArgs$solve.tol)
##empty final model:
if( length(xfinal)==0 ){
result2 <- logitx(vY, intercept=FALSE, vcov.type=vcov.type,
lag.length=x$lag.length, lower=x$lower, upper=x$upper,
control=x$control, eps.tol=x$eps.tol, solve.tol=solve.tol,
plot=plot)
}
##non-empty final model:
if( length(xfinal)>0 ){
mX <- cbind(mX[,xfinal])
colnames(mX) <- x$mXnames[ xfinal ]
mX <- zoo(mX, order.by=x$y.index)
result2 <- logitx(vY, intercept=FALSE, xreg=mX, vcov.type=vcov.type,
lag.length=x$lag.length, lower=x$lower, upper=x$upper,
control=x$control, eps.tol=x$eps.tol, solve.tol=solve.tol,
plot=plot)
}
##return result:
##--------------
result2$call <- NULL
result <- c(result1, result2)
class(result) <- c("logitx", "dlogitx")
return(result)
} #close gets.logitx()
############################################################
## extract log-likelihood
logLik.logitx <- function(object, ...)
{
result <- object$logl
attr(result, "df") <- length(object$coefficients)
attr(result, "nobs") <- object$n
class(result) <- "logLik"
return(result)
}
############################################################
## plot fitted probabilities
plot.logitx <- function(x, ...)
{
plot(x$fit, ylab="probability", xlab="", col="blue")
} #close plot.logitx
############################################################
## print estimation result of model of class "logitx":
print.logitx <- function(x, signif.stars=TRUE, ...)
{
##header:
##-------
##first part:
cat("\n")
cat("Date:", x$date, "\n")
cat("Dependent var.:", x$y.name, "\n")
cat("Method: Maximum Likelihood (logit) \n")
cat("Variance-Covariance:",
switch(x$logit.method, "1" = "Ordinary", "2" = "Ordinary",
"3" = "Robust, Newey and West (1987)"), "\n")
cat("No. of observations:", x$n, "\n")
##sample info:
isRegular <- is.regular(x$fit, strict=TRUE)
isCyclical <- frequency(x$fit) > 1
indexTrimmed <- index(x$fit)
if(isRegular && isCyclical){
cycleTrimmed <- cycle(x$fit)
startYear <- floor(as.numeric(indexTrimmed[1]))
startAsChar <- paste(startYear,
"(", cycleTrimmed[1], ")", sep="")
endYear <- floor(as.numeric(indexTrimmed[length(indexTrimmed)]))
endAsChar <- paste(endYear,
"(", cycleTrimmed[length(indexTrimmed)], ")", sep="")
}else{
startAsChar <- as.character(indexTrimmed[1])
endAsChar <- as.character(indexTrimmed[length(indexTrimmed)])
}
cat("Sample:", startAsChar, "to", endAsChar, "\n")
##if gets modelling:
##------------------
if( !is.null(x$gum.result) ){
cat("\n")
cat("Start model (GUM):\n")
cat("\n")
printCoefmat(x$gum.result, tst.ind=c(1,2), signif.stars=signif.stars)
##paths:
cat("\n")
cat("Paths searched: \n")
cat("\n")
if( is.null(x$paths) || length(x$paths)==0 ){
print(NULL)
}else{
for(i in 1:length(x$paths)){
cat("path",i,":",x$paths[[i]],"\n")
}
} #end if(is.null(x$paths))
##terminal models and results:
if( !is.null(x$terminals) && length(x$terminals)>0 ){
cat("\n")
cat("Terminal models: \n")
if(!is.null(x$terminals)){
cat("\n")
for(i in 1:length(x$terminals)){
cat("spec",i,":",x$terminals[[i]],"\n")
}
}
}
if( !is.null(x$terminals.results) ){
cat("\n")
printCoefmat(x$terminals.results, dig.tst=0, tst.ind=c(3,4),
signif.stars=FALSE)
}
} #close if gets modelling
##results matrix:
##---------------
if( length(x$coefficients)>0 ){
resultsmat <- matrix(NA, length(x$coefficients), 4)
colnames(resultsmat) <- c("coef", "std.error", "t-stat", "p-value")
rownames(resultsmat) <- names(x$coefficients)
resultsmat[,"coef"] <- x$coefficients
resultsmat[,"std.error"] <- sqrt(diag(x$vcov))
resultsmat[,"t-stat"] <- resultsmat[,1]/resultsmat[,2]
resultsmat[,"p-value"] <-
pt(abs(resultsmat[,3]), df=x$df, lower.tail=FALSE)
cat("\n")
if( is.null(x$gum.result) ){
cat("Estimation results:\n")
}else{
cat("Final model:\n")
}
cat("\n")
printCoefmat(resultsmat, signif.stars=signif.stars)
}else{
cat("\n")
cat(" The empty model\n")
}
##goodness-of-fit:
##----------------
gof <- matrix(NA, 1, 1)
rownames(gof) <- paste0("Log-lik.(n=", x$n, ")")
colnames(gof) <- ""
gof[1,1] <- x$logl
printCoefmat(gof, digits=6, signif.stars=signif.stars)
} #end print.logitx
############################################################
## summarise output
summary.logitx <- function(object, ...)
{
summary.default(object)
} #end summary.arx
############################################################
### LaTeX code (equation form)
toLatex.logitx <- function(object, digits=4, gof=TRUE,
nonumber=FALSE, nobs="T", ...) #new argument?: print.info=TRUE
{
##probability equation:
##---------------------
txtAddNonumber <- ifelse(nonumber, " \\nonumber ", "")
probtxt <- paste0(" Pr(", object$y.name,
"_t = 1| ...) &=& \\frac{1}{1 + \\exp(-\\widehat{h}_t)}",
txtAddNonumber, " \\\\[2mm]", " \n")
##logit equation:
##---------------
##y name:
hName <- paste0("\\widehat{h}_t")
##coefs, coef names, std.errors:
coefs <- coef.logitx(object)
coefsNames <- names(coefs)
whereIntercept <- which( coefsNames=="intercept" )
if( whereIntercept > 0 ){ coefsNames[ whereIntercept ] <- "" }
coefs <- as.numeric(coefs)
stderrs <- as.numeric(sqrt(diag(vcov.logitx(object))))
##equation (main part):
eqtxt <- NULL
if(length(coefs) > 0){
for(i in 1:length(coefs) ){
ifpluss <- ifelse(i==1, "", " + ")
eqtxt <- paste(eqtxt,
ifelse(coefs[i] < 0, " - ", ifpluss), "\\underset{(",
format(round(stderrs[i], digits=digits), nsmall=digits), ")}{",
format(round(abs(coefs[i]), digits=digits), nsmall=digits), "}",
coefsNames[i], sep="")
}
}
##equation (put parts together):
txtAddNonumber <- ifelse(nonumber, " \\nonumber ", "")
eqtxt <- paste0(" ", hName, " &=& ", eqtxt, "", txtAddNonumber,
" \\\\[2mm]", " \n")
##goodness of fit:
##----------------
goftxt <- paste(" && LogL=",
format(round(as.numeric(logLik.logitx(object)), digits=digits), nsmall=digits),
" \\qquad ", nobs, " = ", object$n, " \\nonumber \n", sep="")
##print code:
##-----------
# if( print.info ){
cat("% Date:", date(), "\n")
notetxt <- paste0("% LaTeX code generated in R ",
version$major, ".", version$minor, " by the gets package\n")
cat(notetxt)
cat("% Note: The {eqnarray} environment requires the {amsmath} package\n")
# }
cat("\\begin{eqnarray}\n")
cat(probtxt)
cat(eqtxt)
cat(goftxt)
cat("\\end{eqnarray}\n")
} #end toLatex.logitx
############################################################
## variance-covariance extraction function
vcov.logitx <- function(object, ...)
{
object$vcov
} #end vcov.logitx
############################################################
## create alias
dlogitxSim <- function(n, ...){ logitxSim(n, ...) }
############################################################
## create alias
dlogitx <- function(y, ...){ logitx(y, ...) }
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.