jepa-fx-risk/scripts/prepare_hourly.py

"""Prepare EUR/USD hourly OHLCV + realized vol from histdata M1 zips.

Aggregates all M1 bars in data/raw/DAT_ASCII_EURUSD_M1_*.zip to hourly.
Realized vol per hour = sqrt(sum(log-return²)) over the constituent M1 bars.
Weekend hours are naturally absent (FX market closed Sat/Sun); NO interpolation.
Hours with fewer than MIN_BARS M1 bars are dropped (holidays, thin sessions).

Output: data/processed/eurusd_hourly.parquet
Columns: datetime (UTC, tz-naive), close, ret (log), realized_vol

    python scripts/prepare_hourly.py
    RAW=data/raw OUT=data/processed/eurusd_hourly.parquet python scripts/prepare_hourly.py
"""
import glob
import os
import zipfile

import numpy as np
import pandas as pd

RAW_DEFAULT = "data/raw"
OUT_DEFAULT = "data/processed/eurusd_hourly.parquet"
MIN_BARS    = 30   # drop hours thinner than this (holidays, DST boundary artefacts)


# ── Core transformation ──────────────────────────────────────────────────────

def resample_to_hourly(m1: pd.DataFrame) -> pd.DataFrame:
    """Aggregate M1 DataFrame to hourly bars.

    Args:
        m1: DataFrame with columns ['ts', 'open', 'high', 'low', 'close']
            ('open'/'high'/'low' optional — omit for close-only data).

    Returns:
        DataFrame with columns ['datetime', 'close', 'ret', 'realized_vol',
        'hl_range', 'ret_intrabar'] sorted by datetime.
        Hours with fewer than MIN_BARS M1 ticks are dropped.
    """
    m1 = m1.sort_values("ts").copy()
    m1["log_r"] = np.log(m1["close"]).diff()
    m1["hour"]  = m1["ts"].dt.floor("h")

    has_ohlc = all(c in m1.columns for c in ("open", "high", "low"))

    agg_dict = dict(
        close        = ("close", "last"),
        realized_vol = ("log_r", lambda x: np.sqrt(np.nansum(x.values ** 2))),
        n_bars       = ("log_r", "count"),
    )
    if has_ohlc:
        agg_dict["high"]  = ("high",  "max")
        agg_dict["low"]   = ("low",   "min")
        agg_dict["open_"] = ("open",  "first")

    agg = m1.groupby("hour").agg(**agg_dict).reset_index()

    agg = agg[agg["n_bars"] >= MIN_BARS].copy()
    agg["ret"] = np.log(agg["close"]).diff()
    agg = agg.dropna(subset=["ret"]).reset_index(drop=True)
    agg = agg.rename(columns={"hour": "datetime"})

    if has_ohlc:
        agg["hl_range"]    = np.log(agg["high"] / agg["low"])
        agg["ret_intrabar"]= np.log(agg["close"] / agg["open_"])
        cols = ["datetime", "close", "ret", "realized_vol", "hl_range", "ret_intrabar"]
    else:
        cols = ["datetime", "close", "ret", "realized_vol"]

    return agg[cols]


def load_m1_from_zips(raw_dir: str) -> pd.DataFrame:
    """Load and concatenate all M1 zips from raw_dir (histdata format)."""
    pattern = os.path.join(raw_dir, "DAT_ASCII_EURUSD_M1_*.zip")
    zips = sorted(glob.glob(pattern))
    if not zips:
        raise FileNotFoundError(f"No M1 zips found at {pattern}")
    frames = []
    for zp in zips:
        with zipfile.ZipFile(zp) as z:
            csv = [n for n in z.namelist() if n.endswith(".csv")][0]
            with z.open(csv) as f:
                df = pd.read_csv(
                    f, sep=";", header=None,
                    names=["dt", "open", "high", "low", "close", "vol"],
                )
        df["ts"] = pd.to_datetime(df["dt"], format="%Y%m%d %H%M%S")
        frames.append(df[["ts", "open", "high", "low", "close"]])
        print(f"  loaded {os.path.basename(zp)}: {len(df):,} rows")
    return pd.concat(frames).sort_values("ts").reset_index(drop=True)


def build_hourly_parquet(
    raw_dir: str = RAW_DEFAULT,
    out_path: str = OUT_DEFAULT,
) -> pd.DataFrame:
    """Full pipeline: load all M1 zips → hourly parquet. Returns the DataFrame."""
    print(f"Loading M1 zips from {raw_dir}...")
    m1 = load_m1_from_zips(raw_dir)
    print(f"Total M1 bars: {len(m1):,}  ({m1['ts'].min().date()} → {m1['ts'].max().date()})")

    print("Resampling to hourly...")
    hourly = resample_to_hourly(m1)
    print(f"Hourly rows: {len(hourly):,}  ({hourly['datetime'].min()} → {hourly['datetime'].max()})")

    # Sanity: COVID crash (Mar 2020) should show realized vol spike if data covers it
    if hourly["datetime"].dt.year.isin([2020]).any():
        rv = hourly.set_index("datetime")["realized_vol"]
        try:
            mar20 = rv["2020-03-01":"2020-03-31"].max()
            typ   = rv["2019-01-01":"2019-12-31"].median()
            print(f"Sanity — median 2019 RV: {typ:.6f} | max Mar-2020 RV: {mar20:.6f} | spike ×{mar20/typ:.1f}")
        except Exception:
            pass

    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
    hourly.to_parquet(out_path, index=False)
    print(f"Written: {out_path}")
    return hourly


if __name__ == "__main__":
    raw_dir  = os.environ.get("RAW",  RAW_DEFAULT)
    out_path = os.environ.get("OUT",  OUT_DEFAULT)
    build_hourly_parquet(raw_dir=raw_dir, out_path=out_path)