JKSE Index Analysis

Libraries

knitr::opts_chunk$set(
    fig.align = "center",
    message = FALSE,
    warning = FALSE,
    out.width = "100%"
)
library(tidyverse)
library(tidyquant)
library(forecast)
library(ggrepel)
library(scales)

# Save a ggplot as a high-resolution PNG and display it inline
hd_plot <- function(p, name, width = 9, height = 5, dpi = 300) {
  dir.create("figures", showWarnings = FALSE)
  path <- file.path("figures", paste0(name, ".png"))
  ggsave(path, plot = p, width = width, height = height, dpi = dpi,
         device = ragg::agg_png, bg = "white")
  knitr::include_graphics(path)
}

# Colour palette and theme
pal <- c(navy = "#1f4e79", orange = "#e07b39", green = "#2e8b57",
         red = "#c0392b", grey = "grey55")

theme_ihsg <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(
      plot.title           = element_text(face = "bold", size = rel(1.25)),
      plot.subtitle        = element_text(colour = "grey40", margin = margin(b = 8)),
      plot.caption         = element_text(colour = "grey50", size = rel(0.75)),
      plot.title.position  = "plot",
      panel.grid.minor     = element_blank(),
      legend.position      = "top",
      legend.justification = "left",
      legend.title         = element_blank(),
      strip.text           = element_text(face = "bold", hjust = 0)
    )
}

source_label <- "Source: Yahoo Finance (^JKSE)"

Settings

forecast_start    <- as.Date("2026-11-02")
forecast_end      <- as.Date("2026-11-10")
train_years       <- 3       # length of the estimation window (years)
run_slow_backtest <- FALSE   # TRUE adds TBATS and NNAR to the backtest (much slower)
lv                <- c(80, 95)

# IDX holidays that could fall inside a forecast horizon (2026 calendar); extend as needed
idx_holidays <- as.Date(c("2026-12-24", "2026-12-25", "2026-12-31"))

# Trading days in the forecast horizon
future_dates <- seq(forecast_start, forecast_end, by = "day")
future_dates <- future_dates[lubridate::wday(future_dates, week_start = 1) <= 5 &
                               !future_dates %in% idx_holidays]
h <- length(future_dates)

future_dates

[1] "2026-11-02" "2026-11-03" "2026-11-04" "2026-11-05" "2026-11-06" "2026-11-09" "2026-11-10"

Data

Daily data are downloaded from Yahoo Finance. Weekend dates and flat bars (open = high = low = close), which Yahoo occasionally records on IDX holidays, are listed first and then removed.

ihsg_raw <- tq_get("^JKSE", from = "2015-01-01")

ihsg_flagged <- ihsg_raw %>%
  dplyr::filter(!is.na(close), date < forecast_start) %>%
  distinct(date, .keep_all = TRUE) %>%
  arrange(date) %>%
  mutate(
    weekend  = lubridate::wday(date, week_start = 1) > 5,
    flat_bar = coalesce(open == high & high == low & low == close, FALSE)
  )

# Rows to be removed
ihsg_flagged %>% filter(weekend | flat_bar)

ihsg <- ihsg_flagged %>%
  filter(!weekend, !flat_bar) %>%
  select(date, close)

tail(ihsg)

Trading Calendar

The series contains IDX trading days only, so gaps longer than a weekend correspond to national holidays and collective leave days.

# Longest gaps between consecutive trading dates
ihsg %>%
  mutate(previous_date = lag(date),
         gap_days      = as.numeric(date - previous_date)) %>%
  filter(gap_days > 3) %>%
  arrange(desc(gap_days)) %>%
  select(previous_date, date, gap_days) %>%
  head(10)

# Trading days per complete calendar year
trading_days <- ihsg %>%
  count(year = lubridate::year(date), name = "trading_days") %>%
  filter(year < lubridate::year(max(ihsg$date)))

trading_days

# Annual seasonal period for TBATS: mean of the last five complete years
annual_period <- round(mean(tail(trading_days$trading_days, 5)))
annual_period

[1] 241

Closing Price

Full history

key_points <- bind_rows(
  ihsg %>% slice_max(close, n = 1, with_ties = FALSE) %>% mutate(tag = "All-time high"),
  ihsg %>% slice_min(close, n = 1, with_ties = FALSE) %>% mutate(tag = "Lowest close"),
  ihsg %>% slice_tail(n = 1)                          %>% mutate(tag = "Latest")
) %>%
  mutate(label = paste0(tag, "\n", format(date, "%d %b %Y"), "\n",
                        comma(close, accuracy = 1)))

