R/dynadj.R

In nelson16silva/wavcoreinf: Estimation and Evaluation of Wavelet Core Inflation Measures

#' Estimate OLS from Data Frame
#'
#' This function is just a convenient way to estimate an OLS model from a data
#' frame where the first variable is the dependent one of the model and the
#' remaining are the regressors. This function is used internally
#' by other functions.
#'
#' @param df A data frame or a tibble. First column of the data frame must be the \eqn{y}
#' variable in a regression of the type \eqn{y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + ... + \beta_nx_{in} + \epsilon_i}, for example.
#'
#' @return An object of class \code{lm}.
#' @export
#'
#' @examples
#' coreinf_br %>%
#' dplyr::select(ipca, ipcams) %>%
#' comb_regr()
#'
#' coreinf_br %>%
#'  dplyr::select(ipca, ipcadp) %>%
#'  dplyr::mutate(`ipca(-1)` = dplyr::lag(ipca)) %>%
#'  comb_regr()
comb_regr <- function(df) {
  yx_names <- rlang::syms(names(df))
  exogenous <- purrr::reduce(yx_names[-1], ~ rlang::expr(!!.x + !!.y))
  equation <- rlang::expr(!!yx_names[[1]] ~ !!exogenous)
  lm_call <- (rlang::quo(stats::lm(!!equation)))
  rlang::eval_tidy(lm_call, df)
}

#' Estimate Models for Dynamic Adjustment Analysis
#'
#' Use \code{dyn_adj_est} to estimate several models for dynamic
#' adjustment analysis.
#'
#' @details For \eqn{\pi} being the healine inflation and \eqn{\pi^*}
#' a core inflation measure, two specifications of a regression model can be estimated through this
#' function:
#' \enumerate{
#' \item \eqn{\pi_{t + h} - \pi_t =
#'  a_0 + \lambda_h(\pi_t - \pi^*_t) + \sum_{i = 1}^pa_i\pi_{t - i} +
#'  e_{t + h}},
#'  \item  \eqn{\pi^*_{t + h} - \pi^*_t =
#'  a^*_0 + \lambda^*_h(\pi_t - \pi^*_t) + \sum_{i = 1}^pa^*_i\pi^*_{t - i} +
#'  e^*_{t + h}}.
#'}
#'
#' @param y A vector. If the user estimates the first regression described
#' below, \code{y} represents headline inflation or
#' core inflation measure if the regression is the second. The choice depends if interest is in
#' \eqn{\lambda_h} or in \eqn{\lambda^*_h}.
#' @param xreg A vector. Infation or core inflation depending of the \code{y} choice.
#' @param h An integer. The horizon of prediction.
#' @param p An integer to specify the lag p.
#' @param ... Additional parameter to pass to the function
#' \code{\link{lags}} for specifying if it is recursive or not.
#'
#' @return A list where each element of it contains the following:
#' \item{data}{A tibble with the data used in fitting the model.}
#' \item{model}{A \code{lm} object.}
#' \item{n_obs}{Number of observations of \code{y}.}
#' \item{h}{Horizon used in direct estimation.}
#'
#' @seealso \code{\link{lags}, \link{comb_regr}}
#'
#' @export
#'
#' @examples
#' inf_head <- coreinf_br[["ipca"]]
#' inf_corems <- coreinf_br[["ipcams"]]
#' dyn_adj_est(inf_head, inf_corems, 2, 2)
#' dyn_adj_est(inf_corems, inf_head, 2, 2)
dyn_adj_est <- function(y, xreg, h, p, ...) {
  yname <- rlang::enexpr(y)
  if(!is.name(yname)){
    yname <- "y"
  }
  xregname <- rlang::enexpr(xreg)
  if(!is.name(xregname)){
    xregname <- "x"
  }
  lags <- lags(p, q = 1, ...)[, 1:2]
  ydiffname <- rlang::expr(!!paste0(rlang::expr(!!yname), "_diff"))
  yxdiffname <- rlang::expr(!!paste0(rlang::expr(!!yname), "_minus_", rlang::expr(!!xregname)))
  diffy <- y - dplyr::lag(as.vector(y), h)
  diffyx <- dplyr::lag(as.vector(y), h) - dplyr::lag(as.vector(xreg), h)
  models <- purrr::pmap(lags, function(...) {
    lagp <- (..1:..2) + (h - 0)
    names(lagp) <- paste0(rlang::expr(!!yname), "_lag", lagp)
    ylag <- purrr::map_df(lagp, ~dplyr::lag(as.vector(y), .))
    ylag <- tibble::add_column(ylag, !!yxdiffname := diffyx, .before =1 )
    ylag <- tibble::add_column(ylag, !!ydiffname := diffy, .before = 1)
    data <- ylag
    model <- comb_regr(data[(max(lags) + h + 1):(length(y) - (0)), ])
    return(list(data = data, model = model, n_obs = length(y), h = h))
  })
  models
}


