Cryptocurrency Volatility: A Complete Time Series Pipeline

1. Research Setup — Three Questions, One Dataset

This hands-on module is structured as a research investigation. Instead of learning techniques in isolation, you will answer three specific research questions about Bitcoin's price behaviour. Each question is answered by a specific statistical model:

The Dataset

We analyse BTCUSD daily closing prices — 301 trading days from August 2025 to June 2026. Bitcoin dropped from ~$111,000 to ~$64,000 during this period, providing ample volatility for our models to capture.

2. Data Loading and Exploratory Analysis

Step 1: Load and Prepare the Data

# =====================================================================
# HANDS-ON: ARIMA + ARCH + GARCH Research Pipeline
# Dataset: BTCUSD Daily (301 trading days, Aug 2025 – Jun 2026)
# =====================================================================
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

plt.style.use('seaborn-v0_8-darkgrid')

# ── Load CSV ──────────────────────────────────────────────────────
# TO USE YOUR OWN DATA: Change this file path
df = pd.read_csv('../../assets/datasets/btcusd_daily.csv', parse_dates=['date'])
df = df.set_index('date').sort_index()

# Display basic information
print(f"Observations: {len(df)}")
print(f"Period: {df.index[0].strftime('%Y-%m-%d')} to {df.index[-1].strftime('%Y-%m-%d')}")
print(f"Columns: {list(df.columns)}")
print(f"\nFirst 3 rows:")
print(df.head(3))
            

Observations: 301
Period: 2025-08-27 to 2026-06-23
Columns: ['open', 'high', 'low', 'close', 'volume']

First 3 rows:
              open     high      low    close    volume
date
2025-08-27  111720.5  112581.0  110327.0  111206.5  4971592.0
2025-08-28  111200.0  113440.5  110822.0  112519.0  4301404.0
2025-08-29  112519.0  112597.5  107426.0  108329.0  5922968.0

Step 2: Compute Daily Returns

We work with log returns (continuously compounded) rather than raw prices. Returns are stationary — prices are not. The formula: r_t = ln(P_t / P_{t-1}) × 100

# ── Compute returns ───────────────────────────────────────────────
# TO USE YOUR OWN DATA: Change 'close' to your price column name
df['return'] = np.log(df['close'] / df['close'].shift(1)) * 100
ret = df['return'].dropna()  # Remove first row (NaN)

print(f"Return observations: {len(ret)}")
print(f"Mean daily return:   {ret.mean():.4f}%")
print(f"Std daily return:    {ret.std():.4f}%")
print(f"Min return:          {ret.min():.4f}%")
print(f"Max return:          {ret.max():.4f}%")
print(f"Skewness:            {ret.skew():.4f}")
print(f"Excess Kurtosis:     {ret.kurtosis():.4f}")
            

Return observations: 300
Mean daily return:   -0.1570%
Std daily return:    2.3543%
Min return:          -12.4257%
Max return:          9.0721%
Skewness:            -0.2776
Excess Kurtosis:     6.2263

Step 3: Visual Inspection — Price, Returns, and Volatility Clustering

# ── Plot price and returns ────────────────────────────────────────
fig, axes = plt.subplots(3, 1, figsize=(14, 10))

# Panel 1: BTC Price
axes[0].plot(df.index, df['close'], linewidth=1.2, color='#f7931a')
axes[0].set_title('BTCUSD Closing Price', fontweight='bold', fontsize=13)
axes[0].set_ylabel('Price (USD)')

# Panel 2: Daily Returns
axes[1].plot(ret.index, ret.values, linewidth=0.5, color='steelblue')
axes[1].axhline(y=0, color='gray', linestyle='--', linewidth=0.5)
axes[1].set_title('Daily Log Returns (%)', fontweight='bold', fontsize=13)
axes[1].set_ylabel('Return (%)')

# Panel 3: Squared Returns (proxy for volatility)
axes[2].fill_between(ret.index, ret.values**2, alpha=0.6, color='red')
axes[2].set_title('Squared Returns — Volatility Proxy', fontweight='bold', fontsize=13)
axes[2].set_ylabel('Squared Return')
axes[2].set_xlabel('Date')

plt.tight_layout()
plt.savefig('../../assets/images/hands-on-btc-exploration.png', dpi=150, bbox_inches='tight')
plt.show()
            

What We Observe

