R/gen_air_projected_factors.r

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

#' Generate projected seasonal (and holiday) adjustment factors
#'
#' Generates projected adjustment factors from the decomposition of the fractional airline model R routines
#'
#' @param this_est A list object generated by the estimation procedure of the fractional airline model; the regression matrix should not be extended with forecasts.
#' @param this_decomp A list object generated by the decomposition procedure of the fractional airline model; the \code{nfcasts} option should have been specified. 
#' @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 matrix. 
#' @param this_log Logical scalar; set to TRUE if the log transformation is used in the decomposition, FALSE otherwise. Default: TRUE
#' @param nfcasts Integer scalar; number of projected adjustment factors; should match the \code{nfcasts} argument used to generate the decomposition. 
#' @param return_series Character scalar; component returned by the routine. Default is "combined"; other possible choices are "seasonal" and "holiday".
#' @param return_xts Logical scalar; return projected factors of type \code{return_series} as an \code{xts} time series object. Default is FALSE.
#'        If TRUE, \code{this_x} must be an \code{xts} time series object.
#' @return Array of projected adjustment factors from the fractional airline decomposition
#' @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, nfcasts = 104)
#' ic_xtype <- c(rep('hol', 13), rep('ao', 40))
#' ic_proj_seasonal <- 
#'    gen_air_projected_factors(ic_est, ic_decomp,
#'                       this_xtype = ic_xtype,
#'                       this_x     = ic_holiday_matrix_fcst,
#'                       this_log = FALSE, nfcasts = 104)
#' @import stats
#' @export
gen_air_projected_factors <- function(this_est, this_decomp, this_xtype, this_x = NULL,
                                      this_log = TRUE, nfcasts = 104, return_series = "combined",
                                      return_xts = FALSE) {
    # Author: Brian C. Monsell (OEUS), Version 3.3, 3/16/2022
    test_nfcasts <- length(this_decomp$decomposition$y) - this_decomp$likelihood$nobs
    if (test_nfcasts != nfcasts) {
        stop(paste0("number of forecasts specified for decomposition (",test_nfcasts,") not the same as nfcasts (",nfcasts,")"))
    }
    
    if (!(return_series == "combined" | return_series == "seasonal" | return_series == "holiday")) {
        stop("Acceptable entries for return_series are combined, seasonal, or holiday")
    }
    
    this_max_col <- ncol(this_x)    
    this_b <- this_est$model$b[1:this_max_col]
    
    # generate filters for holiday regressors
    filter_hol <- this_xtype[1:this_max_col] == "hol"
    
    n_hol <- sum(filter_hol)

    this_seasonal        <- this_decomp$decomposition$s     
    this_seasonal_length <- length(this_seasonal)
    
    if (this_log) {
        this_default <- rep(1, length.out = this_seasonal_length)
    } else {
        this_default <- rep(0, length.out = this_seasonal_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
    }
    
    if (this_log) {
        this_seasonal <- exp(this_seasonal)
        this_combfac <- this_seasonal * this_hol
    } else {
        this_combfac <- this_seasonal + this_hol
    }
    
    start_factors <- this_seasonal_length - nfcasts + 1
    if (return_series == "combined") {
        this_proj <- this_combfac[start_factors:this_seasonal_length]
    } else {
        if (return_series == "seasonal") {
            this_proj <- this_seasonal[start_factors:this_seasonal_length]
        } else {
            if (return_series == "holiday") {
                this_proj <- this_hol[start_factors:this_seasonal_length]
            }
        }
    }
    
    # return projected factors
    if (return_xts) {
       this_index <- as.Date(zoo::index(this_x), origin = "1970-01-01")
       this_filter <- seq((this_decomp$likelihood$nobs+1),length(this_decomp$decomposition$y)) 
       this_proj <- 
          xts::xts(x = this_proj, order.by = this_index[this_filter])
       return(this_proj)
    } else {
       return(this_proj)
    }
    
}