Time Series Forecasting Pipeline

A comprehensive four-part machine learning pipeline for financial time series forecasting, combining fractional Brownian bridges, XGBoost regression, and Monte Carlo simulation for probabilistic price predictions.

Overview

This project forecasts future stock prices by:

1. Part 1: Estimating volatility and market mean-reversion characteristics
2. Part 2: Generating synthetic price paths using fractional Brownian motion
3. Part 3: Training an XGBoost model on synthetic + real data
4. Part 4: Running Monte Carlo simulations and fitting kernel density estimators

The output is a probabilistic price forecast with confidence intervals, not just point estimates.

---

Quick Start

Installation

# Clone or download this repository
cd time_series_forecasting_pipeline

# Install dependencies
pip install -r requirements.txt

Basic Usage

# Run the entire pipeline with default settings (AAPL, June 2024 - Dec 2024)
python run_all.py

# Run with visualizations
python run_all.py --visualize

# Run with custom parameters
python run_all.py --ticker MSFT --start 2024-01-01 --end 2024-12-31 --horizon 5 --nsim 15000 --visualize

Command-Line Options

--ticker TICKER           Stock ticker symbol (default: AAPL)
--start DATE             Start date in YYYY-MM-DD format (default: 2024-06-01)
--end DATE               End date in YYYY-MM-DD format (default: 2024-12-31)
--horizon DAYS           Forecast horizon in days (default: 3)
--nsim N                 Number of Monte Carlo paths (default: 10000)
--bandwidth METHOD       KDE bandwidth: 'scott', 'silverman', or float (default: scott)
--visualize              Show plots for all stages (default: off)
--quiet                  Suppress verbose output (default: off)

Example Runs

# Minimal output, no plots
python run_all.py --quiet

# Full diagnostics with visualizations, Tesla 1-year forecast
python run_all.py --ticker TSLA --start 2023-12-31 --end 2024-12-31 --visualize

# 5-day forecast with 20k Monte Carlo paths
python run_all.py --horizon 5 --nsim 20000 --visualize

---

Pipeline Architecture

Part 1: Data Preparation & Volatility Estimation

Goal: Clean daily OHLCV data and estimate market volatility and mean-reversion characteristics.

Steps:

1. Data Loading & Cleaning

   • Download historical data from Yahoo Finance
   • Validate OHLC relationships (High ≥ Open, Close, Low; Low ≤ all)
   • Remove days with zero volume or missing data

2. Log Returns

   • Compute daily log-returns: $r_t = \ln(P_t / P_{t-1})$
   • Report mean, std, skewness, kurtosis for sanity checks

3. Yang–Zhang Volatility (parkinson-like, scale-robust)

   • Incorporates High, Low, Open, Close
   • More robust to overnight gaps than simple historical volatility
   • Applied over rolling window (default 20 days)

4. EWMA Volatility (exponential moving average)

   • Reweights recent observations more heavily
   • Decay factor λ = 0.94 (default)
   • Captures time-varying volatility

5. Blended Volatility

   • Equal weight average of Yang–Zhang and EWMA
   • Extended to full price length via spline interpolation
   • Balanced view of volatility across windows

6. Hurst Exponent (R/S analysis)

   • Measures market mean-reversion (H < 0.5) vs momentum (H > 0.5)
   • H ≈ 0.5: random walk; H ≈ 0.36: mean-reverting (AAPL in 2024)
   • Used to parameterize synthetic paths in Part 2

Outputs:

• prices: daily closing prices (n,)
• sigma_extended: daily volatility estimates (n,)
• hurst: scalar mean-reversion measure
• Plots: price vs time, volatility vs time, returns histogram (if --visualize)

Key Insight: Volatility is time-varying and often correlates with price jumps; incorporating both YZ and EWMA smooths noisy single-window estimates.

---

Part 2: Synthetic Path Generation & Feature Extraction

Goal: Generate many synthetic price paths using fractional Brownian motion, then extract lagged features and targets for supervised learning.

Key Concepts:

Fractional Brownian Motion (fBM)