• Price panel: Bitcoin declined from ~$111K to ~$64K — a ~42% drawdown over 10 months
• Returns panel: No obvious trend or seasonality — returns oscillate around zero, suggesting stationarity
• Squared returns panel: Clear clusters of high volatility (large squared returns bunched together). This is the volatility clustering that ARCH/GARCH is designed to model
• Negative skewness (-0.28): More large negative days than large positive days — typical for equities and crypto
• High kurtosis (6.23): Extreme events occur far more often than a normal distribution predicts. The worst day was -12.4% — a 5+ sigma event under normality

These observations motivate our three research questions. The returns look unpredictable (RQ1 → ARIMA), but their magnitude clearly clusters (RQ2 → ARCH), and the clusters appear persistent (RQ3 → GARCH). Let us now test each question formally.

3. RQ1: Can Past Bitcoin Returns Predict Future Returns?

The Research Question

Hypothesis: Bitcoin daily returns contain predictable patterns (momentum, mean-reversion, or cycles).
Null Hypothesis: Returns follow a random walk — they are unpredictable (efficient market).
Test: Fit an ARIMA model. If auto_arima selects ARIMA(0,0,0), the best prediction is just the historical mean — returns are white noise, and we cannot reject market efficiency.

Step 1: Stationarity Test (ADF)

ARIMA requires the series to be stationary. Returns should already be stationary (unlike prices). We confirm this with the ADF test:

from statsmodels.tsa.stattools import adfuller

# ── ADF Test on Returns ───────────────────────────────────────────
adf_stat, adf_p, _, _, cv, _ = adfuller(ret)
print(f"ADF Statistic:     {adf_stat:.4f}")
print(f"p-value:           {adf_p:.6f}")
print(f"5% Critical Value: {cv['5%']:.4f}")
print(f"Conclusion: {'STATIONARY ✓' if adf_p < 0.05 else 'NON-STATIONARY — difference needed'}")
            

ADF Statistic:     -17.8514
p-value:           0.000000
5% Critical Value: -2.8713
Conclusion: STATIONARY ✓

The ADF test decisively rejects the unit root null. Here is how we reached that conclusion:

The logic: H₀ = "returns are non-stationary (have a unit root)." The ADF statistic (-17.85) is far below the 5% critical value (-2.87). If returns truly had a unit root, we would expect the ADF statistic to be near zero or positive. A value of -17.85 is extremely unlikely under H₀ — the p-value is essentially zero. Therefore, we reject H₀ and conclude returns are stationary. No differencing needed — d = 0.

Step 2: Identify p and q Using ACF/PACF

Before running auto_arima, it is important to understand what p and q actually mean and how the ACF and PACF plots help identify them manually. Even though we will use automatic selection, understanding the logic prevents blind trust in the algorithm.

What Are p, d, and q?

How to Read the ACF and PACF Plots

Each vertical bar represents the correlation at that lag. The blue shaded band is the 95% confidence interval — bars that extend outside the blue band are statistically significant (p < 0.05). Bars inside the blue band are not different from zero.

The identification rules:

• For AR(p): The PACF shows significant spikes at lags 1 through p, then drops to zero (inside the blue band) for all lags > p. Count the last significant PACF spike — that is your p.
• For MA(q): The ACF shows significant spikes at lags 1 through q, then drops to zero (inside the blue band) for all lags > q. Count the last significant ACF spike — that is your q.

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# ── ACF and PACF of returns ───────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
plot_acf(ret, lags=20, ax=axes[0])
axes[0].set_title('ACF of BTC Returns — Identifies q (MA order)', fontweight='bold')
plot_pacf(ret, lags=20, ax=axes[1], method='ywm')
axes[1].set_title('PACF of BTC Returns — Identifies p (AR order)', fontweight='bold')
plt.tight_layout()
plt.savefig('../../assets/images/hands-on-acf-pacf.png', dpi=150, bbox_inches='tight')
plt.show()
            

What We See in BTC Returns

Both plots suggest the same thing: Bitcoin daily returns are white noise. There is no significant autocorrelation at any lag. If we were identifying p and q manually, we would choose p = 0, q = 0 — meaning ARIMA(0,0,0), which is just Y_t = constant + ε_t (the mean plus random noise). But let us confirm this with the formal automatic selection procedure.

Step 3: Automatic ARIMA Selection — Let the Algorithm Search

How auto_arima works: It fits dozens of candidate ARIMA(p,d,q) models — trying different combinations of p (0 to max_p), d (fixed at 0 for stationary data), and q (0 to max_q) — and computes the AIC for each one. It then selects the model with the lowest AIC. This automates what we just did manually with the ACF/PACF plots, but more rigorously: instead of eyeballing the plots, it fits every candidate and compares their AIC values.

from pmdarima import auto_arima
from statsmodels.tsa.arima.model import ARIMA

