When forecasting high-volatility financial assets like Bitcoin (BTC-USD), standard point-prediction models—whether classical ARIMA or deep neural networks—often fall short. A point estimate tells you where the model expects the price to go, but it completely ignores market regime shifts and changing risk levels.

In this article, we construct an end-to-end forecasting pipeline in R that combines Machine Learning (Multi-Layer Perceptron via brulee and tidymodels) for directional trend prediction with Financial Econometrics (GJR-GARCH via rugarch) to generate dynamic, volatility-adjusted confidence bounds.

The Architecture Overview

Instead of attempting to forecast level prices directly—which introduces non-stationarity issues—our pipeline operates on log returns and reconstructs price trajectories dynamically.

Key Architectural Decisions & Methodological Rationale

1. Modeling Log Returns Instead of Level Prices

Directly predicting raw prices often creates spurious regression results because price series are non-stationary. By transforming the price series to log returns:

This transformation achieves stationarity, making the target variable suitable for neural network training. Factual price predictions are reconstructed at test time via cumulative exponentiation:

2. Deep Learning for Trend: Multi-Layer Perceptron (brulee)

We use a 2-layer MLP architecture (32, 16) initialized via the tidymodels framework.

• Lag Features: We supply 20 days of lagged log returns to capture short- to medium-term autocorrelation.
• Normalizations & Dropping Zero Variance: Handled within a clean recipes pipeline to prevent data leakage

3. Econometric Volatility Modeling: GJR-GARCH(1,1) on Neural Network Residuals

Machine learning models capture the conditional mean, but their residuals, $\epsilon_t = r_t – \hat{r}_t$ , still contain volatility clustering—periods of high turbulence followed by high turbulence, and calm periods followed by calm.

To model this conditional variance ($\sigma_t^2$), we fit a GJR-GARCH(1,1) model on the MLP training residuals:

• Why GJR-GARCH over Standard GARCH? Crypto markets exhibit leverage effects: price drops often spark higher panic volatility than equivalent price rallies. The indicator term $I_{t-1}$ captures this asymmetry ($\gamma$).
• Why garchOrder = c(1, 1)? The order (1,1) specifies a single lag for the ARCH component ($\epsilon_{t-1}^2$) and a single lag for the GARCH component ($\sigma_{t-1}^2$). Empirical literature demonstrates that a single lag captures the vast majority of volatility persistence without overfitting.
• Why Skewed Student-t (sstd) Distribution? Bitcoin returns violate Gaussian normal distribution assumptions. They display fat tails (extreme events occur frequently) and skewness (asymmetric distribution). The sstd distribution explicitly models both heavy tails and skewness.
• Solver Selection: We use solver = "hybrid" in rugarch to combine optimization algorithms, ensuring convergence even under complex distribution parameters.

When the market is quiet, these confidence bounds contract automatically. During volatile regimes, the bounds expand, giving a realistic probabilistic envelope rather than static, uninformative confidence intervals.

Complete R Implementation

Below is the complete, self-contained R script implementing the full pipeline using tidymodels, brulee, rugarch, and plotly.

# Integrated Forecasting Pipeline with tidymodels, yardstick, and plotly
library(tidymodels)
library(timetk)
library(brulee)
library(modeltime)
library(tidyquant)
library(plotly)
library(scales)
library(glue)
library(ggtext)
library(rugarch) # Added for GJR-GARCH volatility modeling

# 1. Fetch Data & Calculate Log Returns
df_raw <- tq_get("BTC-USD") %>%
  select(date, close) %>%
  filter(date >= last(date) - months(24)) %>%
  drop_na() %>%
  mutate(log_ret = log(close / lag(close))) %>%
  drop_na()

# 2. Add Lag Features on Log Returns
df_prepared <- df_raw %>%
  tk_augment_lags(log_ret, .lags = 1:20) %>%
  drop_na()