• Standard Brownian motion (H=0.5) captures pure randomness
• fBM with H < 0.5 adds mean-reversion: prices drift back to trend
• Covariance structure: $\mathrm{Cov}(B_s^H, B_t^H) = \frac{1}{2}(s^{2H} + t^{2H} - |s-t|^{2H})$

Brownian Bridge Interpolation

• Generate fBM within each daily interval [t, t+1]
• Enforce boundary conditions: path starts at P_t, ends at P_{t+1}
• Apply local volatility scaling: multiply increments by σ(t) to match daily forecasts

Steps:

1. Generate n_paths fBB paths (default 100)

   • For each synthetic path:
     • For each daily interval [t, t+1]:
       • Draw fBM increments with covariance matrix from H estimate
       • Cholesky factorize to generate correlated noise
       • Compute fBB bridge (ensure boundary continuity)
       • Scale increments by local volatility σ(t)
     • Connect intervals to get full price path
   • Densify each day into n_substeps time points (default 10)
   • Result: each path has (n_days × n_substeps + 1) points

2. Include Real Data

   • Linear interpolation of actual prices as one extra path
   • Useful benchmark for visualization

3. Extract Features

   • For each path, every position i:
     • Features: last lag_y prices + last lag_sigma volatilities
     • Target: price at position i+1
   • Aggregate across all paths and positions

4. Diagnostics

   • Print X_train mean/std of first 5 features
   • Print y_train mean/std (should match price level)
   • Plots: sample fBB paths vs real prices, y_train histogram, last-price vs next-price scatter

Outputs:

• X_train: shape (n_total_samples, lag_y + lag_sigma), e.g., (145945, 10)
• y_train: shape (n_total_samples,), next price targets
• sample_types: 0 for synthetic, 1 for real
• Plots: overlaid synthetic/real paths, target distribution (if --visualize)

Key Insight: By generating many synthetic paths with realistic volatility scaling, we create a large, diverse training set without data snooping. The fBM Hurst parameter ensures synthetic paths share mean-reversion characteristics with the real market.

---

Part 3: XGBoost Model Training

Goal: Fit a gradient boosting regression model to predict next-day price from lagged price/vol features.

Steps:

1. Train/Test Split

   • Divide X_train / y_train: 80% train, 20% test (random)
   • Stratified by synthetic vs real (currently not stratified but could be)

2. Feature Standardization

   • Fit StandardScaler on training features only
   • Apply same transformation to test features
   • Prevents data leakage

3. XGBoost Configuration

   • n_estimators: 100 trees
   • max_depth: 6 (moderate complexity)
   • learning_rate: 0.1 (step size per boosting round)
   • subsample: 0.8 (row subsampling)
   • colsample_bytree: 0.8 (column subsampling)
   • Random state: 42 (reproducibility)

4. Training & Evaluation

   • Fit on X_train_scaled, y_train_split
   • Predict on both training and test sets
   • Compute metrics:
     • MSE (mean squared error)
     • MAE (mean absolute error)
     • R² (coefficient of determination)

5. Residual Analysis

   • Compute residuals on test set: $\varepsilon = y_{\mathrm{test}} - \hat{y}_{\mathrm{test}}$
   • Report mean (should be ≈0), std, min, max
   • Plot histogram of residuals (if --visualize)

6. Feature Importance

   • Extract feature_importances_ from XGBoost
   • Print top 5 features

Outputs:

• model: trained XGBRegressor
• scaler: fitted StandardScaler
• Metrics: train/test MSE, MAE, R²
• residuals_test: residuals on held-out test set
• Plots: residual histogram, true vs predicted scatter (if --visualize)

Typical Results (AAPL 2024):

• Train R² ≈ 0.9999, Test R² ≈ 0.9999
• Train MAE ≈ 0.09, Test MAE ≈ 0.09
• Residual std ≈ 0.12 (~0.05% of price level)

Key Insight: XGBoost learns the conditional mean $E\big[P_{t+1} \mid P_t, \ldots, P_{t-\ell_y}, \sigma_t, \ldots, \sigma_{t-\ell_\sigma}\big]$, where $\ell_y$ and $\ell_\sigma$ are the respective lag counts. The tight R² on densified synthetic data means prices are highly predictable from recent history—a consequence of the smooth fBM bridges. In production, test on truly held-out future data.