p <- ggplot(ihsg, aes(date, close)) +
  geom_line(aes(colour = "Daily close"), linewidth = 0.35) +
  geom_ma(aes(colour = "50-day MA"),  ma_fun = SMA, n = 50,
          linetype = "solid", linewidth = 0.8) +
  geom_ma(aes(colour = "200-day MA"), ma_fun = SMA, n = 200,
          linetype = "solid", linewidth = 0.8) +
  geom_point(data = key_points, size = 2.2) +
  geom_label_repel(data = key_points, aes(label = label),
                   size = 3, label.size = 0.15, fill = alpha("white", 0.9),
                   box.padding = 0.6, min.segment.length = 0, seed = 42) +
  scale_colour_manual(values = c("Daily close" = pal[["grey"]],
                                 "50-day MA"   = pal[["orange"]],
                                 "200-day MA"  = pal[["navy"]]),
                      breaks = c("Daily close", "50-day MA", "200-day MA")) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y",
               expand = expansion(mult = c(0.01, 0.03))) +
  scale_y_continuous(labels = comma, expand = expansion(mult = 0.08)) +
  labs(title = "Jakarta Composite Index (JKSE)",
       subtitle = paste("Daily closing price,",
                        format(min(ihsg$date), "%b %Y"), "to",
                        format(max(ihsg$date), "%b %Y")),
       x = NULL, y = "Index level", caption = source_label) +
  theme_ihsg()

hd_plot(p, "01_close_full")

Last 12 months

ihsg_12m <- ihsg %>% filter(date >= max(date) %m-% months(12))

p <- ggplot(ihsg_12m, aes(date, close)) +
  geom_ribbon(aes(ymin = min(close), ymax = close),
              fill = pal[["navy"]], alpha = 0.12) +
  geom_line(colour = pal[["navy"]], linewidth = 0.6) +
  geom_ma(ma_fun = SMA, n = 20, colour = pal[["orange"]],
          linetype = "solid", linewidth = 0.7) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b\n%Y") +
  scale_y_continuous(labels = comma, expand = expansion(mult = c(0, 0.05))) +
  labs(title = "JKSE Closing Price, Last 12 Months",
       subtitle = "Orange line: 20-day simple moving average",
       x = NULL, y = "Index level", caption = source_label) +
  theme_ihsg()

hd_plot(p, "02_close_12m")

Forecasting

Estimation window

train <- ihsg %>% filter(date >= forecast_start %m-% years(train_years))

nrow(train)

[1] 711

range(train$date)

[1] "2023-11-02" "2026-10-30"

Weekday effects in returns

Calendar effects are the only information about the forecast period that is known in advance. If mean returns differ by weekday, weekday dummies can make point forecasts vary from day to day; if they do not, any such variation is sampling noise. Note that Monday returns span the weekend.

ret_dow <- train %>%
  mutate(r   = 100 * log(close / lag(close)),
         dow = factor(lubridate::wday(date, week_start = 1), levels = 1:5,
                      labels = c("Mon", "Tue", "Wed", "Thu", "Fri"))) %>%
  filter(!is.na(r))

# Daily log returns (%) by weekday
ret_dow %>%
  group_by(dow) %>%
  summarise(n = n(), mean_return = mean(r), median_return = median(r), sd = sd(r)) %>%
  mutate(across(where(is.double), ~ round(.x, 3)))

# Rank-based test for differences between weekdays (robust to fat tails)
kruskal.test(r ~ dow, data = ret_dow)

    Kruskal-Wallis rank sum test

data:  r by dow
Kruskal-Wallis chi-squared = 3.0788, df = 4, p-value = 0.5447

Candidate models

Daily index levels behave close to a random walk, so the candidate set combines benchmarks for the level with a volatility model for the uncertainty:

• Random walk and random walk with drift: the benchmarks every model must beat.
• ARIMA and ETS on the log scale (lambda = 0), which keeps prediction intervals positive.
• Theta: a strong benchmark in the M3 and M4 forecasting competitions.
• TBATS with weekly (5) and annual trading-day cycles, the annual period estimated from the data.
• NNAR: a neural network autoregression.
• GJR-GARCH(1,1) with Student-t errors on daily log returns r_t = \log(P_t / P_{t-1}):

r_t = \mu + \varepsilon_t, \qquad \varepsilon_t = \sigma_t z_t, \qquad z_t \sim t_\nu

