R/ols.interval.R

In desk: Didactic Econometrics Starter Kit

ols.interval = function(mod,
                        data = list(),
                        type = c("confidence", "prediction", "acceptance"),
                        which.coef = "all",
                        sig.level = 0.05,
                        q = 0,
                        dir = c("both", "left", "right"),
                        xnew,
                        details = FALSE){

  if (inherits(mod, "formula")) { # Wenn Formel übergeben ...
    mod = ols(mod, data = data)
  }

  if(!(inherits(mod,"desk") | (attr(mod, "type") == "ols"))){
    stop("Model is not compatible. Please use object generated by ols() from package desk.", call. = F)
  }

  type = match.arg(type)
  dir = match.arg(dir)

  # Function to convert percentage
  format.perc = function(probs, digits){
    paste(format(100 * probs, trim = TRUE, scientific = FALSE, digits = digits),
          "%", sep = "")
    }

  co = mod$coefficients # Uebergebene Koeffizienten
  cn = names(co) # Extrahiere Namen

  # Wenn keine Parameter angegeben, wähle alle...
  if(length(which.coef) == 1){if(which.coef == "all"){which.coef = cn}} else {
    if(is.numeric(which.coef)) {which.coef = cn[which.coef]}
  }

  # Case: direction = "both"
  a = c(sig.level/2, 1 - sig.level/2)
  tval = qt(a, mod$df)
  header = c("center", format.perc(a, 3))
  out = array(NA, dim = c(length(which.coef), 3L))

  ##############################
  # Confidence intervals #######
  ##############################
  if (type == "confidence"){
    if (dir != "both") stop("One-sided confidence intervals are not supported.", call. = F)
    if (missing(xnew)){ # CI for model's parameter
      #message("No xnew was specified. Will calculate confidence intervals for the model's parameters ...", call. = F)
      se = sqrt(diag(mod$vcov)[which.coef])
      out = cbind(co[which.coef], co[which.coef] + se %o% tval)
      colnames(out) = header
      c.type = "ci.beta"
    } else { # CI for model's expected y-value
      if (inherits(xnew, "numeric")){xnew = t(xnew)}
      if(mod$ncoef - sum(mod$has.const) == dim(xnew)[2]){
      if(ols.has.const(mod)){
        xnew = cbind(1,xnew) # Füge Vektor von Einsen ein
        }
      ynew = as.vector(xnew %*% co) # Fitted values
      se = sqrt(diag(xnew %*% mod$vcov %*% t(xnew))) # Vektor von quadratischen Formen
      out = cbind(ynew, ynew + se %o% tval)
      dimnames(out) = list(1:dim(xnew)[1], header)
      c.type = "ci.y"
      } else {
        stop("xnew must be (T x K) where K is the number of predictors", call. = F)
      }
    }
  }

  ##############################
  # Prediction intervals #######
  ##############################
  if (type == "prediction"){
    if (dir != "both") stop("One-sided prediction intervals are not supported.", call. = F)
    if (missing(xnew)){
      # warning("No xnew was specified. Will use given observations for prediction.", call. = F)
      xnew = as.matrix(mod$data[,-1])
    }
    if (inherits(xnew, "numeric")){xnew = t(xnew)}
    if(mod$ncoef - sum(mod$has.const) == dim(xnew)[2]){
      if(ols.has.const(mod)){
        xnew = cbind(1,xnew) # Füge Vektor von Einsen ein
      }
      ynew = as.vector(xnew %*% co) # Fitted values
      se = sqrt(diag(xnew %*% mod$vcov %*% t(xnew)) + mod$sig.squ) # Vektor von quadratischen Formen
      out = cbind(ynew, ynew + se %o% tval)
      dimnames(out) = list(1:dim(xnew)[1], header)
      c.type = "pi"
    } else {
      stop("xnew must be (T x K) where K is the number of predictors", call. = F)
    }
  }

  ##############################
  # Acceptance intervals #######
  ##############################
  if (type == "acceptance"){
    se = sqrt(diag(mod$vcov)[which.coef])
    switch(dir,
           both = {
             out = cbind(q, q + se %o% tval)
           },
           right = {
             tval = qt(1-sig.level, mod$df)
             header = c("lower", format.perc(1-sig.level, 3))
             out = cbind(-Inf, q + se %o% tval)
           },
           left = {
             tval = qt(sig.level, mod$df)
             header = c(format.perc(sig.level, 3), "upper")
             out = cbind(q - se %o% tval, Inf)
           }
    )
    colnames(out) = header
    c.type = "ai"
  }

  out = list(results = out)
  out$std.err = se
  out$t.value = tval

  attr(out, "title") = NULL
  attr(out, "details") = if (details) {T} else {F}
  attr(out, "type") = "int"
  attr(out, "c.type") = c.type
  class(out) = c("desk")

return(out)
}