# ── Train-test split (last 30 days for forecast evaluation) ───────
train = ret.iloc[:-30]   # First 270 days
test  = ret.iloc[-30:]   # Last 30 days (holdout)

# ── Auto-select best ARIMA ────────────────────────────────────────
# Searches p=0..5, q=0..5, d=0 (already stationary)
# For each candidate, fits the model and computes AIC
# Returns the order with the LOWEST AIC
auto_model = auto_arima(train, seasonal=False, trace=False,
                        stepwise=True, max_p=5, max_q=5, d=0,
                        information_criterion='aic')

print(f"Best ARIMA order: {auto_model.order}")
print(f"AIC: {auto_model.aic():.2f}")
            

Best ARIMA order: (0, 0, 0)
AIC: 1233.57

auto_arima confirms what the ACF/PACF plots suggested: ARIMA(0,0,0) is the best model. No AR terms (p=0), no differencing (d=0), no MA terms (q=0). The model is simply:

Y_t = c + ε_t

In plain English: Tomorrow's Bitcoin return = the historical average return + random noise. The model says you cannot do better than guessing the average, and even the average is not reliably different from zero.

This means Bitcoin daily returns have no predictable structure whatsoever. The best forecast of tomorrow's return is simply the average of all past returns — and even that average is not statistically different from zero.

Best ARIMA order: (0, 0, 0)
AIC: 1233.57

─── ARIMA(0,0,0) Model Summary ───
                               SARIMAX Results
==============================================================================
Dep. Variable:                 return   No. Observations:                  270
Model:                          ARIMA   Log Likelihood                -615.508
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.1077      0.145     -0.745      0.456      -0.391       0.176
sigma2         5.5925      0.234     23.923      0.000       5.134       6.051
===================================================================================
Ljung-Box (L1) (Q):                   1.49   Jarque-Bera (JB):               486.83
Prob(Q):                              0.22   Prob(JB):                         0.00
===================================================================================

Step 4: Forecast from ARIMA(0,0,0)

An ARIMA(0,0,0) forecast is straightforward: every future prediction equals the historical mean, with prediction intervals based on the residual standard deviation. Let us produce a 30-day forecast and evaluate it against the held-out test data:

# ── 30-day forecast ───────────────────────────────────────────────
forecast_result = arima_model.get_forecast(steps=30)
forecast_mean = forecast_result.predicted_mean
forecast_ci = forecast_result.conf_int(alpha=0.05)

# Calculate forecast accuracy on test set
from sklearn.metrics import mean_squared_error, mean_absolute_error
rmse = np.sqrt(mean_squared_error(test, forecast_mean))
mae  = mean_absolute_error(test, forecast_mean)

print(f"─── 30-Day Forecast Evaluation ───")
print(f"RMSE: {rmse:.4f}%")
print(f"MAE:  {mae:.4f}%")
print(f"(The forecast is flat at {forecast_mean[0]:.4f}% — the historical mean)")
print(f"Test set mean: {test.mean():.4f}%")
            

─── 30-Day Forecast Evaluation ───
RMSE: 2.2530%
MAE:  1.6923%
(The forecast is flat at -0.1077% — the historical mean)
Test set mean: -0.0534%

The RMSE of 2.25% is essentially the standard deviation of Bitcoin returns — the model predicts "average" every day, and the errors are just the natural daily volatility. An RMSE close to the return standard deviation (2.35%) confirms that the model has no predictive power beyond the mean. The flat forecast is the honest answer: you cannot predict Bitcoin's direction.

4. RQ2: Does Yesterday's Volatility Predict Today's Volatility?

The Research Question

Hypothesis: Bitcoin exhibits volatility clustering — large price moves tend to be followed by more large moves.
Null Hypothesis: Volatility is constant over time (homoscedastic — variance never changes).
Test: ARCH-LM test on squared returns. If p < 0.05, volatility clustering exists. Then fit ARCH(1) to estimate how much of yesterday's shock carries into today's variance.

Step 1: Formal ARCH-LM Test

import statsmodels.api as sm
from scipy.stats import chi2

# ── ARCH-LM Test: Regress squared returns on 10 lags ──────────────
sq_ret = ret.values ** 2
y = sq_ret[10:]  # Dependent: squared returns from t=11 onward
X = np.column_stack([sq_ret[10-i-1 : len(sq_ret)-i-1] for i in range(10)])
X = sm.add_constant(X)

arch_lm_model = sm.OLS(y, X).fit()
n = len(y)
lm_stat = n * arch_lm_model.rsquared
lm_pval = 1 - chi2.cdf(lm_stat, 10)

