## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 3.5
)
set.seed(2024)
library(bridgr)

## ----bootstrap-model----------------------------------------------------------
gdp_growth <- suppressMessages(tsbox::ts_na_omit(tsbox::ts_pc(gdp)))

boot_model <- mf_model(
  target = gdp_growth,
  indic = baro,
  indic_predict = "auto.arima",
  indic_aggregators = "mean",
  indic_lags = 1,
  target_lags = 1,
  h = 2,
  se = TRUE,
  bootstrap = list(N = 20, block_length = NULL)
)

## ----forecast-output----------------------------------------------------------
fc <- forecast(boot_model)

fc
fc$bootstrap

## ----width-tolerance----------------------------------------------------------
forecast_table <- dplyr::tibble(
  time = fc$time,
  mean = as.numeric(fc$mean),
  lower_95 = fc$lower[, "95%"],
  upper_95 = fc$upper[, "95%"]
) |>
  dplyr::mutate(
    half_width_95 = (.data$upper_95 - .data$lower_95) / 2
  )

tolerance <- 1.5
dplyr::filter(forecast_table, .data$half_width_95 <= tolerance)

## ----summary-output-----------------------------------------------------------
summary(boot_model)

## ----scenarios----------------------------------------------------------------
# `model.frame(which = "forecast")` returns the regressor path the baseline
# forecast is built from; shifting it is all a scenario needs.
baseline_xreg <- model.frame(boot_model, which = "forecast")
xreg_names <- variable.names(boot_model, which = "xreg")

make_xreg <- function(level_shift) {
  dplyr::tibble(
    id = rep(xreg_names, each = nrow(baseline_xreg)),
    time = rep(baseline_xreg$time, times = length(xreg_names)),
    value = c(
      baseline_xreg$baro + level_shift,
      baseline_xreg$baro_lag1 + level_shift
    )
  )
}

fc_baseline <- forecast(boot_model, xreg = make_xreg(0))
fc_shock    <- forecast(boot_model, xreg = make_xreg(-5))

scenario_df <- dplyr::bind_rows(
  dplyr::tibble(
    scenario = "baseline",
    time = fc_baseline$time,
    mean = as.numeric(fc_baseline$mean),
    lower = fc_baseline$lower[, "95%"],
    upper = fc_baseline$upper[, "95%"]
  ),
  dplyr::tibble(
    scenario = "shock (-5)",
    time = fc_shock$time,
    mean = as.numeric(fc_shock$mean),
    lower = fc_shock$lower[, "95%"],
    upper = fc_shock$upper[, "95%"]
  )
)

ggplot2::ggplot(
  scenario_df,
  ggplot2::aes(x = .data$time, color = .data$scenario, fill = .data$scenario)
) +
  ggplot2::geom_ribbon(
    ggplot2::aes(ymin = .data$lower, ymax = .data$upper),
    alpha = 0.2, color = NA
  ) +
  ggplot2::geom_line(ggplot2::aes(y = .data$mean), linewidth = 0.8) +
  ggplot2::geom_point(ggplot2::aes(y = .data$mean), size = 2) +
  ggplot2::labs(
    title = "Forecast under baseline and shock scenarios",
    x = NULL, y = "GDP growth forecast"
  ) +
  theme_bridgr()

## ----full-system-bootstrap----------------------------------------------------
full_model <- mf_model(
  target = gdp_growth,
  indic = baro,
  indic_predict = "auto.arima",
  indic_aggregators = "mean",
  indic_lags = 1,
  target_lags = 1,
  h = 2,
  se = TRUE,
  full_system_bootstrap = TRUE,
  bootstrap = list(N = 20, block_length = NULL)
)

forecast(full_model)$bootstrap

## ----point-only---------------------------------------------------------------
point_model <- mf_model(
  target = gdp_growth,
  indic = baro,
  indic_predict = "auto.arima",
  indic_aggregators = "mean",
  indic_lags = 1,
  target_lags = 1,
  h = 1,
  se = FALSE,
  bootstrap = list(N = 20)
)

forecast(point_model)