# 3. Time Series Split (10 Assessment Window)
split <- time_series_split(
  df_prepared,
  assess = 10,
  cumulative = TRUE
)

df_train <- training(split)
df_test  <- testing(split)

# 4. Recipe Setup (Predict Log Return, Not Level Price)
recipe_mlp <- recipe(log_ret ~ ., data = df_train) %>%
  step_rm(date, close) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_zv(all_predictors()) %>%
  step_normalize(all_numeric_predictors())

# modeltime::naive_reg() requires the date column as a predictor to parse the time index
recipe_naive <- recipe(log_ret ~ date, data = df_train)

# 5. Model Specifications
model_spec_mlp <- mlp(
  hidden_units = c(32, 16),
  epochs = 100,
  dropout = 0.05,
  learn_rate = 0.005
) %>%
  set_engine("brulee") %>%
  set_mode("regression")

model_spec_naive <- naive_reg() %>%
  set_engine("naive")

# 6. Workflow Construction & Model Fitting
workflow_mlp <- workflow() %>%
  add_recipe(recipe_mlp) %>%
  add_model(model_spec_mlp) %>%
  fit(df_train)

workflow_naive <- workflow() %>%
  add_recipe(recipe_naive) %>%
  add_model(model_spec_naive) %>%
  fit(df_train)

# 7. Predictions & Price Reconstruction (Reference Code Logic)
preds_mlp_ret   <- predict(workflow_mlp, df_test)
preds_naive_ret <- predict(workflow_naive, df_test)

df_results <- df_test %>%
  mutate(
    yesterday_close = lag(close, default = first(close)),
    actual      = close,
    pred_brulee = yesterday_close * exp(preds_mlp_ret$.pred),
    pred_naive  = yesterday_close * exp(preds_naive_ret$.pred)
  )

# ==============================================================================
# 7.5. GJR-GARCH (Skewed Student-t) Dynamic Volatility & Confidence Bands
# ==============================================================================
# Calculate MLP residuals on the training set
train_preds_mlp <- predict(workflow_mlp, df_train)$.pred
train_residuals <- df_train$log_ret - train_preds_mlp

# GJR-GARCH(1,1) specification with Skewed Student-t distribution for fat-tails and asymmetry
garch_spec <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 0), include.mean = FALSE),
  distribution.model = "sstd"
)

# Fit GJR-GARCH model on training residuals
garch_fit <- ugarchfit(spec = garch_spec, data = train_residuals, solver = "hybrid")

# Forecast conditional volatility for the 10-day test horizon
garch_fore <- ugarchforecast(garch_fit, n.ahead = 10)
dynamic_sigmas <- sigma(garch_fore) # 10-day dynamic conditional standard deviations

# Attach dynamic standard deviations to test results and compute dynamic confidence bounds
df_results <- df_results %>%
  mutate(
    conditional_sigma = as.vector(dynamic_sigmas),
    # Dynamically scaled confidence bounds in price space (± 1.96 * conditional_sigma)
    conf_hi = pred_brulee * exp(1.96 * conditional_sigma),
    conf_lo = pred_brulee * exp(-1.96 * conditional_sigma)
  )

# 8. Accuracy Table Construction (Yardstick on Original Price Scale)
accuracy_table <- tibble(
  .model_id   = c(1, 2),
  .model_desc = c("BRULEE MLP (32, 16) + GJR-GARCH", "NAIVE"),
  rmse        = c(rmse_vec(df_results$actual, df_results$pred_brulee),
                  rmse_vec(df_results$actual, df_results$pred_naive)),
  rsq         = c(rsq_vec(df_results$actual, df_results$pred_brulee),
                  rsq_vec(df_results$actual, df_results$pred_naive)),
  mape        = c(mape_vec(df_results$actual, df_results$pred_brulee),
                  mape_vec(df_results$actual, df_results$pred_naive)),
  mae         = c(mae_vec(df_results$actual, df_results$pred_brulee),
                  mae_vec(df_results$actual, df_results$pred_naive))
)