print(f"─── ARCH-LM Test (10 lags) ───")
print(f"LM Statistic: {lm_stat:.1f}")
print(f"p-value:      {lm_pval:.6f}")
print(f"Conclusion:   {'ARCH EFFECTS PRESENT — volatility is predictable' if lm_pval < 0.05 else 'No ARCH effects'}")
            

─── ARCH-LM Test (10 lags) ───
LM Statistic: 54.2
p-value:      0.000000
Conclusion:   ARCH EFFECTS PRESENT — volatility is predictable

What the ARCH-LM Test Actually Tells Us

The ARCH-LM result (LM = 54.2, p < 0.000001) confirms three interconnected facts. Here they are, explained from the ground up:

1. "Squared returns are autocorrelated" — What does this mean?
Take Bitcoin's daily returns: +2.3%, -1.5%, +0.8%, -4.2%, +3.1%, -3.8%... Now square them: 5.29, 2.25, 0.64, 17.64, 9.61, 14.44... Notice how large squared values (17.64) tend to be near other large values (9.61, 14.44), and small squared values (0.64, 2.25) tend to be near other small values. The squared returns are correlated with their own past — a large squared return today makes a large squared return tomorrow more likely. This is autocorrelation of squared returns. It means the magnitude of moves clusters, even though the direction does not.

2. Why this contradicts OLS
Standard OLS regression assumes homoscedasticity — the error variance is constant across all observations. In Module 1, we tested this with Breusch-Pagan. If variance is constant, then knowing yesterday's squared return tells you nothing about today's variance. But the ARCH-LM test proves it DOES. The OLS assumption is violated. This means OLS standard errors, t-statistics, and p-values are all unreliable for this data. You need models (ARCH/GARCH) that explicitly allow variance to change over time.

3. "Bitcoin's risk is time-varying and predictable" — What does this mean practically?
On Monday, you can look at last week's price movements and estimate whether this week will be calm or turbulent. If last week had ±4% daily swings, this week is likely to also have large swings. If last week was quiet (±1%), this week is likely quiet too. You cannot predict direction (RQ1 proved that), but you CAN predict risk level. This is the difference between predicting a coin flip's outcome (impossible) and predicting whether the coin is being flipped in a hurricane or a calm room (possible, and useful).

Step 2: Fit ARCH(1) Model

ARCH(1) models the conditional variance as: σ²_t = ω + α · ε²_{t-1}. The parameter α tells us what fraction of yesterday's squared shock transmits to today's variance.

from arch import arch_model

# ── ARCH(1) ───────────────────────────────────────────────────────
arch_res = arch_model(ret, vol='ARCH', p=1).fit(disp=False)
print("─── ARCH(1) Model ───")
print(arch_res.summary())

# Extract key parameters
omega = arch_res.params['omega']
alpha = arch_res.params['alpha[1]']
print(f"\n─── ARCH(1) Key Parameters ───")
print(f"ω (omega) = {omega:.4f}  — baseline daily variance")
print(f"α (alpha) = {alpha:.4f}  — shock transmission rate")
print(f"Baseline daily volatility (σ) = {np.sqrt(omega):.4f}%")
print(f"If yesterday's shock was ±3%: today's σ² = {omega:.4f} + {alpha:.4f}×9 = {omega + alpha*9:.4f}")
print(f"  → today's σ = {np.sqrt(omega + alpha*9):.4f}%")
            

─── ARCH(1) Model ───
                      Constant Mean - ARCH Model Results
==============================================================================
Dep. Variable:                 return   R-squared:                       0.000
Mean Model:             Constant Mean   Adj. R-squared:                  0.000
Vol Model:                       ARCH   Log-Likelihood:               -670.385
Distribution:                  Normal   AIC:                           1346.77
==============================================================================
                 coef    std err          t      P>|t|  95.0% Conf. Int.
------------------------------------------------------------------------------
mu            -0.0177      0.190     -0.093      0.926 [ -0.389,  0.354]
omega          4.3219      0.607      7.123  1.057e-12 [  3.133,  5.511]
alpha[1]       0.1954      0.204      0.958      0.338 [ -0.204,  0.595]
==============================================================================

─── ARCH(1) Key Parameters ───
ω (omega) = 4.3219  — baseline daily variance
α (alpha) = 0.1954  — shock transmission rate
Baseline daily volatility (σ) = 2.0789%
If yesterday's shock was ±3%: today's σ² = 4.3219 + 0.1954×9 = 6.0805
  → today's σ = 2.4659%

5. RQ3: How Persistent Is Bitcoin Volatility After a Crash?

The Research Question

