R/midas_sp_methods.R

In mpiktas/midasr: Mixed Data Sampling Regression

##' Fitted values for semi-parametric MIDAS regression model
##'
##' Returns the fitted values of a fitted semi-parametric MIDAS regression object
##' @param object a \code{\link{midas_r}} object
##' @param ... currently nothing
##' @return the vector of fitted values
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
##' @method fitted midas_sp
fitted.midas_sp <- function(object, ...) {
  fit <- midas_sp_fit(object, Xeval = NULL, Zeval = NULL)
  as.numeric(fit$fitted.values)
}

##' @export
##' @method print midas_sp
print.midas_sp <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  # Code adapted from dynln:::print.dynlm code
  model_string <- "\nSemi-parametric MIDAS regression model with \""
  cat(paste(model_string, class(x$lhs)[1],
    "\" data:\n",
    sep = ""
  ))
  cat(paste("Start = ", x$lhs_start,
    ", End = ", x$lhs_end,
    "\n",
    sep = ""
  ))
  cat(" model:", deparse(x$call$formula), "\n")
  print(coef(x), digits = digits, ...)
  cat("\n")
  cat("Function", x$argmap_opt$Ofunction, "was used for fitting\n")
  invisible(x)
}


##' @export
##' @importFrom stats deviance pt pnorm residuals printCoefmat cor
##' @method summary midas_sp
summary.midas_sp <- function(object, df = NULL, ...) {
  r <- as.vector(residuals(object))
  param <- coef(object)
  pnames <- names(param)
  n <- length(r)
  p <- length(param)
  rdf <- n - p

  resvar <- deviance(object) / rdf

  dg <- jacobian(object$rhs, param)
  V <- resvar * ginv(crossprod(dg) / nrow(dg)) / n

  colnames(V) <- pnames
  rownames(V) <- pnames

  se <- sqrt(diag(V))

  f <- as.vector(object$fitted.values)
  mss <- sum((f - mean(f))^2)
  rss <- sum(r^2)

  n <- length(r)
  p <- length(coef(object))
  rdf <- n - p
  df.int <- 0L

  r_squared <- R2.np(object$model[, 1], f)
  tval <- param / se

  # Code stolen from coeftest function
  if (is.null(df)) {
    df <- rdf
  }
  if (is.finite(df) && df > 0) {
    pval <- 2 * pt(abs(tval), df = df, lower.tail = FALSE)
    cnames <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
    mthd <- "t"
  }
  else {
    pval <- 2 * pnorm(abs(tval), lower.tail = FALSE)
    cnames <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
    mthd <- "z"
  }

  param <- cbind(param, se, tval, pval)
  dimnames(param) <- list(pnames, c(
    "Estimate", "Std. Error",
    "t value", "Pr(>|t|)"
  ))
  ans <- list(
    formula = formula(object$call$formula), residuals = r, sigma = sqrt(resvar),
    df = c(p, rdf), cov.unscaled = V / resvar, call = object$call,
    coefficients = param,
    r_squared = r_squared, lhs_start = object$lhs_start, lhs_end = object$lhs_end, class_lhs = class(object$lhs)
  )
  class(ans) <- "summary.midas_sp"
  ans
}

##' Semi-parametric MIDAS regression model deviance
##'
##' Returns the deviance of a fitted MIDAS regression object
##' @param object a \code{\link{midas_r}} object
##' @param ... currently nothing
##' @return The sum of squared residuals
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @rdname deviance.midas_sp
##' @method deviance midas_sp
##' @export
deviance.midas_sp <- function(object, ...) {
  sum(residuals(object)^2, na.rm = TRUE)
}

