R/gen_air_components.r

In bcmonsell/airutilities: Support routines for fractional airline decomposition routines in the rjd3highfreq package

#' Generate signal extraction components
#'
#' Generates components from the decomposition of the fractional airline model R routines, partioning regression effects to different components as in \code{X-13ARIMA-SEATS}
#'
#' @param this_est A list object generated by the estimation procedure of the fractional airline model
#' @param this_decomp A list object generated by the decomposition procedure of the fractional airline model
#' @param this_xtype Character vector; Type of regressor used for including regression effects into components (\code{'hol'} for holiday, \code{'ao'} for point outliers, \code{'ls'} for level change outliers).
#' @param this_x Character vector; regression. Default: use the regression matrix in \code{this_est$model$X}.
#' @param this_log Logical scalar; set to TRUE if the log transformation is used in the decomposition, FALSE otherwise. Default: TRUE
#' @param this_stde Logical scalar; component standard error included in \code{this_decomp}. Default: FALSE.
#' @param tc2trend Logical scalar; TC outliers included in trend component if TRUE. Default: FALSE, TC outliers included in irregular.
#' @param remove_user Logical scalar; remove user-defined regression factors from seasonally adjusted series if TRUE. Default: FALSE, user defined regression factors included in seasonally adjusted series.
#' @return List of components from the fractional airline decomposition, partioning regression effects to different components as in \code{X-13ARIMA-SEATS}
#' @examples
#' ic_est <-
#'    rjd3highfreq::fractionalAirlineEstimation(ic_obs, periods=c(365.25/7),
#'                                       x=ic_default_matrix)
#' ic_decomp <-
#'    rjd3highfreq::fractionalAirlineDecomposition(ic_est$model$linearized,
#'                                          365.25/7, stde = TRUE)
#' ic_xtype <- c(rep('hol', 13), rep('ao', 40))
#' ic_comp <- 
#'    gen_air_components(ic_est, ic_decomp,
#'                       this_xtype = ic_xtype,
#'                       this_log = FALSE, 
#'                       this_stde = TRUE)
#' ic_comp_tis <- lapply(ic_comp, function(x)
#'    try(tis::tis(x, start = ic_start, tif = 'wsaturday')))
#' @import stats
#' @export
gen_air_components <- function(this_est, this_decomp, this_xtype, this_x = NULL,
                               this_log = TRUE, this_stde = FALSE, tc2trend = FALSE, 
                               remove_user = FALSE) {
    # Author: Brian C. Monsell (OEUS), Version 4.0, 2/2/2022
    # set variables for the regression matrix and coefficient estimates
    if (is.null(this_x)) {
      this_x <- this_est$model$X
    }
    this_b <- this_est$model$b
    
    # generate filters for each type of regressors
    filter_hol <- this_xtype == "hol"
    filter_ao <- this_xtype == "ao"
    filter_ls <- this_xtype == "ls"
    filter_tc <- this_xtype == "tc"
    filter_user <- this_xtype == "user"
    
    n_hol <- sum(filter_hol)
    n_ao <- sum(filter_ao)
    n_ls <- sum(filter_ls)
    n_tc <- sum(filter_tc)
    n_user <- sum(filter_user)
    
    # generate default factor based on log vs no log
    this_y <- this_est$model$y
    this_length <- length(this_est$model$y)
    if (length(this_decomp$decomposition$y) > this_length) {
        this_y <- 
            c(this_est$model$y, 
              this_decomp$decomposition$y[(this_length+1):length(this_decomp$decomposition$y)])
        this_length <- length(this_y)
    }

    if (this_log) {
        this_default <- rep(1, length.out = this_length)
    } else {
        this_default <- rep(0, length.out = this_length)
    }
    
    # generate holiday factors
    if (n_hol > 0) {
        this_hol <- as.vector(this_x[, filter_hol] %*% matrix(this_b[filter_hol], ncol = 1))
        if (this_log) {
            this_hol <- exp(this_hol)
        }
    } else {
        this_hol <- this_default
    }
    
    # generate ao factors
    if (n_ao > 0) {
        this_ao <- 
           as.vector(this_x[, filter_ao] %*% matrix(this_b[filter_ao], ncol = 1))
        if (this_log) {
            this_ao <- exp(this_ao)
        }
    } else {
        this_ao <- this_default
    }
    
    # generate ls factors
    if (n_ls > 0) {
        this_ls <- 
           as.vector(this_x[, filter_ls] %*% matrix(this_b[filter_ls], ncol = 1))
        if (this_log) {
            this_ls <- exp(this_ls)
        }
    } else {
        this_ls <- this_default
    }
    
    # generate tc factors
    if (n_tc > 0) {
        this_tc <- 
           as.vector(this_x[, filter_tc] %*% matrix(this_b[filter_tc], ncol = 1))
        if (this_log) {
            this_tc <- exp(this_tc)
        }
    } else {
        this_tc <- this_default
    }
    
    this_outlier <- this_default
    if (this_log) {
      if (n_ao > 0) { this_outlier <- this_outlier * this_ao }
      if (n_ls > 0) { this_outlier <- this_outlier * this_ls }
      if (n_tc > 0) { this_outlier <- this_outlier * this_tc }
    } else {
      if (n_ao > 0) { this_outlier <- this_outlier + this_ao }
      if (n_ls > 0) { this_outlier <- this_outlier + this_ls }
      if (n_tc > 0) { this_outlier <- this_outlier + this_tc }
    }
    
    # generate user factors
    if (n_user > 0) {
        this_user <- 
           as.vector(this_x[, filter_user] %*% matrix(this_b[filter_user], ncol = 1))
        if (this_log) {
            this_user <- exp(this_user)
        }
    } else {
        this_user <- this_default
    }
    
    # generate trend component
    this_trend <- this_decomp$decomposition$t
    if (this_log) {
        this_trend <- exp(this_trend) * this_ls
    } else {
        this_trend <- this_trend + this_ls
    }
    
    # generate irregular component
    this_irr <- this_decomp$decomposition$i
    if (this_log) {
        this_irr <- exp(this_irr) * this_ao
    } else {
        this_irr <- this_irr + this_ao
    }
    
    # assign tc outliers to a component based on \code{tc2trend}
    if (tc2trend) {
      if (this_log) {
        this_trend <- exp(this_trend) * this_tc
      } else {
        this_trend <- this_trend + this_tc
      }
    } else {
      if (this_log) {
        this_irr <- exp(this_irr) * this_tc
      } else {
        this_irr <- this_irr + this_tc
      }
    }
    
    # generate seasonally adjusted series
    this_sa <- this_decomp$decomposition$sa
    if (this_log) {
        this_sa <- exp(this_sa) * this_ao * this_ls * this_tc
        if (remove_user) this_sa <- this_sa / this_user
    } else {
        this_sa <- this_sa + this_ao + this_ls + this_tc
        if (remove_user) this_sa <- this_sa - this_user
    }
    
    
    # generate seasonal and combinded factor component
    this_seasonal <- this_decomp$decomposition$s
    if (this_log) {
        this_seasonal <- exp(this_seasonal)
        this_combfac <- this_seasonal * this_hol
    } else {
        this_combfac <- this_seasonal + this_hol
    }
    
    # generate a list of components to return
    if (this_stde) {
        this_return <- list(y = this_est$model$y, trend = this_trend, seasonal = this_seasonal, irr = this_irr, 
            adjfac = this_combfac, sa = this_sa, holiday = this_hol, ao = this_ao, ls = this_ls, tc = this_tc, outlier = this_outlier,
            user = this_user, trend_stde = this_decomp$decomposition$t.stde, seasonal_stde = this_decomp$decomposition$s.stde, 
            irr_stde = this_decomp$decomposition$i.stde)
    } else {
        this_return <- list(y = this_est$model$y, trend = this_trend, seasonal = this_seasonal, irr = this_irr, 
            adjfac = this_combfac, sa = this_sa, holiday = this_hol, ao = this_ao, ls = this_ls, tc = this_tc, outlier = this_outlier,
            user = this_user)
    }
    
    # return list of components
    return(this_return)
    
}