rjdqa
Quality Assessment for Seasonal Adjustment

Package index
Search the rjdqa package
Vignettes
Functions
42
Source code
10
Man pages
4
Home
/
CRAN
/
rjdqa
/
R/sc_d_net_effects.R

R/sc_d_net_effects.R
In rjdqa: Quality Assessment for Seasonal Adjustment

Defines functions plot_net_graph plot_net_effect_text plot.net_effect net_effect.TRAMO_SEATS net_effect.X13 net_effect 

net_effect <- function(x){
    UseMethod("net_effect", x)
}
#'@exportS3Method NULL
net_effect.X13 <- function(x){
    sa <- x$final$series[,"sa"]
    y <- x$final$series[,"y"]
    s <- x$final$series[,"s"]
    
    tde <- x$regarima$model$effects[,"tde"]
    mhe <- x$regarima$model$effects[,"ee"]
    is_multiplicative <- x$regarima$model$spec_rslt[,"Log transformation"]
    # if (is_multiplicative) {
    #     tde <- exp(tde)
    #     mhe <- exp(mhe)
    # }
    
    observed_sa <- tail(sa, 2)
    observed_y <- tail(y, 2)
    
    observed_tde <- tail(tde, 2)
    observed_mhe <- tail(mhe, 2)
    
    observed_sa_evol <- (observed_sa[2] / observed_sa[1] - 1) * 100
    observed_y_evol <- (observed_y[2] / observed_y[1] - 1) * 100
    
    neutral_line_sa <- rep(observed_sa[1], 2)
    
    if (is_multiplicative) {
        s <- s / (exp(tde + mhe))
        observed_td <- observed_y / exp(observed_tde)
        observed_cal <- observed_y / exp(observed_tde + observed_mhe)
        
        neutral_line_y <- observed_y * c(1, observed_sa[1] / observed_sa[2])
        neutral_line_td <- observed_td * c(1, observed_sa[1] / observed_sa[2])
        neutral_line_cal <- observed_cal * c(1, observed_sa[1] / observed_sa[2])
    }else{
        s <- s - tde - mhe
        observed_td <- observed_y - observed_tde
        observed_cal <- observed_y - observed_tde - observed_mhe
        
        neutral_line_y <- observed_y + c(0, observed_sa[1] - observed_sa[2])
        neutral_line_td <- observed_td + c(0, observed_sa[1] - observed_sa[2])
        neutral_line_cal <- observed_cal + c(0, observed_sa[1] - observed_sa[2])
    }
    observed_s <- tail(s, 2)
    observed_s_evol <- (observed_s[2] / observed_s[1] - 1) * 100
    
    # if (is_multiplicative) {
    #     s <- s / (tde * mhe)
    #     observed_s <- tail(s, 2)
    #     observed_s_evol <- (observed_s[2] / observed_s[1] - 1) * 100
    #     
    #     observed_td <- observed_y / observed_tde
    #     observed_cal <- observed_y / (observed_tde * observed_mhe)
    #     
    #     neutral_line_y <- observed_y * c(1, observed_sa[1] / observed_sa[2])
    #     neutral_line_td <- observed_td * c(1, observed_sa[1] / observed_sa[2])
    #     neutral_line_cal <- observed_cal * c(1, observed_sa[1] / observed_sa[2])
    #     
    # }else{
    #     s <- s - tde - mhe
    #     observed_s <- tail(s, 2)
    #     observed_s_evol <- (observed_s[2] / observed_s[1] - 1) * 100
    #     
    #     observed_td <- observed_y - observed_tde
    #     observed_cal <- observed_y - observed_tde - observed_mhe
    #     
    #     neutral_line_y <- observed_y + c(0, observed_sa[1] - observed_sa[2])
    #     neutral_line_td <- observed_td + c(0, observed_sa[1] - observed_sa[2])
    #     neutral_line_cal <- observed_cal + c(0, observed_sa[1] - observed_sa[2])
    # }
    
    EM <- sd(tail(ev(sa) / 100, 5 * frequency(sa)), na.rm = TRUE)
    
    upper_bound_sa <- neutral_line_sa * c(1, 1 + EM)
    lower_bound_sa <- neutral_line_sa * c(1, 1 - EM)
    upper_bound_y <- neutral_line_y * c(1, 1 + EM)
    lower_bound_y <- neutral_line_y * c(1, 1 - EM)
    
    y_movement <- c("decrease", "increase")[(observed_y_evol > 0) + 1]
    
    neutral_line_y_evol <- (neutral_line_y[2] / neutral_line_y[1] - 1) * 100
    neutral_line_td_evol <- (neutral_line_td[2] / neutral_line_td[1] - 1) * 100
    neutral_line_cal_evol <- (neutral_line_cal[2] / neutral_line_cal[1] - 1) * 100
    
    
    tde_contrib <- neutral_line_y_evol - neutral_line_td_evol
    mhe_contrib <- neutral_line_td_evol - neutral_line_cal_evol
    observed_s_evol <- neutral_line_cal_evol
    # if(!is_multiplicative){
    #     # observed_sa + observed_mhe + observed_tde + observed_s - observed_y
    #     # neutral_line_sa + observed_mhe + observed_tde + observed_s - neutral_line_y
    #     # observed_mhe + observed_tde + observed_s + observed_sa
    #     # observed_y
    #     # observed_td <- neutral_line_y - neutral_line_td
    #     # observed_mhe <- neutral_line_td - neutral_line_cal
    #     # observed_td + observed_mhe + observed_s_evol + observed_sa[1]
    #     # neutral_line_y
    #     # tde_contrib <- (observed_tde[2]/ observed_tde[1] - 1) * observed_tde[1] / observed_y[1] * 100
    #     # mhe_contrib <- (observed_mhe[2]/ observed_mhe[1] - 1) * observed_mhe[1] / observed_y[1] * 100
    #     # observed_s_evol <- observed_s_evol * observed_s[1] / observed_y[1]
    #     if(is.nan(mhe_contrib))
    #         mhe_contrib <- observed_mhe[2] / observed_y[1] * 100
    #     if(is.nan(tde_contrib))
    #         tde_contrib <- tde_contrib[2] / observed_y[1] * 100
    #     if(is.nan(observed_s_evol))
    #         observed_s_evol <- observed_s[2] / observed_y[1] * 100
    #     # c(tde_contrib,
    #     # mhe_contrib,
    #     # observed_s_evol)
    #     # tde_contrib + mhe_contrib + observed_s_evol
    #     # tde_contrib + observed_s_evol
    #     #
    #     # neutral_line_y_evol
    # }
 
    
    neutral_y_movement <- c("decrease", "increase")[(neutral_line_y_evol > 0) + 1]
    s_movement <- c("decrease", "increase")[(observed_s_evol > 0) + 1]
    
    neutral_cal_movement <- c("decrease", "increase")[(neutral_line_cal_evol > 0) + 1]
    
    line1 <- sprintf("Observed %.1f%% raw %s from last month",
                     abs(observed_y_evol), y_movement)
    line2 <- sprintf("Neutral result requires %.1f%% raw %s",
                     abs(neutral_line_y_evol), neutral_y_movement)
    line2 <- c(line2, "from last month:")
    if(round(tde_contrib, 1) == 0){
        line3 <- "No trading day effect"
    }else{
        tde_contrib_movement <- c("decrease", "increase")[(tde_contrib > 0) + 1]
        line3 <- sprintf("Trading day effects represent %.1f%% raw %s",
                         abs(tde_contrib), tde_contrib_movement)
    }
    if (round(mhe_contrib, 1) == 0) {
        line4 <- "No moving holiday effect"
    }else{
        mhe_contrib_movement <- c("decrease", "increase")[(mhe_contrib > 0) + 1]
        line4 <- sprintf("Moving holiday effects represent %.1f%% raw %s",
                         abs(mhe_contrib), mhe_contrib_movement)
    }
    line5 <- sprintf("Seasonal effects represent %.1f%% raw %s",
                     abs(observed_s_evol), s_movement)
    line6 <- sprintf("SA movement of %+.1f%% from last month",
                     observed_sa_evol)
    res <- list(raw_data = list(observed = observed_y, neutral_line = neutral_line_y,
                                lower_bound = lower_bound_y, upper_bound = upper_bound_y,
                                frequency = frequency(y)),
                sa_data = list(observed = observed_sa, neutral_line = neutral_line_sa,
                               lower_bound = lower_bound_sa, upper_bound = upper_bound_sa,
                               frequency = frequency(sa)),
                text = c(line1, line2,
                         line3, line4,
                         line5, line6)
    )
    class(res) <- c("net_effect", class(res))
    
    res
}
#'@exportS3Method NULL
net_effect.TRAMO_SEATS <- function(x){
    sa <- x$final$series[,"sa"]
    y <- x$final$series[,"y"]
    s <- x$final$series[,"s"]
    
    tde <- x$regarima$model$effects[,"tde"]
    mhe <- x$regarima$model$effects[,"ee"]
    is_multiplicative <- x$regarima$model$spec_rslt[,"Log transformation"]
    
    observed_sa <- tail(sa, 2)
    observed_y <- tail(y, 2)
    
    if (is_multiplicative) {
        tde <- exp(tde)
        mhe <- exp(mhe)
    }
    
    observed_tde <- tail(tde, 2)
    observed_mhe <- tail(mhe, 2)
    
    observed_sa_evol <- (observed_sa[2] / observed_sa[1] - 1) * 100
    observed_y_evol <- (observed_y[2] / observed_y[1] - 1) * 100
    
    neutral_line_sa <- rep(observed_sa[1], 2)
    if (is_multiplicative) {
        s <- s / (tde * mhe)
        observed_s <- tail(s, 2)
        observed_s_evol <- (observed_s[2] / observed_s[1] - 1) * 100
        
        observed_td <- observed_y / observed_tde
        observed_cal <- observed_y / (observed_tde * observed_mhe)
        
        neutral_line_y <- observed_y * c(1, observed_sa[1] / observed_sa[2])
        neutral_line_td <- observed_td * c(1, observed_sa[1] / observed_sa[2])
        neutral_line_cal <- observed_cal * c(1, observed_sa[1] / observed_sa[2])
    }else{
        s <- s - tde - mhe
        observed_s <- tail(s, 2)
        observed_s_evol <- (observed_s[2] / observed_s[1] - 1) * 100
        
        observed_td <- observed_y - observed_tde
        observed_cal <- observed_y - observed_tde - observed_mhe
        
        neutral_line_y <- observed_y + c(0, observed_sa[1] - observed_sa[2])
        neutral_line_td <- observed_td + c(0, observed_sa[1] - observed_sa[2])
        neutral_line_cal <- observed_cal + c(0, observed_sa[1] - observed_sa[2])
    }
    
    EM <- sd(tail(ev.ts(sa) / 100, 5 * frequency(sa)), na.rm = TRUE)
    
    upper_bound_sa <- neutral_line_sa * c(1, 1 + EM)
    lower_bound_sa <- neutral_line_sa * c(1, 1 - EM)
    upper_bound_y <- neutral_line_y * c(1, 1 + EM)
    lower_bound_y <- neutral_line_y * c(1, 1 - EM)
    
    y_movement <- c("decrease", "increase")[(observed_y_evol > 0) + 1]
    
    neutral_line_y_evol <- (neutral_line_y[2] / neutral_line_y[1] - 1) * 100
    neutral_line_td_evol <- (neutral_line_td[2] / neutral_line_td[1] - 1) * 100
    neutral_line_cal_evol <- (neutral_line_cal[2] / neutral_line_cal[1] - 1) * 100
    
    
    tde_contrib <- neutral_line_y_evol - neutral_line_td_evol
    mhe_contrib <- neutral_line_td_evol - neutral_line_cal_evol
    observed_s_evol <- neutral_line_cal_evol
    
    
    neutral_y_movement <- c("decrease", "increase")[(neutral_line_y_evol > 0) + 1]
    s_movement <- c("decrease", "increase")[(observed_s_evol > 0) + 1]
    
    neutral_cal_movement <- c("decrease", "increase")[(neutral_line_cal_evol > 0) + 1]
    
    line1 <- sprintf("Observed %.1f%% raw %s from last month",
                     abs(observed_y_evol), y_movement)
    line2 <- sprintf("Neutral result requires %.1f%% raw %s",
                     abs(neutral_line_y_evol), neutral_y_movement)
    line2 <- c(line2, "from last month:")
    if(round(tde_contrib, 1) == 0){
        line3 <- "No trading day effect"
    }else{
        tde_contrib_movement <- c("decrease", "increase")[(tde_contrib > 0) + 1]
        line3 <- sprintf("Trading day effects represent %.1f%% raw %s",
                         abs(tde_contrib), tde_contrib_movement)
    }
    if (round(mhe_contrib, 1) == 0) {
        line4 <- "No moving holiday effect"
    }else{
        mhe_contrib_movement <- c("decrease", "increase")[(mhe_contrib > 0) + 1]
        line4 <- sprintf("Moving holiday effects represent %.1f%% raw %s",
                         abs(mhe_contrib), mhe_contrib_movement)
    }
    line5 <- sprintf("Seasonal effects represent %.1f%% raw %s",
                     abs(observed_s_evol), s_movement)
    line6 <- sprintf("SA movement of %+.1f%% from last month",
                     observed_sa_evol)
    res <- list(raw_data = list(observed = observed_y, neutral_line = neutral_line_y,
                                lower_bound = lower_bound_y, upper_bound = upper_bound_y,
                                frequency = frequency(y)),
                sa_data = list(observed = observed_sa, neutral_line = neutral_line_sa,
                               lower_bound = lower_bound_sa, upper_bound = upper_bound_sa,
                               frequency = frequency(sa)),
                text = c(line1, line2,
                         line3, line4,
                         line5, line6)
    )
    class(res) <- c("net_effect", class(res))
    
    res
}
# net_effect <- function(x){
#     sa <- x$final$series[,"sa"]
#     y <- x$final$series[,"y"]
#     s <- x$final$series[,"s"]
#     
#     tde <- x$regarima$model$effects[,"tde"]
#     mhe <- x$regarima$model$effects[,"ee"]
#     is_multiplicative <- x$regarima$model$spec_rslt["Log transformation"] == TRUE
#     
#     observed_sa <- tail(sa, 2)
#     observed_y <- tail(y, 2)
#     
#     observed_tde <- tail(tde, 2)
#     observed_mhe <- tail(mhe, 2)
#     
#     observed_sa_evol <- (observed_sa[2]/ observed_sa[1] - 1) * 100
#     observed_y_evol <- (observed_y[2]/ observed_y[1] - 1) * 100
#     tde_evol <- (observed_tde[2]/ observed_tde[1] - 1) * 100
#     mhe_evol <- (observed_mhe[2]/ observed_mhe[1] - 1) *100
#     # observed_s_evol <- observed_s_evol * observed_s[1] / observed_y[1]
#     
#     neutral_line_sa <- rep(observed_sa[1], 2)
#     if(is_multiplicative){
#         s <- s / (tde * mhe)
#         observed_s <- tail(s, 2)
#         observed_s_evol <- (observed_s[2]/ observed_s[1] - 1) * 100
#         observed_s_contrib <-  observed_s_evol
#         
#         neutral_line_y <- observed_y * c(1, observed_sa[1] / observed_sa[2])
#         tde_contrib <- tde_evol
#         mhe_contrib <- mhe_evol
#     }else{
#         s <- s - tde - mhe
#         observed_s <- tail(s, 2)
#         observed_s_evol <- (observed_s[2]/ observed_s[1] - 1) * 100
#         observed_s_contrib <-  observed_s_evol * observed_s[1] / observed_y[1]
#         
#         neutral_line_y <- observed_y + c(0, observed_sa[1] - observed_sa[2])
#         tde_contrib <- tde_evol * observed_tde[1] / observed_y[1] * 100
#         mhe_contrib <- mhe_evol * observed_mhe[1] / observed_y[1] * 100
#         observed_s_evol <- observed_s_evol * observed_s[1] / observed_y[1]
#         if(is.nan(mhe_contrib))
#             mhe_contrib <- observed_mhe[2] / observed_y[1] * 100
#         if(is.nan(tde_contrib))
#             tde_contrib <- tde_contrib[2] / observed_y[1] * 100
#         if(is.nan(observed_s_evol))
#             observed_s_evol <- observed_s[2] / observed_y[1] * 100
#     }
#     
#     EM <- sd(tail(ev(sa) / 100, 5 * frequency(sa)), na.rm = TRUE)
#     
#     upper_bound_sa <- neutral_line_sa * c(1, 1 + EM)
#     lower_bound_sa <- neutral_line_sa * c(1, 1 - EM)
#     upper_bound_y <- neutral_line_y * c(1, 1 + EM)
#     lower_bound_y <- neutral_line_y * c(1, 1 - EM)
#     
#     y_movement <- c("decrease", "increase")[(observed_y_evol > 0) + 1]
#     
#     neutral_line_y_evol <- (neutral_line_y[2]/ neutral_line_y[1] - 1) * 100
#     
#     
#     neutral_y_movement <- c("decrease", "increase")[(neutral_line_y_evol > 0) + 1]
#     s_movement <- c("decrease", "increase")[(observed_s_evol > 0) + 1]
#     
# 
#     line1 <- sprintf("Observed %.1f%% raw %s from last month",
#                      abs(observed_y_evol), y_movement)
#     line2 <- sprintf("Neutral result requires %.1f%% raw %s from last month",
#                      abs(neutral_line_y_evol), neutral_y_movement)
#     if(round(tde_contrib,1) == 0){
#         line3 <- "No trading day effect"
#     }else{
#         tde_contrib_movement <- c("decrease", "increase")[(tde_contrib > 0) + 1]
#         line3 <- sprintf("Trading day effects represent %.1f%% raw %s",
#                          abs(tde_contrib), tde_contrib_movement)
#     }
#     if(round(mhe_contrib,1) == 0){
#         line4 <- "No moving holiday effect"
#     }else{
#         mhe_contrib_movement <- c("decrease", "increase")[(mhe_contrib > 0) + 1]
#         line4 <- sprintf("Moving holiday effects represent %.1f%% raw %s",
#                          abs(mhe_contrib), mhe_contrib_movement)
#     }
#     line5 <- sprintf("Seasonal effects represent %.1f%% raw %s",
#                      abs(observed_s_evol), s_movement)
#     line6 <- sprintf("SA movement of %+.1f%% from last month",
#                      observed_sa_evol)
#     res <- list(raw_data = list(observed = observed_y, neutral_line = neutral_line_y,
#                                 lower_bound = lower_bound_y, upper_bound = upper_bound_y,
#                                 frequency = frequency(y)),
#                 sa_data = list(observed = observed_sa, neutral_line = neutral_line_sa,
#                                lower_bound = lower_bound_sa, upper_bound = upper_bound_sa,
#                                frequency = frequency(sa)),
#                 text = c(line1, line2,
#                          line3, line4,
#                          line5, line6)
#     )
#     class(res) <- c("net_effect", class(res))
#     
#     res
# }
#'@exportS3Method NULL
plot.net_effect <- function(x,
                            raw_color = "#33A02C",
                            sa_color = "#E31A1C", ...){
    mai <- par('mai')
    
    par(mai = c(mai[1], mai[2], mai[3], 0.1))
    plot_net_graph(x$raw_data, title = "Unadjusted", subtitle = "Raw",
                   arrow_col = raw_color)
    par(mai = c(mai[1], 0, 0.1, 0.1))
    plot_net_effect_text(x$text)
    par(mai = c(mai[1], 0.2, mai[3], mai[4]))
    
    plot_net_graph(x$sa_data, title = "Seasonally Adjusted", subtitle = "SA",
                   arrow_col = sa_color)
    par(mai = mai)
}
plot_net_effect_text <- function(x){
    base_inset_x <- -0.00
    plot.new()
    box()
    legend("topleft", legend = x[1], 
           bty = "n", text.font =  2, inset = c(base_inset_x, 0))
    legend("topleft", legend = c(rep(NA, 3), x[2:3]), 
           bty = "n", text.font =  2, inset = c(base_inset_x, 0))
    legend("topleft", legend = c(rep(NA, 5), x[4:6]), 
           bty = "n", text.font =  1,inset = c(base_inset_x + 0.05,0))
    legend("topleft", legend = c(rep(NA, 10), x[7]), 
           bty = "n", text.font =  2, inset = c(base_inset_x, 0))
}
plot_net_graph <- function(x, title = "", subtitle = "", arrow_col = "black"){
    
    main <- sprintf("%s (%s)", title, subtitle)
    legend_text <- c(sprintf("%s, %+.1f%%",subtitle,
                             (x$observed[2]/ x$observed[1] - 1)*100), 
                     sprintf("%s %s",subtitle, c("Neutral", "Boundaries")))
    if(x$frequency == 12){
        x_lab <- c("m-1","m")
    }else{
        x_lab <- c("Q-1","Q")
    }
    plot(x = 1:2, y = x$neutral_line, type = "l",
         ylim = range(c(x$upper_bound, x$lower_bound)),
         lty = c("dotted"),xlim = c(0.8,2.1),
         xlab = "", ylab = "", main = main,
         xaxt = "n")
    lines(x = 1:2, y = x$upper_bound, lty = c("dashed"),lwd = 2)
    lines(x = 1:2, y = x$lower_bound, lty = c("dashed"),lwd = 2)
    arrows(1, x$observed[1], 2, x$observed[2], col = arrow_col,lwd = 2)
    axis(1, at = 1:2,
         labels = x_lab,
         las = 0)
    legend("topleft", legend = legend_text, pch = c(NA, NA,NA),
           lty = c("solid","dotted","dashed"), col = c(arrow_col, "black","black"),
           lwd = 3,cex = 0.8, bty="n") 
}

Try the rjdqa package in your browser

Any scripts or data that you put into this service are public.

rjdqa documentation built on June 22, 2024, 10:20 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close