##' @export
##' @method print summary.midas_sp
print.summary.midas_sp <- function(x, digits = max(3, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"), ...) {
  cat(paste("\nSemi-parametric MIDAS regression model with \"", x$class_lhs[1],
    "\" data:\n",
    sep = ""
  ))
  cat(paste("Start = ", x$lhs_start,
    ", End = ", x$lhs_end,
    "\n",
    sep = ""
  ))
  cat("\n Formula:", deparse(x$formula), "\n")
  df <- x$df
  rdf <- df[2L]
  cat("\n Parameters:\n")
  printCoefmat(x$coefficients, digits = digits, signif.stars = signif.stars, ...)
  cat("\n Residual standard error:", format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
  #   cat(" Multiple R-squared: ", formatC(x$r_squared, digits = digits))
  #   cat(",\tAdjusted R-squared: ", formatC(x$adj_r_squared, digits = digits),"\n")
  invisible(x)
}


##' @include midas_r_methods.R
setMethod("extract", signature = className("midas_sp", "midasr"), definition = extract.midas_r)

##' @method update midas_sp
##' @importFrom stats getCall setNames
##' @export
update.midas_sp <- function(object, formula., ..., evaluate = TRUE) {
  if (is.null(call <- getCall(object))) {
    stop("need an object with call component")
  }
  extras <- match.call(expand.dots = FALSE)$...
  if (!missing(formula.)) {
    call$formula <- update(Formula(object), formula.)
  }

  if (length(extras)) {
    existing <- !is.na(match(names(extras), names(call)))
    for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
    if (any(!existing)) {
      call <- c(as.list(call), extras[!existing])
      call <- as.call(call)
    }
  }

  ## 1. If no start is supplied update the start from the call
  ## 2. If start is supplied intersect it with already fitted values.

  cf <- coef(object)
  ustart <- lapply(object$term_info, function(x) cf[x[["coef_index"]]])

  redo <- FALSE
  if (!("start" %in% names(extras))) {
    if (!("start" %in% names(call) && is.null(call$start))) {
      call$start <- ustart
      object$start_opt <- cf
    }
  } else {
    cstart <- eval(call$start, object$Zenv)
    ustart[names(cstart)] <- cstart
    call$start <- ustart
    object$start_opt <- unlist(ustart)
  }

  if (evaluate) {
    if (!missing(formula.) || "data" %in% names(extras) || "weight_gradients" %in% names(extras) || redo) {
      eval(call, parent.frame())
    } else {
      ## If we got here, we assume that we do not need to reevaluate terms.
      if (!is.null(extras$Ofunction)) {
        Ofunction <- eval(extras$Ofunction)
        extras$Ofunction <- NULL
      } else {
        Ofunction <- object$argmap_opt$Ofunction
      }
      dotargnm <- names(extras)
      if (length(dotargnm) > 0) {
        offending <- dotargnm[!dotargnm %in% names(formals(Ofunction))]
        if (length(offending) > 0) {
          stop(paste("The function ", Ofunction, " does not have the following arguments: ",
            paste(offending, collapse = ", "),
            sep = ""
          ))
        }
      }
      else {
        extras <- NULL
      }
      if (Ofunction != object$argmap_opt$Ofunction) {
        argmap <- c(list(Ofunction = Ofunction), extras)
      }
      else {
        argmap <- object$argmap_opt
        argmap$Ofunction <- NULL
        argnm <- union(names(argmap), names(extras))
        marg <- vector("list", length(argnm))
        names(marg) <- argnm
        marg[names(extras)] <- extras
        oldarg <- setdiff(names(argmap), names(extras))
        marg[oldarg] <- argmap[oldarg]
        argmap <- c(list(Ofunction = Ofunction), marg)
      }
      object$call <- call
      object$argmap_opt <- argmap
      midas_nlpr.fit(object)
    }
  }
  else {
    call
  }
}

##' Predicted values based on \code{midas_sp} object.
##'
##' \code{predict.midas_sp} produces predicted values, obtained by evaluating regression function in the frame \code{newdata}. This means that the appropriate model matrix is constructed using only the data in \code{newdata}. This makes this function not very convenient for forecasting purposes. If you want to supply the new data for forecasting horizon only use the function \link{forecast.midas_r}. Also this function produces only static predictions, if you want dynamic forecasts use the \link{forecast.midas_r}.
##'
##' @title Predict method for semi-parametric MIDAS regression fit
##' @param object \code{\link{midas_nlpr}} object
##' @param newdata a named list containing data for mixed frequencies. If omitted, the in-sample values are used.
##' @param na.action function determining what should be done with missing values in \code{newdata}. The most likely cause of missing values is the insufficient data for the lagged variables. The default is to omit such missing values.
##' @param ... additional arguments, not used
##' @return a vector of predicted values
##' @author Virmantas Kvedaras, Vaidotas Zemlys-Balevičius
##' @method predict midas_sp
##' @rdname predict.midas_sp
##' @export
predict.midas_sp <- function(object, newdata, na.action = na.omit, ...) {
  Zenv <- new.env(parent = parent.frame())

  if (missing(newdata)) {
    return(as.vector(fitted(object)))
  } else {
    ee <- data_to_env(newdata)
    ZZ <- object$Zenv
    parent.env(ee) <- ZZ
  }

  assign("ee", ee, Zenv)
  cll <- match.call()
  mf <- match.call(expand.dots = FALSE)
  m <- match(c("object", "newdata"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")

  formula <- Formula(formula(object))
  Terms <- delete.response(formula)
  mf[[2L]] <- Terms
  mf[[3L]] <- as.name("ee")
  mf[[4L]] <- na.action
  names(mf)[c(2, 3, 4)] <- c("formula", "data", "na.action")

  mf <- eval(mf, Zenv)
  mt <- attr(mf, "terms")



  Xeval <- model.matrix(formula, data = mf, rhs = 1)
  if (length(attr(formula, "rhs")) > 1) {
    Zeval <- model.matrix(formula, data = mf, rhs = 2)
    if (attr(object$terms2, "intercept") == 1) {
      Zeval <- Zeval[, -1, drop = FALSE]
    }
    if (attr(object$terms, "intercept") == 1) {
      Xeval <- Xeval[, -1, drop = FALSE]
    }
  }
  else {
    if (attr(object$terms2, "intercept") == 1) {
      Zeval <- Xeval[, -1, drop = FALSE]
    }
    Xeval <- NULL
  }

  as.vector(midas_sp_fit(object, Xeval, Zeval)$fitted.values)
}

##' Extracts various coefficients of MIDAS regression
##'
##' MIDAS regression has two sets of cofficients. The first set is the coefficients associated with the parameters
##' of weight functions associated with MIDAS regression terms. These are the coefficients of the NLS problem associated with MIDAS regression.
##' The second is the coefficients of the linear model, i.e  the values of weight
##' functions of terms, or so called MIDAS coefficients. By default the function returns the first set of the coefficients.
##'
##' @title Extract coefficients of MIDAS regression
##' @param object \code{midas_nlpr} object
##' @param type one of plain, midas, or nlpr. Returns appropriate coefficients.
##' @param term_names a character vector with term names. Default is \code{NULL}, which means that coefficients of all the terms are returned
##' @param ... not used currently
##' @return a vector with coefficients
##' @author Vaidotas Zemlys
##' @method coef midas_sp
##' @rdname coef.midas_sp
##' @export
coef.midas_sp <- function(object, type = c("plain", "midas", "bw"), term_names = NULL, ...) {
  type <- match.arg(type)
  if (is.null(term_names) & type == "midas") stop("Please provide a term name to get midas type coefficients")
  if (is.null(term_names)) {
    if (type == "bw") {
      object$coefficients[1:length(object$bws)]
    }
    else {
      return(object$coefficients)
    }
  } else {
    if (type == "bw") {
      warning("Bandwith setting is the same for all terms")
      coef(object, type = "bw", term_names = NULL)
    }
    if (length(setdiff(term_names, names(object$term_info))) > 0) stop("Some of the term names are not present in estimated MIDAS regression")
    if (type == "plain") {
      res <- lapply(object$term_info[term_names], function(x) {
        if (!is.null(x$coef_index)) object$coefficients[x$coef_index]
      })
    }
    if (type == "midas") {
      res <- lapply(object$term_info[term_names], function(x) {
        if (!is.null(x$coef_index)) x$weight(object$coefficients[x$coef_index])
      })
    }
    names(res) <- NULL
    if (length(res) == 1) {
      return(res[[1]])
    } else {
      return(unlist(res))
    }
  }
}

##' @include midas_r_methods.R
##' @include midas_nlpr_methods.R
##' @method plot_midas_coef midas_sp
##' @export
plot_midas_coef.midas_sp <- plot_midas_coef.midas_nlpr

##' Plot non-parametric part of the single index MIDAS regression
##'
##' Plot non-parametric part of the single index MIDAS regression
##' of unrestricted MIDAS regression
##' @param x \code{midas_r} object
##' @param term_name the term name for which the coefficients are plotted. Default is \code{NULL}, which selects the first MIDAS term
##' @param title the title string of the graph. The default is \code{NULL} for the default title.
##' @param compare the parameters for weight function to compare with the model, default is NULL
##' @param ... not used
##' @return a data frame with restricted MIDAS coefficients, unrestricted MIDAS coefficients and lower and upper confidence interval limits. The data
##' frame is returned invisibly.
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @importFrom graphics plot points
##' @importFrom numDeriv jacobian
##' @importFrom stats na.omit approx density
##' @export
plot_sp <- function(x, term_name, title = NULL, compare = NULL, ...) {
  if (length(x$bws) > 1) stop("Plotting is currently supported for univariate non-parametric functions only")
  ti <- x$term_info[[term_name]]

  if (ti$term_type != "Z") stop("The term name supplied does not correspond to non-parametric term")

  gfun <- function(p) {
    if (is.null(x$model_matrix_map$X)) {
      xi <- 0
    } else {
      xi <- x$model[, x$model_matrix_map$X] %*% x$coef_X(p)
    }
    Z <- x$model[, x$model_matrix_map$Z, drop = FALSE]
    if (is.null(ti$coef_index)) {
      zi <- Z[, ti$midas_coef_index, drop = FALSE]
    } else {
      zi <- Z[, ti$midas_coef_index] %*% ti$weight(p[ti$coef_index])
    }
    if (ncol(zi) == 1) zi <- as.numeric(zi)
    u <- x$model[, 1] - xi
    list(xi = xi, zi = zi, u = u, g = cv_np(u, zi, p[1:length(x$bws)], x$degree))
  }


  gg <- gfun(coef(x))

  kappa <- 0.5 / pi # for standard Gaussian kernel
  dns <- density(gg$zi, na.rm = T) # choose a better bandwith here !!!
  f_z <- approx(dns$x, dns$y, xout = gg$zi)$y
  se_np <- sqrt(kappa * deviance(x) / f_z / (exp(coef(x)[1]) * x$nobs))


  ozi <- order(gg$zi)

  se_npo <- se_np[ozi]
  xio <- gg$zi[ozi]
  trmo <- gg$g[ozi]
  pd <- data.frame(xi = xio, term = trmo, compare = NA, lower = trmo - 1.96 * se_npo, upper = trmo + 1.96 * se_npo)

  if (!is.null(compare)) {
    pd$compare <- compare(pd$xi)
  }

  ylim <- range(na.omit(c(pd$term, pd$compare, pd$lower, pd$upper)))

  plot(pd$xi, pd$term, col = "blue", ylab = "Estimated non-parametric function", xlab = "MIDAS aggregate", ylim = ylim, type = "l")

  points(pd$xi, pd$compare, type = "l", col = "black")
  points(pd$xi, pd$lower, type = "l", col = "grey", lty = 2)
  points(pd$xi, pd$upper, type = "l", col = "grey", lty = 2)

  if (is.null(title)) {
    title(main = paste0("Non-parametric estimate for term ", term_name, " with weight: ", ti$weight_name))
  } else {
    title(main = title)
  }

  invisible(pd)
}