R/imidasreg.R

In mpiktas/midasr: Mixed Data Sampling Regression

##' Restricted MIDAS regression with I(1) regressors
##'
##' Estimate restricted MIDAS regression using non-linear least squares, when the regressor is I(1)
##'
##' @param formula formula for restricted MIDAS regression. Formula must include \code{\link{fmls}} function
##' @param data a named list containing data with mixed frequencies
##' @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 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 ... 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{z_env}{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}
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @rdname imidas_r
##' @seealso midas_r.midas_r
##' @examples
##' 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)
##' }
##'
##' theta0 <- theta.h0(c(-0.1,10,-10,-10),4*12)
##'
##' xx <- ts(cumsum(rnorm(600*12)), frequency = 12)
##'
##' ##Simulate the response variable
##' y <- midas_sim(500, xx, theta0)
##'
##' x <- window(xx, start=start(y))
##'
##' imr <- imidas_r(y~fmls(x,4*12-1,12,theta.h0)-1,start=list(x=c(-0.1,10,-10,-10)))
##'
##' @details Given MIDAS regression:
##'
##' \deqn{y_t=\sum_{j=0}^k\sum_{i=0}^{m-1}\theta_{jm+i} x_{(t-j)m-i}+\mathbf{z_t}\beta+u_t}
##'
##' estimate the parameters of the restriction
##'
##' \deqn{\theta_h=g(h,\lambda),}
##' where \eqn{h=0,...,(k+1)m}, together with coefficients \eqn{\beta} corresponding to additional
##' low frequency regressors.
##'
##' It is assumed that \eqn{x} is a I(1) process, hence the special transformation is made.
##' After the transformation \link{midas_r} is used for estimation.
##'
##' MIDAS regression involves times series with different frequencies.
##'
##' The restriction function must return the restricted coefficients of
##' the MIDAS regression.
##'
##' @export
imidas_r <- function(formula, data, start, Ofunction = "optim", weight_gradients = NULL, ...) {
  z_env <- new.env(parent = environment(formula))

  mt <- terms(formula(formula), specials = "fmls")


  if (missing(data)) {
    nms <- all.vars(mt)
    insample <- lapply(nms, function(nm) eval(as.name(nm), z_env))
    names(insample) <- nms
    insample <- insample[!sapply(insample, is.function)]
  }
  else {
    insample <- data_to_list(data)
  }

  vl <- as.list(attr(mt, "variables"))
  vl <- vl[-1]

  pl <- attr(mt, "specials")$fmls
  if (length(pl) > 1) stop("Only one high frequency term is supported currently")
  fr <- vl[[pl]]

  ## Do an ugly hack
  xname <- as.character(fr[[2]])
  insample <- c(list(insample[[xname]]), insample)
  names(insample)[1] <- ".Level"

  assign(".IMIDASdata", insample, envir = z_env)

  wf <- fr[-4:-5]
  wf[[1]] <- fr[[5]]
  for (j in 3:length(wf)) {
    wf[[j]] <- eval(wf[[j]], z_env)
  }
  wf[[3]] <- wf[[3]] + 1
  pp <- function(p, d) {
    wf[[2]] <- p
    r <- eval(wf, z_env)
    cumsum(r)[1:d]
  }
  mfg <- wf
  mfg[[1]] <- as.name(paste(as.character(wf[[1]]), "gradient", sep = "."))
  pp.gradient <- function(p, d) {
    mfg[[2]] <- p
    r <- eval(mfg, z_env)
    apply(r, 2, cumsum)[1:d, ]
  }

  formula <- expandfmls(formula(formula), "pp", z_env, 0)
  cl <- match.call(expand.dots = TRUE)
  cl <- cl[names(cl) != "model"]

  assign("pp", pp, z_env)
  assign("pp.gradient", pp.gradient, z_env)
  environment(formula) <- z_env
  cl[[2]] <- formula
  cl[[1]] <- as.name("midas_r")
  cl$data <- as.name(".IMIDASdata")
  res <- eval(cl, z_env)
  class(res) <- c(class(res), "imidas_r")
  return(res)
}

## Function for expanding the formula in I(1) case
expandfmls <- function(expr, wfun, z_env, diff = 0, truncate = FALSE) {
  if (length(expr) == 3) {
    expr[[2]] <- expandfmls(expr[[2]], wfun, z_env, diff, truncate)
    expr[[3]] <- expandfmls(expr[[3]], wfun, z_env, diff, truncate)
  }
  if (length(expr) == 5) {
    if (expr[[1]] == as.name("fmls")) {
      rr <- modifyfmls(expr, wfun, z_env, diff)
      if (truncate) {
        return(rr[[2]])
      }
      else {
        return(rr)
      }
    }
    else {
      return(expr)
    }
  }
  return(expr)
}

## Substitute fmls with dmls
modifyfmls <- function(expr, wfun, z_env, diff) {
  res <- expression(a + b)[[1]]
  t2 <- expression(mls(x, a, b))[[1]]
  nol <- eval(expr[[3]], z_env)
  expr[[3]] <- nol + diff
  m <- eval(expr[[4]], z_env)
  expr[[1]] <- as.name("dmls")
  t2[[2]] <- as.name(".Level")
  t2[[3]] <- nol + 1 + diff
  t2[[4]] <- m
  expr[[5]] <- as.name(wfun)
  res[[2]] <- expr
  res[[3]] <- t2
  res
}