R/midasqr.R

In mpiktas/midasr: Mixed Data Sampling Regression

##' Restricted MIDAS quantile regression
##'
##' Estimate restricted MIDAS quantile regression using nonlinear quantile regression
##'
##' @param formula formula for restricted MIDAS regression or \code{midas_qr} object. Formula must include \code{\link{mls}} function
##' @param data a named list containing data with mixed frequencies
##' @param tau quantile
##' @param start the starting values for optimisation. Must be a list with named elements.
##' @param Ofunction the list with information which R function to use for optimisation. The list must have element named \code{Ofunction} which contains character string of chosen R function. Other elements of the list are the arguments passed to this function.  The default optimisation function is \code{\link{optim}} with argument \code{method="BFGS"}. Other supported functions are \code{\link{nls}}
##' @param weight_gradients a named list containing gradient functions of weights. The weight gradient function must return the matrix with dimensions
##' \eqn{d_k \times q}, where \eqn{d_k} and \eqn{q} are the number of coefficients in unrestricted and restricted regressions correspondingly.
##' The names of the list should coincide with the names of weights used in formula.
##' The default value is NULL, which means that the numeric approximation of weight function gradient is calculated. If the argument is not NULL, but the
##' name of the weight used in formula is not present, it is assumed that there exists an R function which has
##' the name of the weight function appended with \code{_gradient}.
##' @param guess_start, logical, if TRUE tries certain strategy to improve starting values
##' @param ... additional arguments supplied to optimisation function
##' @return a \code{midas_r} object which is the list with the following elements:
##'
##' \item{coefficients}{the estimates of parameters of restrictions}
##' \item{midas_coefficients}{the estimates of MIDAS coefficients of MIDAS regression}
##' \item{model}{model data}
##' \item{unrestricted}{unrestricted regression estimated using \code{\link{midas_u}}}
##' \item{term_info}{the named list. Each element is a list with the information about the term, such as its frequency, function for weights, gradient function of weights, etc.}
##' \item{fn0}{optimisation function for non-linear least squares problem solved in restricted MIDAS regression}
##' \item{rhs}{the function which evaluates the right-hand side of the MIDAS regression}
##' \item{gen_midas_coef}{the function which generates the MIDAS coefficients of MIDAS regression}
##' \item{opt}{the output of optimisation procedure}
##' \item{argmap_opt}{the list containing the name of optimisation function together with arguments for optimisation function}
##' \item{start_opt}{the starting values used in optimisation}
##' \item{start_list}{the starting values as a list}
##' \item{call}{the call to the function}
##' \item{terms}{terms object}
##' \item{gradient}{gradient of NLS objective function}
##' \item{hessian}{hessian of NLS objective function}
##' \item{gradD}{gradient function of MIDAS weight functions}
##' \item{Zenv}{the environment in which data is placed}
##' \item{use_gradient}{TRUE if user supplied gradient is used, FALSE otherwise}
##' \item{nobs}{the number of effective observations}
##' \item{convergence}{the convergence message}
##' \item{fitted.values}{the fitted values of MIDAS regression}
##' \item{residuals}{the residuals of MIDAS regression}
##'
##' @examples
##' ##Take the same example as in midas_r
##'
##' theta_h0 <- function(p, dk, ...) {
##'    i <- (1:dk-1)/100
##'    pol <- p[3]*i + p[4]*i^2
##'    (p[1] + p[2]*i)*exp(pol)
##' }
##'
##' ##Generate coefficients
##' theta0 <- theta_h0(c(-0.1,10,-10,-10),4*12)
##'
##' ##Plot the coefficients
##' plot(theta0)
##'
##' ##Generate the predictor variable
##' xx <- ts(arima.sim(model = list(ar = 0.6), 600 * 12), frequency = 12)
##'
##' ##Simulate the response variable
##' y <- midas_sim(500, xx, theta0)
##'
##' x <- window(xx, start=start(y))
##'
##' ##Fit quantile regression. All the coefficients except intercept should be constant.
##' ##Intercept coefficient should correspond to quantile function of regression errors.
##' mr <- midas_qr(y~fmls(x,4*12-1,12,theta_h0), tau = c(0.1, 0.5, 0.9),
##'               list(y=y,x=x),
##'               start=list(x=c(-0.1,10,-10,-10)))
##'
##' mr
##' @author Vaidotas Zemlys-Balevicius
##' @rdname midas_qr
##' @import quantreg
##' @export
midas_qr <- function(formula, data, tau = 0.5, start, Ofunction = "nlrq", weight_gradients = NULL, guess_start = TRUE, ...) {
  Zenv <- new.env(parent = environment(formula))

  if (missing(data)) {
    ee <- NULL
  }
  else {
    ee <- data_to_env(data)
    parent.env(ee) <- parent.frame()
  }

  if (missing(start)) {
    stop("Please supply starting values.")
  }

  assign("ee", ee, Zenv)
  formula <- as.formula(formula)
  cl <- match.call()
  mf <- match.call(expand.dots = FALSE)
  mf$formula <- formula
  m <- match(c("formula", "data"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf[[1L]] <- as.name("model.frame")
  mf[[3L]] <- as.name("ee")
  mf[[4L]] <- as.name("na.omit")
  names(mf)[c(2, 3, 4)] <- c("formula", "data", "na.action")

  ## Add check for dynamic terms.
  ## They are not supported as they do not make sense in quantile regression
  itr <- checkARstar(terms(eval(mf[[2]], Zenv)))
  if (!is.null(itr$lagsTable)) mf[[2]] <- itr$x

  mf <- eval(mf, Zenv)
  mt <- attr(mf, "terms")
  args <- list(...)
  y <- model.response(mf, "numeric")
  X <- model.matrix(mt, mf)

  if (is.null(ee)) {
    yy <- eval(formula[[2]], Zenv)
  } else {
    yy <- eval(formula[[2]], ee)
  }

  y_index <- 1:length(yy)
  if (!is.null(attr(mf, "na.action"))) {
    y_index <- y_index[-attr(mf, "na.action")]
  }
  if (length(y_index) > 1) {
    if (sum(abs(diff(y_index) - 1)) > 0) warning("There are NAs in the middle of the time series")
  }

  ysave <- yy[y_index]

  if (inherits(yy, "ts")) {
    class(ysave) <- class(yy)
    attr(ysave, "tsp") <- c(time(yy)[range(y_index)], frequency(yy))
  }

  if (inherits(yy, c("zoo", "ts"))) {
    y_start <- index2char(index(ysave)[1], frequency(ysave))
    y_end <- index2char(index(ysave)[length(ysave)], frequency(ysave))
  } else {
    y_start <- y_index[1]
    y_end <- y_index[length(y_index)]
  }

  eps <- .Machine$double.eps^(2 / 3)
  if (any(tau < 0) || any(tau > 1)) {
    stop("invalid tau:  taus should be >= 0 and <= 1")
  }

  if (any(tau == 0)) {
    tau[tau == 0] <- eps
  }
  if (any(tau == 1)) {
    tau[tau == 1] <- 1 - eps
  }

  prepmd <- prepmidas_r(y, X, mt, Zenv, cl, args, start, Ofunction, weight_gradients, itr$lagsTable, guess_start = guess_start, tau = tau)

  prepmd <- c(prepmd, list(lhs = ysave, lhs_start = y_start, lhs_end = y_end))

  class(prepmd) <- c("midas_qr")

  midas_qr.fit(prepmd)
}

midas_qr.fit <- function(x) {
  args <- x$argmap_opt
  function.opt <- args$Ofunction
  args$Ofunction <- NULL
  if (function.opt == "nlrq") {
    rhs <- x$rhs
    if (x$use_gradient) {
      orhs <- rhs
      rhs <- function(p) {
        res <- orhs(p)
        attr(res, "gradient") <- x$model[, -1] %*% x$gradD(p)
        res
      }
    }
    z <- x$model[, 1]
    args$formula <- formula(z ~ rhs(p))
    args$start <- list(p = x$start_opt)
    res <- vector(mode = "list", length = length(x$tau))
    for (i in 1:length(x$tau)) {
      args$tau <- x$tau[i]
      opt <- try(do.call("nlrq", args), silent = TRUE)
      if (inherits(opt, "try-error")) {
        stop("The optimisation algorithm of MIDAS regression failed with the following message:\n", opt, "\nPlease try other starting values or a different optimisation function")
      }
      par <- coef(opt)
      names(par) <- names(coef(x))
      res[[i]] <- list(opt = opt, par = par, convergence = opt$convInfo$stopCode)
    }
  }
  if (function.opt == "dry_run") {
    opt <- NULL
    res_template <- list(opt = opt, par = x$start_opt, convergence = "Dry run, no optimisation done")
    res <- vector(mode = "list", length = length(x$tau))
    for (i in 1:length(x$tau)) {
      res[[i]] <- res_template
    }
  }
  if (length(res) == 1) {
    x$opt <- res[[1]]$opt
    x$coefficients <- res[[1]]$par
    names(res[[1]]$par) <- NULL
    x$midas_coefficients <- x$gen_midas_coef(res[[1]]$par)
    x$fitted.values <- as.vector(x$model[, -1] %*% x$midas_coefficients)
    x$residuals <- as.vector(x$model[, 1] - x$fitted.values)
  } else {
    x$opt <- lapply(res, "[[", "opt")
    x$convergence <- sapply(res, "[[", "convergence")
    x$coefficients <- sapply(res, "[[", "par")
    colnames(x$coefficients) <- x$tau
    par <- sapply(res, "[[", "par")
    rownames(par) <- NULL
    x$midas_coefficients <- apply(par, 2, x$gen_midas_coef)
    colnames(x$midas_coefficients) <- x$tau
    x$fitted.values <- x$model[, -1] %*% x$midas_coefficients
    x$residuals <- x$model[, 1] - x$fitted.values
  }
  x
}