Hypothesis: Bitcoin volatility has long memory — after a crash, elevated risk persists for days or weeks, not hours.
Test: Fit GARCH(1,1): σ²_t = ω + α·ε²_{t-1} + β·σ²_{t-1}. Compute persistence α+β and half-life.

Step 1: Fit GARCH(1,1)

# ── GARCH(1,1) ────────────────────────────────────────────────────
garch_res = arch_model(ret, vol='GARCH', p=1, q=1).fit(disp=False)
print("─── GARCH(1,1) Model ───")
print(garch_res.summary())

# Extract key parameters
gparams = garch_res.params
omega_g = gparams['omega']
alpha_g = gparams['alpha[1]']
beta_g  = gparams['beta[1]']
persistence = alpha_g + beta_g
half_life = np.log(0.5) / np.log(persistence)

print(f"\n─── GARCH(1,1) Key Parameters ───")
print(f"ω  = {omega_g:.4f}   (baseline daily variance)")
print(f"α  = {alpha_g:.4f}   (shock transmission)")
print(f"β  = {beta_g:.4f}    (variance persistence)")
print(f"α+β = {persistence:.4f} (total persistence)")
print(f"Half-life = {half_life:.1f} days ({half_life/7:.1f} weeks)")
            

─── GARCH(1,1) Model ───
                     Constant Mean - GARCH Model Results
==============================================================================
Dep. Variable:                 return   R-squared:                       0.000
Model:                         GARCH   Log-Likelihood:               -663.900
Vol Model:                      GARCH   AIC:                           1335.80
Distribution:                  Normal   BIC:                           1350.61
==============================================================================
                 coef    std err          t      P>|t|  95.0% Conf. Int.
------------------------------------------------------------------------------
mu            -0.0234      0.126     -0.186      0.852 [ -0.270,  0.223]
omega          0.6489      0.864      0.751      0.452 [ -1.044,  2.342]
alpha[1]       0.1612      0.105      1.541      0.123 [-4.381e-02,  0.366]
beta[1]        0.7264      0.229      3.176  1.492e-03 [  0.278,  1.175]
==============================================================================

─── GARCH(1,1) Key Parameters ───
ω  = 0.6489   (baseline daily variance)
α  = 0.1612   (shock transmission)
β  = 0.7264   (variance persistence)
α+β = 0.8876 (total persistence)
Half-life = 5.8 days (0.8 weeks)

Step 2: Model Comparison — ARCH(1) vs GARCH(1,1)

6. Research Summary and Using Your Own Data

What We Found — Three Research Questions Answered

Research Methods Paragraph (Copy-Paste Ready)

The following is a template methods section suitable for a research paper or thesis:

METHODS — Time Series Analysis of Cryptocurrency Returns

We analyse daily Bitcoin (BTCUSD) closing prices from August 2025 to
June 2026 (N = 301 trading days). Log returns are computed as
r_t = ln(P_t / P_{t-1}) × 100.

Return Predictability: An augmented Dickey-Fuller test confirms returns
are stationary (ADF = -17.85, p < 0.001). Automatic ARIMA selection via
AIC minimisation identifies ARIMA(0,0,0) as the best-fitting model,
indicating that Bitcoin daily returns follow a white noise process
(Ljung-Box Q(1) = 1.49, p = 0.22). The constant term is not
significantly different from zero (μ = -0.11, p = 0.456), consistent
with the efficient market hypothesis.

Volatility Clustering: An ARCH-LM test on squared returns strongly
rejects the null of constant variance (LM = 54.2, p < 0.001),
confirming significant volatility clustering.

Volatility Persistence: A GARCH(1,1) model reveals moderate volatility
persistence (α+β = 0.888). The GARCH β coefficient (0.726, p = 0.001)
dominates the ARCH α coefficient (0.161, p = 0.123), indicating that
volatility is driven primarily by its own history rather than by
individual shocks. The half-life of a volatility shock is approximately
5.8 trading days. GARCH(1,1) outperforms ARCH(1) on both AIC
(1335.80 vs 1346.77) and log-likelihood (-663.90 vs -670.39),
confirming that the persistence term substantially improves model fit.
All models were estimated in Python using statsmodels and the arch
library.

Using Your Own Data — Template Guide

This entire module is designed as a reusable template. To analyse your own stock, crypto, forex, or any financial time series, change only these lines:

Your CSV must have:

• A date column (any format that pandas can parse)
• A price column (closing price — the code computes returns automatically)
• At least 100 observations for reliable ARIMA/GARCH estimation
• Rows sorted by date (oldest first — the code sorts automatically)