\sigma_t^2 = \omega + \left(\alpha + \gamma\, \mathbb{1}\{\varepsilon_{t-1} < 0\}\right)\varepsilon_{t-1}^2 + \beta\,\sigma_{t-1}^2

where \gamma captures the leverage effect (negative shocks raise volatility more than positive shocks of the same size). Price forecasts and intervals are obtained from 5,000 simulated return paths.

The ARIMA + weekday model regresses daily log returns (in %) on weekday dummies, with Monday as the baseline and ARMA errors:

r_t = \beta_0 + \sum_{d = \text{Tue}}^{\text{Fri}} \beta_d D_{d,t} + u_t, \qquad u_t \sim \text{ARMA}(p, q)

Its price forecasts and intervals come from 2,000 bootstrapped simulated paths.

library(rugarch)   

garch_spec <- ugarchspec(
  variance.model     = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model         = list(armaOrder = c(0, 0), include.mean = TRUE),
  distribution.model = "std"
)

# GJR-GARCH price forecasts from simulated log-return paths
fc_garch <- function(y, h, nsim = 5000) {
  fit <- ugarchfit(garch_spec, data = diff(log(y)), solver = "hybrid")
  if (convergence(fit) != 0) stop("GARCH estimation did not converge")
  sim   <- ugarchsim(fit, n.sim = h, m.sim = nsim, startMethod = "sample")
  ret   <- matrix(fitted(sim), nrow = h)
  paths <- matrix(tail(y, 1) * exp(apply(ret, 2, cumsum)), nrow = h)
  q     <- function(p) apply(paths, 1, quantile, probs = p)
  list(method = "GJR-GARCH(1,1), Student-t",
       mean   = rowMeans(paths),
       lower  = cbind(q(0.100), q(0.025)),
       upper  = cbind(q(0.900), q(0.975)),
       model  = fit)
}

# Weekday dummies (Monday is the baseline)
dow_matrix <- function(dates) {
  dow <- factor(lubridate::wday(dates, week_start = 1), levels = 1:5,
                labels = c("Mon", "Tue", "Wed", "Thu", "Fri"))
  model.matrix(~ dow)[, -1, drop = FALSE]
}

# Weekday regression with ARMA errors on log returns (%); needs the dates of the
# estimation window and of the trading days to be forecast
fc_weekday <- function(y, dates, new_dates, nsim = 2000) {
  h     <- length(new_dates)
  r     <- 100 * diff(log(y))
  fit   <- auto.arima(r, xreg = dow_matrix(dates[-1]), d = 0)
  sims  <- replicate(nsim, as.numeric(simulate(fit, nsim = h, xreg = dow_matrix(new_dates),
                                               future = TRUE, bootstrap = TRUE)))
  paths <- matrix(tail(y, 1) * exp(apply(matrix(sims, nrow = h) / 100, 2, cumsum)), nrow = h)
  q     <- function(p) apply(paths, 1, quantile, probs = p)
  list(method = sprintf("Weekday + ARMA(%d,%d) errors", fit$arma[1], fit$arma[2]),
       mean   = rowMeans(paths),
       lower  = cbind(q(0.100), q(0.025)),
       upper  = cbind(q(0.900), q(0.975)),
       model  = fit)
}

# Each function takes a numeric price vector y and forecasts h trading days ahead
fc_funs <- list(
  "Random walk"            = function(y, h) naive(ts(y), h = h, level = lv),
  "Random walk with drift" = function(y, h) rwf(ts(y), drift = TRUE, h = h, level = lv),
  "ARIMA"                  = function(y, h) forecast(auto.arima(ts(y), lambda = 0, biasadj = TRUE),
                                                     h = h, level = lv),
  "ETS"                    = function(y, h) forecast(ets(ts(y), lambda = 0, biasadj = TRUE),
                                                     h = h, level = lv),
  "Theta"                  = function(y, h) thetaf(ts(y), h = h, level = lv),
  "TBATS"                  = function(y, h) forecast(tbats(msts(y, seasonal.periods = c(5, annual_period)),
                                                           use.box.cox = TRUE, use.parallel = FALSE),
                                                     h = h, level = lv),
  "NNAR"                   = function(y, h) forecast(nnetar(ts(y), lambda = 0),
                                                     h = h, PI = TRUE, level = lv),
  "GJR-GARCH"              = function(y, h) fc_garch(y, h)
)

