Volatility Forecasting in Crypto Markets Using Neural Networks as Nonlinear Autoregressive Models

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt

from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error

# =========================
# Configuration
# =========================
TICKER = "ETH-USD"
START = "2020-01-01"
LAGS = 3
ANN_VOL_TARGET = 0.25
MAX_LEVERAGE = 2.0
REFIT_EVERY_DAYS = 7
TRAIN_MIN_YEARS = 1
TC_BPS_PER_TURNOVER = 10.0
OUT_DIR_PREFIX = "nnar_tv_"

# =========================
# Data and volatility
# =========================
def load_data(ticker, start):
    df = yf.download(ticker, start=start, auto_adjust=True, progress=False)
    if df.empty:
        raise ValueError("No data downloaded; check ticker/date.")
    return df

def garman_klass_vol(df):
    rs = 0.5 * (np.log(df['High'] / df['Low']))**2 - \
         (2*np.log(2)-1) * (np.log(df['Close'] / df['Open']))**2
    vol = np.sqrt(np.maximum(rs, 0))
    vol.name = "rv"
    return vol

def daily_log_returns(prices):
    r = np.log(prices).diff()
    r.name = "ret"
    return r

# =========================
# Features
# =========================
def build_feature_table(rv, ret, lags):
    df = pd.concat([rv, ret], axis=1)
    for k in range(1, lags+1):
        df[f"rv_lag_{k}"] = df["rv"].shift(k)
        df[f"ret_lag_{k}"] = df["ret"].shift(k)
    df = df.dropna()
    return df

# =========================
# NN-AR forecaster
# =========================
class NNARVolForecaster:
    def __init__(self, hidden=(32,16), max_iter=800, random_state=42):
        self.model = MLPRegressor(hidden_layer_sizes=hidden, activation="relu",
                                  solver="adam", max_iter=max_iter, random_state=random_state)
        self.scaler = StandardScaler()
        self.fitted = False

    def fit(self, X, y):
        Xs = self.scaler.fit_transform(X)
        self.model.fit(Xs, y)
        self.fitted = True

    def predict(self, X):
        if not self.fitted:
            raise RuntimeError("Model not fitted.")
        Xs = self.scaler.transform(X)
        return self.model.predict(Xs)

def walk_forward_forecast(features, first_train_end_idx):
    cols_X = [c for c in features.columns if c != "rv"]
    vhat = pd.Series(index=features.index, dtype=float)

    model = NNARVolForecaster()
    last_fit_loc = None

    for loc in range(first_train_end_idx, len(features)):
        if (last_fit_loc is None) or ((loc - last_fit_loc) >= REFIT_EVERY_DAYS):
            train = features.iloc[:loc]
            X_train = train[cols_X]
            y_train = train["rv"]
            model.fit(X_train, y_train)
            last_fit_loc = loc

        row = features.iloc[[loc]][cols_X]
        vhat.iloc[loc] = model.predict(row)[0]

    out = features.copy()
    out["vhat"] = vhat
    return out

# =========================
# Target-vol overlay
# =========================
def target_vol_weights(vhat_daily, ann_vol_target=0.10, max_leverage=1.5):
    sigma_target_daily = ann_vol_target / np.sqrt(252.0)
    w = sigma_target_daily / vhat_daily.replace(0, np.nan)
    w = w.clip(lower=0, upper=max_leverage).fillna(0.0)
    return w

def apply_transaction_costs(ret_series, weights, bps_per_turnover=0.0):
    turnover = weights.diff().abs().fillna(0.0)
    daily_cost = turnover * (bps_per_turnover / 1e4)
    return ret_series - daily_cost

# =========================
# Performance metrics
# =========================
def annualized_vol(returns):
    return returns.std() * np.sqrt(252.0)

def sharpe(returns, rf=0.0):
    ex = returns - rf/252.0
    vol = ex.std()
    return (ex.mean() / vol * np.sqrt(252.0)) if vol > 0 else np.nan

def max_drawdown(series):
    cum = (1 + series).cumprod()
    peak = cum.cummax()
    dd = (cum / peak - 1.0).min()
    return dd

# =========================
# Main
# =========================
def main():
    df = load_data(TICKER, START)
    rv = garman_klass_vol(df)
    ret = daily_log_returns(df["Close"])
    feat = build_feature_table(rv, ret, LAGS).dropna()

    first_train_end_idx = feat.index.get_indexer_for(
        [feat.index[0] + pd.DateOffset(years=TRAIN_MIN_YEARS)]
    )
    if len(first_train_end_idx) == 0 or first_train_end_idx[0] <= 0:
        raise ValueError("Not enough history for initial training window.")
    first_train_end_loc = first_train_end_idx[0]

    wf = walk_forward_forecast(feat, first_train_end_loc).dropna(subset=["vhat"])

    w = target_vol_weights(wf["vhat"], ANN_VOL_TARGET, MAX_LEVERAGE)
    strat_ret_gross = w * wf["ret"]
    strat_ret = apply_transaction_costs(strat_ret_gross, w, TC_BPS_PER_TURNOVER)
    bh_ret = wf["ret"]

    metrics = pd.DataFrame({
        "Total Return": [(1+bh_ret).prod()-1, (1+strat_ret).prod()-1],
        "Ann. Vol": [annualized_vol(bh_ret), annualized_vol(strat_ret)],
        "Sharpe": [sharpe(bh_ret), sharpe(strat_ret)],
        "Max Drawdown": [max_drawdown(bh_ret), max_drawdown(strat_ret)]
    }, index=["Buy & Hold", "NNAR Target-Vol"])

    print("\nVolatility Forecast MSE:", mean_squared_error(wf["rv"], wf["vhat"]))
    print("\nPerformance Summary:")
    print(metrics)

    eq_bh = (1 + bh_ret).cumprod()
    eq_st = (1 + strat_ret).cumprod()

    fig, ax = plt.subplots(figsize=(12,4))
    ax.plot(eq_bh.index, eq_bh.values, label="Buy & Hold")
    ax.plot(eq_st.index, eq_st.values, label="NNAR Target-Vol")
    ax.set_title("Equity Curves")
    ax.legend()
    ax.grid(True, alpha=0.3)
    fig.savefig(f"{OUT_DIR_PREFIX}{TICKER}_equity.png", dpi=150)

if __name__ == "__main__":
    main()