[This article was first published on DataGeeek, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.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 OverviewInstead 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 Rationale1. Modeling Log Returns Instead of Level PricesDirectly predicting raw prices often creates spurious regression results because price series are non-stationary. By transforming the price series to log returns:rt=ln(PtPt−1)r_t = \ln\left(\frac{P_t}{P_{t-1}}\right)This transformation achieves stationarity, making the target variable suitable for neural network training. Factual price predictions are reconstructed at test time via cumulative exponentiation:P^t=Pt−1×exp(r^t)\hat{P}_t = P_{t-1} \times \exp(\hat{r}_t)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 leakage3. Econometric Volatility Modeling: GJR-GARCH(1,1) on Neural Network ResidualsMachine learning models capture the conditional mean, but their residuals, ϵt=rt−r^t \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 (σt2\sigma_t^2), we fit a GJR-GARCH(1,1) model on the MLP training residuals:σt2=ω+αϵt−12+γIt−1ϵt−12+βσt−12\sigma_t^2 = \omega + \alpha \epsilon_{t-1}^2 + \gamma I_{t-1} \epsilon_{t-1}^2 + \beta \sigma_{t-1}^2Why GJR-GARCH over Standard GARCH? Crypto markets exhibit leverage effects: price drops often spark higher panic volatility than equivalent price rallies. The indicator term It−1I_{t-1} captures this asymmetry (γ\gamma).Why garchOrder = c(1, 1)? The order (1,1) specifies a single lag for the ARCH component (ϵt−12\epsilon_{t-1}^2) and a single lag for the GARCH component (σt−12\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.Upper Boundt=P^t×exp(+1.96×σt)\text{Upper Bound}_t = \hat{P}_t \times \exp(+1.96 \times \sigma_t)Lower Boundt=P^t×exp(−1.96×σt)\text{Lower Bound}_t = \hat{P}_t \times \exp(-1.96 \times \sigma_t)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 ImplementationBelow is the complete, self-contained R script implementing the full pipeline using tidymodels, brulee, rugarch, and plotly.# Integrated Forecasting Pipeline with tidymodels, yardstick, and plotlylibrary(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 Returnsdf_raw % 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 Returnsdf_prepared % tk_augment_lags(log_ret, .lags = 1:20) %>% drop_na()# 3. Time Series Split (10 Assessment Window)split % step_normalize(all_numeric_predictors())# modeltime::naive_reg() requires the date column as a predictor to parse the time indexrecipe_naive % set_mode("regression")model_spec_naive % set_engine("naive")# 6. Workflow Construction & Model Fittingworkflow_mlp % add_recipe(recipe_mlp) %>% add_model(model_spec_mlp) %>% fit(df_train)workflow_naive % add_recipe(recipe_naive) %>% add_model(model_spec_naive) %>% fit(df_train)# 7. Predictions & Price Reconstruction (Reference Code Logic)preds_mlp_ret % mutate(text_hi = glue("GJR-GARCH Upper Bound: {scales::dollar(conf_hi, accuracy = 0.01)}Date: {format(date, '%b %d, %Y')}"))df_plot_lo % select(date, conf_lo) %>% mutate(text_lo = glue("GJR-GARCH Lower Bound: {scales::dollar(conf_lo, accuracy = 0.01)}Date: {format(date, '%b %d, %Y')}"))# Build Base GGPlot Layerp