# Convert a forecast object (or a simulation-based list) into a tidy tibble
tidy_fc <- function(fc, model, h) {
  if (is.null(fc)) {
    return(tibble(model = model, method = NA_character_, horizon = seq_len(h),
                  mean = NA_real_, lo80 = NA_real_, hi80 = NA_real_,
                  lo95 = NA_real_, hi95 = NA_real_))
  }
  tibble(model   = model,
         method  = if (is.null(fc$method)) model else fc$method,
         horizon = seq_len(h),
         mean    = as.numeric(fc$mean),
         lo80    = as.numeric(fc$lower[, 1]), hi80 = as.numeric(fc$upper[, 1]),
         lo95    = as.numeric(fc$lower[, 2]), hi95 = as.numeric(fc$upper[, 2]))
}

# Fit the selected price-only models on price vector y; failed fits return NA rows
run_models <- function(y, h, models = names(fc_funs)) {
  fcs <- map(fc_funs[models], ~ tryCatch(.x(y, h), error = function(e) NULL))
  imap_dfr(fcs, ~ tidy_fc(.x, .y, h))
}

Forecasts for 2 to 10 November 2026

set.seed(123)
fc_objects <- map(fc_funs, ~ tryCatch(.x(train$close, h), error = function(e) NULL))

# The weekday model also needs dates; list() keeps a failed fit as NULL
fc_objects["ARIMA + weekday"] <- list(
  tryCatch(fc_weekday(train$close, train$date, future_dates), error = function(e) NULL)
)

fc_final <- imap_dfr(fc_objects, ~ tidy_fc(.x, .y, h)) %>%
  mutate(date = future_dates[horizon])

# Specification selected for each model
fc_final %>% distinct(model, method)

garch_fit <- fc_objects[["GJR-GARCH"]]$model
wd_fit    <- fc_objects[["ARIMA + weekday"]]$model

# GJR-GARCH estimates, persistence, and half-life of volatility shocks (trading days)
if (!is.null(garch_fit)) {
  print(round(garch_fit@fit$matcoef, 4))
  print(c(persistence = persistence(garch_fit), half_life_days = halflife(garch_fit)))
}

        Estimate  Std. Error  t value Pr(>|t|)
mu        0.0005      0.0003   1.6035   0.1088
omega     0.0000      0.0000  38.2196   0.0000
alpha1    0.0000      0.0140   0.0000   1.0000
beta1     0.7942      0.0267  29.7512   0.0000
gamma1    0.2553      0.0585   4.3632   0.0000
shape     5.3398      0.8927   5.9814   0.0000
   persistence half_life_days 
     0.9218104      8.5136758 

# Weekday regression with ARMA errors (returns in %)
if (!is.null(wd_fit)) print(wd_fit)

Series: r 
Regression with ARIMA(3,0,0) errors 

Coefficients:
         ar1      ar2     ar3  dowTue  dowWed   dowThu   dowFri
      0.0669  -0.0799  0.0829  0.0722  0.0650  -0.0363  -0.2137
s.e.  0.0374   0.0377  0.0376  0.1127  0.1132   0.1132   0.1146

sigma^2 = 1.847:  log likelihood = -1221.75
AIC=2459.51   AICc=2459.71   BIC=2496.03

# Point forecasts by model
fc_final %>%
  select(date, model, mean) %>%
  pivot_wider(names_from = model, values_from = mean) %>%
  mutate(across(-date, ~ round(.x, 1)))

# Point forecasts with 80% and 95% prediction intervals
fc_final %>%
  select(model, date, mean, lo80, hi80, lo95, hi95) %>%
  mutate(across(mean:hi95, ~ round(.x, 1)))