---

Part 4: Monte Carlo Simulation & KDE

Goal: Forecast a full predictive distribution (not just a point estimate) by running many forward simulations and fitting kernel density estimators.

Key Concepts:

Monte Carlo Simulation

• Use the trained XGBoost as a one-step-ahead conditional mean model
• For each simulation s = 1, ..., n_sim:
  • Start from current state: last lag_y prices + lag_sigma volatilities
  • For each horizon h = 1, ..., steps_ahead:
    • Build feature vector from lagged prices/vols
    • Standardize with fitted scaler
    • Predict mean: $\hat{\mu}{t+h} = f(X{t+h})$
    • Sample noise: $\varepsilon \sim \mathrm{Residuals}_{\mathrm{test}}$ (residual bootstrap)
    • Simulated price: $P_{t+h}^{(s)} = \hat{\mu}_{t+h} + \varepsilon$
    • Update lagged prices/vols for next step

Kernel Density Estimation (KDE)

• For each horizon h:
  • Collect terminal prices: ${P_{t+h}^{(1)}, \ldots, P_{t+h}^{(n_{\mathrm{sim}})}}$
  • Fit univariate Gaussian KDE: $f_h(x) = \frac{1}{n_{\mathrm{sim}}} \sum_{i=1}^{n_{\mathrm{sim}}} K(x - P_{t+h}^{(i)})$
  • Bandwidth options:
    • Scott's rule: $h \approx 1.06 \sigma n^{-1/5}$ (default, good for smooth distributions)
    • Silverman's rule: $h \approx 0.9 \sigma n^{-1/5}$ (slightly more conservative)
    • Custom: pass a float as scaling factor

Steps:

1. Compute Residuals

   • From Part 3 test set (or training, depending on philosophy)
   • Used for noise injection during simulation

2. Simulate Paths

   • Initialize state from last lag_y prices + lag_sigma vols in real data
   • For each MC path:
     • Iteratively apply model + noise
     • Save terminal price at each horizon

3. Fit KDE

   • Per horizon: take all simulated terminal prices
   • Fit Gaussian KDE with bandwidth method
   • Print statistics: mean, std, percentiles (1, 5, 25, 50, 75, 95, 99)

4. Diagnostics & Visualization

   • Compute terminal price stats and returns relative to last real price
   • Plot histogram + KDE curve per horizon
   • Vertical lines: last real price (green), mean simulated price (red)
   • Interpret: where is the predicted price likely to be?

Outputs:

• sims: shape (n_sim, steps_ahead), all simulated prices
• kde_funcs: list of KDE objects per horizon
• residuals: residuals used for noise
• Plots: histogram + KDE per horizon (if --visualize)

Example Output (AAPL, 3-day forecast, n_sim=10000):

Last real price: 251.0645
Horizon 1:
  Mean: 252.0168, Std: 0.1234
  5th pct: 251.8248, 95th pct: 252.2130

Horizon 2:
  Mean: 251.4394, Std: 0.1267
  5th pct: 251.2373, 95th pct: 251.6355

Horizon 3:
  Mean: 251.3182, Std: 0.1541
  5th pct: 251.0873, 95th pct: 251.5602
  Mean return: 0.1010%

Key Insight: Monte Carlo simulates the trajectory distribution. KDE smooths the empirical samples into a continuous density. Together, they answer: "Given the model and historical noise, what is the probability distribution of prices in 1, 2, 3 days?" This enables risk quantification (e.g., "5% chance of drop below 251") without parametric assumptions.

---

File Structure

time_series_forecasting_pipeline/
├── README.md                   # This file
├── requirements.txt            # Dependencies
├── run_all.py                  # Main orchestrator
├── data_gen.py                 # Yahoo Finance data loader
├── part1.py                    # Data + volatility estimation
├── part2.py                    # fBB path generation + features
├── part3.py                    # XGBoost training
└── part4.py                    # Monte Carlo + KDE

---

Dependencies

numpy>=1.21.0
pandas>=1.3.0
yfinance>=0.2.0
scikit-learn>=1.0.0
xgboost>=1.5.0
scipy>=1.7.0
matplotlib>=3.4.0

Install with:

