ur_forecast_error.R

# Estimate confidence intervals around RBA forecasts
# As per https://www.rba.gov.au/publications/rdp/2012/pdf/rdp2012-07.pdf
library(tidyverse)
library(readabs)
library(readrba)
library(tsibble)

raw_forecasts <- read_forecasts()

raw_actual_ur <- read_abs_series("A84423050A")

actual_ur <- raw_actual_ur |>
  mutate(year_qtr = yearquarter(date)) |>
  group_by(year_qtr) |>
  summarise(actual_value = mean(value))

forecast_ur <- raw_forecasts |>
  filter(series == "unemp_rate")  |>
  mutate(forecast_qtr = yearquarter(forecast_date),
         year_qtr = yearquarter(date)) |>
  mutate(horizon = year_qtr - forecast_qtr) |>
  select(year_qtr,
         forecast_qtr,
         forecast_value = value,
         horizon) |>
  filter(horizon >= 0,
         year(forecast_qtr) >= 1993)

ur_errors <- forecast_ur |>
  left_join(actual_ur,
            by = "year_qtr") |>
  filter(!is.na(actual_value)) |>
  mutate(forecast_error = actual_value - forecast_value)

ur_errors |>
  filter(horizon <= 8) |>
  ggplot(aes(x = forecast_error,
             y = horizon,
             fill = factor(horizon),
             group = horizon)) +
  ggridges::geom_density_ridges(#fill = "light blue",
                                quantile_lines = TRUE,
                                quantiles = 0.5,
                                rel_min_height = 0.005) +
  scale_y_continuous(breaks = 0:8) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none") +
  labs(y = "Forecast horizon (quarters)",
       x = "Forecast error (percentage points)")

ur_errors |>
  filter(year(year_qtr) >= 1993,
         horizon <= 8) |>
  mutate(abs_forecast_error = abs(forecast_error)) |>
  group_by(horizon) |>
  summarise(ave_error = list(quantile(abs_forecast_error,
                                  probs = c(0.5, 0.7, 0.9))),
            .groups = "keep")  |>
  unnest(ave_error) |>
  mutate(q = factor(c("Median",
               "70th percentile",
               "90th percentile"),
               levels = c("Median",
                          "70th percentile",
                          "90th percentile"))) |>
  ggplot(aes(x = horizon,
             y = ave_error,
             fill = q)) +
  geom_col(position = "dodge") +
  coord_flip() +
  theme_minimal(base_size = 14) +
  theme(legend.title = element_blank(),
       legend.position = c(0.8, 0.1)) +
  labs(x = "Forecast horizon (quarters)",
       y = "Absolute forecast error (percentage points)")