R/hpj.R

In jumps: Hodrick-Prescott Filter with Jumps

# This script contains the main function to be run by the user.
# It is a wrapper for the various Hodrick-Prescott filters with jumps functions.
# Moreover, it contains the methods for the hpj object generated by the hpj() function.


#' Hodrick-Prescott filter with jumps
#' 
#' @description
#' This function is a wrapper for the various Hodrick-Prescott filters with jumps functions.
#' The end user should use this: it is simpler and more flexible than the other functions.
#' 
#' @param y either a numeric vector or a time series object containing the time series to filter;
#' @param maxsum maximum sum of additional level standard deviations,
#' if \code{NULL} the value is selected automatically on a grid using the specified information criterion (ic);
#' @param lambda smoothing constant of the HP filter,
#' if \code{NULL} the value is estimated by maximum likelihood;
#' @param xreg matrix of regressors;
#' @param ic string with information criterion for the automatic choice of \code{maxsum}:
#' the default is "bic" (simulations show this is the best choice), but also "hq", "aic" and "aicc"
#' are available.
#' @param scl scaling factor for the time series (default is 1000): the time series
#' is rescaled as (y-min(y))/(max(y)-min(y))*scl. This is done since the default
#' starting values for the optimization seem to work well in this scale);
#' If `scl` is set equal to the string `"original"` the time series is not rescaled.
#' 
#' @returns S3 object of class hpj with the following slots:
#' \itemize{
#'  \item y: the input time series;
#'  \item maxsum: the maximum sum of additional standard deviations;
#'  \item lambda: the smoothing constant of the HP filter;
#'  \item pars: vector of estimated parameters (sigma_slope, sigma_noise, gamma);
#'  \item hpj: the time series of the HP filter with jumps;
#'  \item hpj_std: the time series of the HP filter with jumps standard deviations;
#'  \item std_devs: vector of additional standard deviations of the level disturbance;
#'  \item breaks: vector of indices of the breaks;
#'  \item xreg: matrix of regressors;
#'  \item df: model's degrees of freedom;
#'  \item loglik: value of the log-likelihood at maximum;
#'  \item ic: vector of information criteria (aic, aicc, bic, hq);
#'  \item opt: the output of the optimization function (nloptr);
#'  \item call: the call to the function.
#' }
#' @examples 
#' set.seed(202311)
#' n <- 100
#' mu <- 100*cos(3*pi/n*(1:n)) - ((1:n) > 50)*n - c(rep(0, 50), 1:50)*10
#' y <- mu + rnorm(n, sd = 20)
#' hp <- hpj(y, 50)
#' plot(hp)
#' @export
hpj <- function(y, maxsum = NULL, lambda = NULL, xreg = NULL,
                ic = c("bic", "hq", "aic", "aicc"), scl = 1000) {
  if (scl == "original") scl <- max(y) - min(y)
  y_num <- as.numeric(y)
  y_min <- min(y_num)
  y_max <- max(y_num)
  y_rng <- y_max - y_min
  y_scl <- (y_num - y_min) / y_rng * scl # rescaled y
  if (!is.null(xreg)) {
    stop("xreg is not supported yet")
  }
  if (!is.null(lambda)) {  # lambda is specified
    if (is.null(maxsum)) { # no maxsum is provided
      out <- auto_hpfj_fix(y = y_scl, lambda = lambda, ic = ic)
    } else { # maxsum is provided
      out <- hpfj_fix(y = y_scl, lambda = lambda, maxsum = maxsum)
    }
  } else {
    if (is.null(maxsum)) { # no maxsum is provided
      out <- auto_hpfj(y = y_scl, ic = ic)
    } else { # maxsum is provided
      out <- hpfj(y = y_scl, maxsum = maxsum)
    }
  }
  # antitransform
  smo <- y
  smo[] <- out$smoothed_level/scl*y_rng + y_min
  smo_std <- y
  smo_std[] <- sqrt(out$var_smoothed_level)/scl*y_rng
  sigmas <- y
  sigmas[] <- out$sigmas/scl*y_rng
  sigmas[1] <- 0
  pars <- out$pars
  pars[1:2] <- pars[1:2]/scl*y_rng
  # names(pars) <- c("sigma_slope", "sigma_noise", "gamma")
  if (is.null(maxsum)) maxsum = out$maxsum/scl*y_rng
  if (is.null(lambda)) lambda = unname(pars[2]/pars[1])^2
  structure(
    list(
      y = y,
      maxsum = maxsum,
      lambda = lambda,
      pars = pars,
      hpj = smo,
      hpj_std = smo_std,
      breaks = which(out$sigmas > scl/100000),
      sigmas = sigmas,
      xreg = xreg,
      df = out$df,
      loglik = out$loglik,
      ic = out$ic,
      opt = out$opt,
      call = deparse(sys.call())
    ),
    class = "hpj"
  )
}