#' Select Best Dynamic Adjust Model
#'
#' Given models estimated by the function \code{\link{dyn_adj_est}}, \code{dyn_adj_best} selects
#' the best one.
#'
#' @param models An object generated by the funcion \code{\link{dyn_adj_est}}.
#' @param ic Information criterion BIC or AIC.
#'
#' @return A list containing the following elements:
#' \item{data}{A tibble with the data used in fitting the model.}
#' \item{model}{An \code{lm} object.}
#' \item{n_obs}{Number of observations.}
#' \item{h}{Horizon used in direct estimation. See \code{\link{dyn_adj_est}}.}
#'
#' @seealso \code{\link{dyn_adj_est}}
#'
#' @export
#'
#' @examples
#' inf_head <- coreinf_br[["ipca"]]
#' inf_corems <- coreinf_br[["ipcams"]]
#' models <- dyn_adj_est(inf_head, inf_corems, 2, 2)
#' dyn_adj_best(models, ic = AIC)
dyn_adj_best <- function(models, ic = BIC) {
  inf_crit <- purrr::map(purrr::map(models, ~`[[`(.x, "model")), ~ic(.x))
  best_model <- models[[which.min(inf_crit)]]
  rm(models)
  best_model
}


#' Extract the Dynamic Adjustment Parameter
#'
#' This function is a convenient way of extracting both the dynamic adjustment parameter
#' and the associated p-value from the best dynamic model estimated by the
#' function \code{\link{dyn_adj_best}}.
#'
#' @param best_model An object from function \code{\link{dyn_adj_best}}.
#'
#' @return A list of two elements: dynamic adjustment parameter and
#' p-value.
#'
#' @seealso \code{\link{dyn_adj_est}, \link{dyn_adj_best}}
#'
#' @export
#'
#' @examples
#' inf_head <- coreinf_br[["ipca"]]
#' inf_corems <- coreinf_br[["ipcams"]]
#' models <- dyn_adj_est(inf_head, inf_corems, 2, 2)
#' best_model <- dyn_adj_best(models, ic = AIC)
#' dyn_adj_pred(best_model)
dyn_adj_pred <- function(best_model) {
  lambda_adj <-  eval(best_model$model[["call"]], best_model[["data"]])
  lambda <- lambda_adj$coefficients[2]
  p_value <- stats::coef(summary(lambda_adj))[2, 4]
  return(list(alpha = lambda, val_p = p_value))
}


#' Map over Several Horizons to Extract Dynamic Adjustment Parameter
#'
#' \code{dyn_adj} is a fast way to extract the parameter
#' of the dynamic adjustment for a range of horizons.
#'
#' @details For a model \eqn{\pi_{t + h} - \pi_t =
#' a_0 + \lambda_h(\pi_t - \pi^*_t) + \sum_{i = 1}^pa_i\pi_{t - i} +
#' e_{t + h}}, where \eqn{\pi} is the headline inflation and
#' \eqn{\pi^*} a measure of core inflation,
#' this function gives \eqn{\lambda_h} for \eqn{h \in 1, 2 ..., H}.
#'
#' In core inflation's literature, parameters \eqn{\lambda_h}
#' and \eqn{\lambda^*_h} of the following regressions are important:
#' \itemize{
#' \item \eqn{\pi_{t + h} - \pi_t = a_0 + \lambda_h(\pi_t - \pi^*_t) +
#' \sum_{i = 1}^pa_i\pi_{t - i} + e_{t + h}},
#' \item \eqn{\pi^*_{t + h} - \pi^*_t = a^*_0 + \lambda^*_h(\pi_t - \pi^*_t) +
#' \sum_{i = 1}^pa^*_i\pi^*_{t - i} + e^*_{t + h}}.
#' }
#' A good core inflation measure should imply \eqn{\lambda_h < 0} and
#' \eqn{\lambda^*_h = 0}.
#'
#' @param inf A vector containing headline inflation.
#' @param core A vector containing a measure of core inflation.
#' @param H An integer that gives the maximum horizon.
#' @param p An integer, lag in the regression (see details).
#' @param ... Additional parameters.
#'
#' @return A tibble with the dynamic adjustment parameter and p-value for the t test
#' of \eqn{\lambda^j_h = 0} for \eqn{j = ( , *)}. Row 1 of
#' the tibble is \eqn{h = 1}, row 2 is \eqn{h = 2} and so on. For estimating
#' \eqn{\lambda^*_h}, put \code{core} as the first argument of the function
#' and headline inflation as the second.
#'
#' @seealso \code{\link{dyn_adj_est}, \link{dyn_adj_best},
#' \link{dyn_adj_pred}}
#'
#' @export
#'
#' @examples
#' inf_head <- coreinf_br[["ipca"]]
#' inf_corems <- coreinf_br[["ipcams"]]
#' inf_coredp <- coreinf_br[["ipcadp"]]
#' dyn_adj(inf_head, inf_corems, 4, 5, recursive = TRUE)
#' out <- purrr::map(list(inf_corems, inf_coredp), ~ dyn_adj(inf_head, .x, 4, 5, recursive = TRUE))
#' purrr::map(out, dplyr::summarise_all, mean)
dyn_adj <- function(inf, core, H, p, ...) {
  purrr::map_df(1:H, ~
           dyn_adj_pred(
             dyn_adj_best(
               dyn_adj_est(inf, core, .x, p, ...)
             )
           )
  )
}