plot_fc <- function(fc, hist, ncol = 2) {
  fc       <- fc %>% mutate(model = factor(model, levels = unique(model)))
  last_obs <- slice_tail(hist, n = 1)

  # Connect the last observation to each forecast line
  bridge <- distinct(fc, model) %>%
    mutate(date = last_obs$date, mean = last_obs$close,
           lo80 = mean, hi80 = mean, lo95 = mean, hi95 = mean)

  fc_line  <- bind_rows(fc, bridge) %>% arrange(model, date)
  hist_rep <- cross_join(distinct(fc, model), select(hist, date, close))

  ggplot() +
    geom_ribbon(data = fc_line, aes(date, ymin = lo95, ymax = hi95, fill = "95% interval"),
                alpha = 0.35) +
    geom_ribbon(data = fc_line, aes(date, ymin = lo80, ymax = hi80, fill = "80% interval"),
                alpha = 0.55) +
    geom_line(data = hist_rep, aes(date, close, colour = "Actual"), linewidth = 0.5) +
    geom_line(data = fc_line, aes(date, mean, colour = "Point forecast"), linewidth = 0.7) +
    geom_point(data = fc, aes(date, mean, colour = "Point forecast"), size = 1.2) +
    facet_wrap(~ model, ncol = ncol) +
    scale_fill_manual(values = c("95% interval" = "#c6d9ec", "80% interval" = "#8fb3d9"),
                      breaks = c("80% interval", "95% interval")) +
    scale_colour_manual(values = c("Actual" = "grey25", "Point forecast" = pal[["orange"]])) +
    scale_x_date(date_breaks = "2 weeks", date_labels = "%d %b") +
    scale_y_continuous(labels = comma) +
    labs(x = NULL, y = "Index level") +
    theme_ihsg()
}

# Use the selected specification as the panel title where it is informative
fc_plot_data <- fc_final %>%
  mutate(model = if_else(model %in% c("ARIMA", "ETS", "TBATS", "NNAR",
                                      "GJR-GARCH", "ARIMA + weekday"),
                         coalesce(method, model), model))

p <- plot_fc(fc_plot_data, slice_tail(train, n = 40), ncol = 3) +
  labs(title = "JKSE Closing Price Forecasts, 2 to 10 November 2026",
       subtitle = "Last 40 trading days with point forecasts and 80% / 95% prediction intervals",
       caption = source_label)

hd_plot(p, "03_forecasts_by_model", width = 10, height = 9)

Simulated price paths

The figure shows 30 of 2,000 paths simulated from the fitted GJR-GARCH model. Individual paths move up and down, but because the direction of each move is unpredictable, their average is nearly flat. The bands reflect the current level of volatility.

if (is.null(garch_fit)) stop("The GJR-GARCH fit is not available")

set.seed(123)
n_paths    <- 2000
last_close <- tail(train$close, 1)

sim   <- ugarchsim(garch_fit, n.sim = h, m.sim = n_paths, startMethod = "sample")
ret   <- matrix(fitted(sim), nrow = h)
paths <- rbind(rep(last_close, n_paths),
               last_close * exp(apply(ret, 2, cumsum)))
colnames(paths) <- paste0("path_", seq_len(n_paths))

path_tbl <- as_tibble(paths) %>%
  mutate(date = c(max(train$date), future_dates)) %>%
  pivot_longer(-date, names_to = "path", values_to = "close")

band_tbl <- path_tbl %>%
  group_by(date) %>%
  summarise(mean = mean(close),
            lo80 = quantile(close, 0.100), hi80 = quantile(close, 0.900),
            lo95 = quantile(close, 0.025), hi95 = quantile(close, 0.975),
            .groups = "drop")

shown <- sample(colnames(paths), 30)

p <- ggplot() +
  geom_ribbon(data = band_tbl, aes(date, ymin = lo95, ymax = hi95),
              fill = "#c6d9ec", alpha = 0.7) +
  geom_ribbon(data = band_tbl, aes(date, ymin = lo80, ymax = hi80),
              fill = "#8fb3d9", alpha = 0.7) +
  geom_line(data = filter(path_tbl, path %in% shown),
            aes(date, close, group = path),
            colour = pal[["navy"]], alpha = 0.4, linewidth = 0.35) +
  geom_line(data = slice_tail(train, n = 40), aes(date, close),
            colour = "grey25", linewidth = 0.6) +
  geom_line(data = band_tbl, aes(date, mean), colour = pal[["orange"]], linewidth = 1) +
  scale_x_date(date_breaks = "2 weeks", date_labels = "%d %b") +
  scale_y_continuous(labels = comma) +
  labs(title = "Simulated JKSE Price Paths, 2 to 10 November 2026",
       subtitle = "30 of 2,000 GJR-GARCH paths (blue), their mean (orange), and 80% / 95% bands",
       x = NULL, y = "Index level", caption = source_label) +
  theme_ihsg()

hd_plot(p, "04_simulated_paths")

Rolling-origin evaluation

A single holdout of a few days is too noisy to rank models. Each model is re-estimated on a rolling window of the same length as the final estimation window, at forecast origins spaced 5 trading days apart over the last 250 trading days. The ARIMA + weekday model uses the dates of each window and the actual trading dates that follow each origin, which are known in advance from the IDX calendar. Errors are averaged over all origins and horizons. Coverage close to 0.80 and 0.95 indicates well-calibrated prediction intervals. TBATS and NNAR are only included when run_slow_backtest = TRUE.