pip install -r requirements.txt

---

Interpretation & Usage

Interpreting the Output

1. Part 1 Hurst Exponent

   • H < 0.5: Mean-reverting market (prices bounce around trend)
   • H ≈ 0.5: Random walk (no predictability from history)
   • H > 0.5: Trending market (momentum)
   • Example: AAPL 2024 H ≈ 0.36 → mean-reverting

2. Part 3 R² Score

   • R² ≈ 1: Model fits training data almost perfectly
   • R² ≈ 0.5: Model explains 50% of variance
   • R² < 0: Model worse than simply predicting mean
   • Note: High R² on synthetic data doesn't guarantee future performance; always test on real held-out data

3. Part 4 KDE Plots

   • Histogram: empirical distribution of MC samples (noisy but unbiased)
   • KDE curve: smooth estimate (may over/under-smooth depending on bandwidth)
   • Green line: last real price (reference)
   • Red line: mean simulated price (directional forecast)
   • Spread: uncertainty range

Forecast Confidence

Use KDE to compute tail probabilities:

# In a Jupyter notebook after running the pipeline:
import numpy as np
from scipy.integrate import quad

kde_func = part4_results['kde_funcs'][-1]  # Last horizon KDE
last_price = part1_results['prices'][-1]

# P(price > last_price)
prob_up, _ = quad(kde_func, last_price, np.inf)
print(f"Prob of increase: {prob_up:.2%}")

# 5% confidence interval
from scipy.optimize import fminbound
x_min = np.percentile(part4_results['sims'][:, -1], 2.5)
x_max = np.percentile(part4_results['sims'][:, -1], 97.5)
print(f"95% CI: [{x_min:.2f}, {x_max:.2f}]")

Customization

Increase forecast horizon:

python run_all.py --horizon 10  # 10-day forecast instead of 3

More simulations for smoother KDE:

python run_all.py --nsim 50000  # 50k paths instead of 10k

Try different KDE bandwidth:

python run_all.py --bandwidth silverman
python run_all.py --bandwidth 0.2

Different stock or date range:

python run_all.py --ticker MSFT --start 2023-01-01 --end 2023-12-31

---

Caveats & Limitations

1. Synthetic Data Overfitting: Part 3 trains on smooth fBB paths + real data. The high R² is partly due to unrealistic smoothness. Always validate on truly held-out future data.

2. No Model Uncertainty: XGBoost provides point estimates. Monte Carlo adds sampling noise but not epistemic (model) uncertainty. Consider ensembles for robustness.

3. Residual Bootstrap: Part 4 assumes residuals from training are representative of future noise. In regime change, this fails.

4. Lagged Features Only: Model uses only recent prices/vols. No sentiment, macroeconomic, or microstructure features.

5. No Transaction Costs or Slippage: Forecasts are "ideal" prices; real trading incurs costs.

6. Past Performance: Historical Hurst exponent and volatility may not persist. Always re-estimate on rolling windows in production.

---

Example Workflow

# 1. Forecast MSFT for the next 5 trading days with diagnostics
python run_all.py --ticker MSFT --horizon 5 --visualize

# 2. Output shows:
# - MSFT daily prices and volatility (Part 1)
# - Sample fBB paths + feature distributions (Part 2)
# - XGBoost test metrics + residual plots (Part 3)
# - KDE plots for days 1-5 ahead (Part 4)

# 3. Interpret:
# If Hurst < 0.5 and model predicts upward, market may be mean-reverting → take contrarian view
# If KDE spread is wide, uncertainty is high → smaller position size
# If mean simulated price far from last real price, model may be extrapolating → check residuals

---

References & Further Reading

• Hurst Exponent: Peters, E.E. (1991). Chaos and Order in the Capital Markets.
• fBM & Bridges: Nualart, D. (2006). The Malliavin Calculus and Related Topics.
• XGBoost: Chen, T. & Guestrin, C. (2016). "XGBoost: A Scalable Tree Boosting System". KDD.
• KDE: Silverman, B.W. (1986). Density Estimation for Statistics and Data Analysis.

---

License

MIT License (adjust as needed)

---

Questions or Issues?

Open an issue on GitHub or contact the author.