# Print Model Comparison Table
cat("=== MODEL ACCURACY EVALUATION TABLE (PRICE SCALE) ===\n")
cat(knitr::kable(accuracy_table, digits = 4, format = "simple"), sep = "\n")

# 9. Modern Interactive Plotly Visualization (Reference Style)
brulee_mape <- accuracy_table$mape[1]
df_eval <- df_results

# Custom Hover Text Data Frames
df_plot_actual <- df_eval %>% select(date, actual)  %>% mutate(text_actual = glue("<b>Actual Price:</b> {scales::dollar(actual, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_pred   <- df_eval %>% select(date, pred_brulee) %>% mutate(text_pred   = glue("<b>AI Prediction Price:</b> {scales::dollar(pred_brulee, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_hi     <- df_eval %>% select(date, conf_hi) %>% mutate(text_hi     = glue("<b>GJR-GARCH Upper Bound:</b> {scales::dollar(conf_hi, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))
df_plot_lo     <- df_eval %>% select(date, conf_lo) %>% mutate(text_lo     = glue("<b>GJR-GARCH Lower Bound:</b> {scales::dollar(conf_lo, accuracy = 0.01)}<br><b>Date:</b> {format(date, '%b %d, %Y')}"))

# Build Base GGPlot Layer
p <- ggplot() +
  geom_ribbon(data = df_eval, aes(x = date, ymin = conf_lo, ymax = conf_hi), fill = "#808080", alpha = 0.18) +
  geom_point(data = df_plot_hi, aes(x = date, y = conf_hi, text = text_hi), color = "transparent", alpha = 0, size = 3) +
  geom_point(data = df_plot_lo, aes(x = date, y = conf_lo, text = text_lo), color = "transparent", alpha = 0, size = 3) +
  geom_line(data = df_plot_actual, aes(x = date, y = actual), color = "#2c3e50", linewidth = 1.2) +
  geom_point(data = df_plot_actual, aes(x = date, y = actual, text = text_actual), color = "#2c3e50", size = 2) +
  geom_line(data = df_plot_pred, aes(x = date, y = pred_brulee), color = "#e74c3c", linetype = "dashed", linewidth = 1.2) +
  geom_point(data = df_plot_pred, aes(x = date, y = pred_brulee, text = text_pred), color = "#e74c3c", size = 2) +
  scale_y_continuous(labels = dollar_format(accuracy = 1)) +
  labs(
    x = "", y = "",
    title = paste0(
      "<span style='color:#2c3e50'>BTC-USD Prices</span> vs ",
      "<span style='color:#e74c3c'>MLP + GJR-GARCH Dynamic Forecast</span><br>",
      "<span style='font-size:12px; color:#555555;'>10-Day Horizon | Out-of-Sample MAPE: ", round(brulee_mape, 2), "</span>"
    )
  ) +
  theme_minimal() +
  theme(
    plot.title = element_markdown(hjust = 0.5, face = "bold"),
    plot.background = element_rect(fill = "#ffffff", color = NA),
    panel.background = element_rect(fill = "#ffffff", color = NA),
    panel.grid.minor = element_blank()
  )

# Typography Settings & Plotly Output
font_family <- list(family = "Roboto Slab, Sans-Serif", size = 16)
label_font  <- list(font = list(family = "Roboto Slab, Sans-Serif", size = 13))

ggplotly(p, tooltip = "text") %>%
  style(hoverlabel = label_font) %>%
  layout(font = font_family)

Conclusion

Combining machine learning with traditional econometrics produces a hybrid system that overcomes the limitations of both fields:

1. Machine Learning captures non-linear relationships and trend movements.
2. Econometrics (GARCH) supplies explicit statistical parameters for uncertainty and volatility clustering.

The resulting output provides realistic, dynamic risk envelopes that evolve with changing market regimes—giving quantitative traders and analysts a robust framework for volatility-aware decision-making.