#' Plot method for the class hpj
#' 
#' @param x an object of class hpj;
#' @param prob coverage probability for the confidence interval of the filter;
#' @param show_breaks logical, if \code{TRUE} vertical lines are drawn at the jumps;
#' @param main title of the plot;
#' @param use_ggplot logical, if \code{TRUE} the plot is done with ggplot2;
#' @param ... additional arguments passed to the plot function (if no ggplot2 is used).
#' @returns If \code{use_ggplot} is \code{TRUE} a ggplot object is returned, otherwise
#' a base plot is shown and nothing is returned.
#' @rdname plot
#' @exportS3Method base::plot
#' @export
plot.hpj <- function(x, prob = NULL, show_breaks = TRUE, main = "original + filter",
                     use_ggplot = TRUE, ...) {
  if (use_ggplot & ("timeSeries" %in% class(x$y))) {
    warning("timeSeries objects are not supported for ggplot2")
    use_ggplot <- FALSE
  }
  value <- x$y
  filtered <- x$hpj
  if (use_ggplot) { # if use_ggplot is TRUE
    dt <- data.frame(time = time(value),
                     value = as.numeric(value),
                     filtered = as.numeric(filtered))
    gg <- ggplot2::ggplot(dt, ggplot2::aes(x = time, y = value)) +
      ggplot2::geom_line() +
      ggplot2::geom_line(ggplot2::aes(y = filtered), color = "red") +
      ggplot2::ggtitle(main)
    if (!is.null(prob)) {
      gg <- gg +
        ggplot2::geom_ribbon(ggplot2::aes(ymax = filtered - qnorm((1-prob)/2)*x$hpj_std,
                        ymin = filtered + qnorm((1-prob)/2)*x$hpj_std),
                        fill = "red", alpha = 0.1)
    }
    if (show_breaks & length(x$breaks)) {
      gg <- gg + ggplot2::geom_vline(xintercept = time(value)[x$breaks], linetype = "dashed")
    }
    return(gg)
  }
  # if use_ggplot is FALSE
  if (is.null(prob)) {
    plot(value, type = "l", main = main, ...)
    lines(filtered, col = "red", ...)
  } else {
    z <- abs(qnorm((1-prob)/2))
    up <- filtered + z*x$hpj_std
    lo <- filtered - z*x$hpj_std
    rng <- range(c(as.numeric(value), as.numeric(up), as.numeric(lo)))
    plot(value, type = "l", ylim = rng, main = main, ...)
    lines(filtered, col = "red", ...)
    lines(up, col = "red", lty = 2, ...)
    lines(lo, col = "red", lty = 2, ...)
  }
  if (show_breaks & length(x$breaks)) {
    if ("xts" %in% class(value)) {
      evt <- xts::xts(rep("", length(x$breaks)), time(value)[x$breaks])
      xts::addEventLines(evt, lty = 2)
    } else {
      abline(v = time(value)[x$breaks], lty = 2)
    }
  }
}


#' Print method for the class hpj
#' 
#' @param x an object of class hpj;
#' @param ... not used: for consistency with generic function.
#' @returns No return value, called for side effects
#' @rdname print
#' @exportS3Method base::print
#' @export
print.hpj <- function(x, ...) {
  cat("Hodrick-Prescott filter with jumps\n")
  cat("Call:\n")
  cat(" ", x$call, "\n")
  cat("Parameters:\n")
  cat("  sd(slope) = ", x$pars[1], "\n", sep = "")
  cat("  sd(noise) = ", x$pars[2], "\n", sep = "")
  cat("  gamma     = ", x$pars[3], "\n", sep = "")
  cat("  lambda    = ", x$lambda, "\n", sep = "")
  cat("  maxsum    = ", x$maxsum, "\n", sep = "")
  cat("Model's degrees of freedom: ", x$df, "\n")
  cat("Log-likelihood: ", x$loglik, "\n")
  cat("Information criteria:\n")
  cat("  AIC  = ", x$ic[1], "\n", sep = "")
  cat("  AICc = ", x$ic[2], "\n", sep = "")
  cat("  BIC  = ", x$ic[3], "\n", sep = "")
  cat("  HQ   = ", x$ic[4], "\n", sep = "")
  cat("Break dates:")
  if (length(x$breaks) == 0) {
    cat("  no breaks were detected\n")
  } else {
    M <- as.matrix(as.character(time(x$y)[x$breaks]))
    rownames(M) <- 1:length(x$breaks)
    colnames(M) <- ""
    print(M, quote = FALSE)
  }
}

#' logLik method for the class hpj
#' 
#' @param object an object of class hpj;
#' @param ... not used: for consistency with generic function.
#' @returns An object of logLik class with the log-likelihood value and two
#' attributes: df, the number of degrees of freedom, and nobs, the number of
#' observations.
#' @rdname logLik
#' @exportS3Method base::logLik
#' @export
logLik.hpj <- function(object, ...) {
  structure(
    object$loglik,
    class = "logLik",
    df = object$df,
    nobs = sum(!is.na(object$y))
  )
}

#' nobs method for the class hpj
#' 
#' @param object an object of class hpj;
#' @param ... not used: for consistency with generic function.
#' @returns The number of (non missing) observations of the time series on
#' which the HP filter with jumps has been applied.
#' @rdname nobs
#' @exportS3Method base::nobs
#' @export
nobs.hpj <- function(object, ...) {
  sum(!is.na(object$y))
}

#' BIC method for the class hpj
#' 
#' @param object an object of class hpj;
#' @param ... additional objects of class hpj.
#' @returns If just one object is provided,
#' a numeric value with the corresponding BIC.
#' If multiple objects are provided, a data.frame with rows corresponding to
#' the objects and columns representing the number of parameters in the
#' model (df) and the BIC.
#' @rdname BIC
#' @exportS3Method base::BIC
#' @export
BIC.hpj <- function(object, ...) {
  dots <- list(...)
  if (length(dots) == 0) return(unname(object$ic["bic"]))
  ndx <- which(sapply(dots, function(x) class(x) == "hpj"))
  if (length(ndx) == 0) return(unname(object$ic["bic"]))
  df <- data.frame(df = c(object$df, sapply(dots[ndx], function(x) x$df)),
                   BIC = c(object$ic["bic"], sapply(dots[ndx], function(x) x$ic["bic"])))
  arg_names <- as.character(substitute(list(object, ...)))[-1]
  rownames(df) <- c(arg_names[1], arg_names[-1][ndx])
  df
}