bt_models <- c("Random walk", "Random walk with drift", "ARIMA", "ETS", "Theta",
               "GJR-GARCH", "ARIMA + weekday")
if (run_slow_backtest) bt_models <- c(bt_models, "TBATS", "NNAR")

y_all   <- ihsg$close
n_win   <- nrow(train)                               # same window length as the final models
origins <- seq(length(y_all) - 250, length(y_all) - h, by = 5)

set.seed(123)
bt <- map_dfr(origins, function(o) {
  idx <- (o - n_win + 1):o

  # Models that use prices only
  fc_price <- run_models(y_all[idx], h, models = intersect(bt_models, names(fc_funs)))

  # Weekday model: also needs the window dates and the next h trading dates
  wd_obj <- tryCatch(fc_weekday(y_all[idx], ihsg$date[idx], ihsg$date[o + seq_len(h)],
                                nsim = 1000),
                     error = function(e) NULL)

  bind_rows(fc_price, tidy_fc(wd_obj, "ARIMA + weekday", h)) %>%
    mutate(origin      = o,
           origin_date = ihsg$date[o],
           actual      = y_all[o + horizon])
})

saveRDS(bt, "backtest_results.rds")   

bt %>%
  group_by(model) %>%
  summarise(
    RMSE           = sqrt(mean((actual - mean)^2, na.rm = TRUE)),
    MAE            = mean(abs(actual - mean), na.rm = TRUE),
    MAPE           = 100 * mean(abs((actual - mean) / actual), na.rm = TRUE),
    `Coverage 80%` = mean(actual >= lo80 & actual <= hi80, na.rm = TRUE),
    `Coverage 95%` = mean(actual >= lo95 & actual <= hi95, na.rm = TRUE),
    `Failed fits`  = sum(is.na(mean)) / h,
    .groups = "drop"
  ) %>%
  mutate(`RMSE vs RW` = RMSE / RMSE[model == "Random walk"]) %>%
  relocate(`RMSE vs RW`, .after = RMSE) %>%
  arrange(RMSE) %>%
  mutate(across(where(is.numeric), ~ round(.x, 3)))

rmse_h <- bt %>%
  group_by(model, horizon) %>%
  summarise(RMSE = sqrt(mean((actual - mean)^2, na.rm = TRUE)), .groups = "drop") %>%
  group_by(horizon) %>%
  mutate(rel_rmse = RMSE / RMSE[model == "Random walk"]) %>%
  ungroup() %>%
  mutate(model = factor(model, levels = bt_models))

# Nine distinct colours, enough for all models when run_slow_backtest = TRUE
model_cols <- c("#1f4e79", "#e07b39", "#2e8b57", "#c0392b", "#7b3294",
                "#8c6d31", "#17becf", "#e377c2", "grey40")

p <- ggplot(rmse_h, aes(horizon, rel_rmse, colour = model)) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey40") +
  geom_line(linewidth = 0.8) +
  geom_point(size = 1.8) +
  scale_colour_manual(values = model_cols) +
  scale_x_continuous(breaks = seq_len(h)) +
  guides(colour = guide_legend(nrow = 2)) +
  labs(title = "Forecast Accuracy Relative to the Random Walk",
       subtitle = "RMSE ratio by horizon across all origins; values below 1 beat the random walk",
       x = "Forecast horizon (trading days)", y = "Relative RMSE",
       caption = source_label) +
  theme_ihsg()

hd_plot(p, "05_relative_rmse", height = 4.5)

Diebold–Mariano tests compare one-step-ahead squared errors of each model with the random walk. A negative statistic means the model has smaller errors; a p-value below 0.05 indicates that the difference is statistically significant.

e1 <- bt %>%
  filter(horizon == 1) %>%
  transmute(origin, model, error = actual - mean) %>%
  pivot_wider(names_from = model, values_from = error)

map_dfr(setdiff(bt_models, "Random walk"), function(m) {
  ok   <- complete.cases(e1[[m]], e1[["Random walk"]])
  test <- dm.test(e1[[m]][ok], e1[["Random walk"]][ok], h = 1, power = 2)
  tibble(model = m, DM_statistic = unname(test$statistic), p_value = test$p.value)
}) %>%
  mutate(across(where(is.numeric), ~ round(.x, 3)))