Feature Engineering Gold Volatility Signals That Actually Predict Price Moves

import numpy as np
from datetime import datetime

# Personal note: Learned this after bad data killed a live model
def load_gold_data(filepath):
    """Load 1-minute gold futures OHLCV data"""
    # Expected format: timestamp,open,high,low,close,volume
    data = np.genfromtxt(
        filepath, 
        delimiter=',', 
        skip_header=1,
        dtype=[('time', 'U19'), ('open', 'f8'), ('high', 'f8'),
               ('low', 'f8'), ('close', 'f8'), ('volume', 'i8')]
    )
    
    # Watch out: Missing data will propagate through all features
    assert not np.any(np.isnan(data['close'])), "Found NaN prices"
    assert np.all(data['high'] >= data['low']), "High < Low detected"
    
    return data

# Load your data
gold_data = load_gold_data('gold_1min_2024.csv')
print(f"Loaded {len(gold_data)} bars")
print(f"Date range: {gold_data['time'][0]} to {gold_data['time'][-1]}")

def engineer_price_features(data):
    """Build 5 core price features"""
    close = data['close']
    high = data['high']
    low = data['low']
    
    features = {}
    
    # 1. Log returns - percentage change in log space
    # Personal note: Log returns are additive and handle large moves better
    features['log_return'] = np.diff(np.log(close), prepend=np.log(close[0]))
    
    # 2. True Range - actual volatility including gaps
    prev_close = np.roll(close, 1)
    prev_close[0] = close[0]
    
    tr1 = high - low
    tr2 = np.abs(high - prev_close)
    tr3 = np.abs(low - prev_close)
    features['true_range'] = np.maximum(tr1, np.maximum(tr2, tr3))
    
    # 3. High-Low Range Ratio - intrabar volatility
    features['hl_ratio'] = (high - low) / close
    
    # 4. Close Position in Range - where did price settle?
    # 0 = closed at low, 1 = closed at high
    features['close_position'] = np.where(
        high != low,
        (close - low) / (high - low),
        0.5
    )
    
    # 5. Price acceleration - rate of change of returns
    features['return_accel'] = np.diff(
        features['log_return'], 
        prepend=features['log_return'][0]
    )
    
    return features

price_features = engineer_price_features(gold_data)
print(f"Built {len(price_features)} price features")
print(f"True Range mean: {price_features['true_range'].mean():.4f}")

def rolling_volatility_features(returns, true_range):
    """Calculate multi-timeframe volatility without lookahead"""
    features = {}
    
    # Windows: 5min, 15min, 60min, 240min (4hr)
    windows = [5, 15, 60, 240]
    
    for window in windows:
        # Standard deviation of returns
        features[f'vol_std_{window}'] = rolling_std(returns, window)
        
        # Average True Range - volatility with gaps
        features[f'vol_atr_{window}'] = rolling_mean(true_range, window)
        
        # Realized volatility - proper annualization
        # sqrt(252 * 390) = sqrt(trading mins per year)
        features[f'vol_realized_{window}'] = (
            rolling_std(returns, window) * np.sqrt(252 * 390)
        )
    
    return features

def rolling_mean(data, window):
    """Vectorized rolling mean - no lookahead"""
    # Personal note: cumsum trick is 40x faster than loop
    cumsum = np.cumsum(np.insert(data, 0, 0))
    result = (cumsum[window:] - cumsum[:-window]) / window
    
    # Pad start with expanding mean
    pad = np.array([np.mean(data[:i+1]) for i in range(window-1)])
    return np.concatenate([pad, result])

def rolling_std(data, window):
    """Vectorized rolling std - numerically stable"""
    # Watch out: naive variance calc fails with large numbers
    mean = rolling_mean(data, window)
    squared_diff = (data - mean) ** 2
    variance = rolling_mean(squared_diff, window)
    return np.sqrt(variance)

vol_features = rolling_volatility_features(
    price_features['log_return'],
    price_features['true_range']
)
print(f"Built {len(vol_features)} volatility features")

def volatility_regime_features(vol_features):
    """Detect volatility regime changes"""
    features = {}
    
    # Use 60-minute ATR as baseline volatility
    baseline_vol = vol_features['vol_atr_60']
    
    # 1. Volatility percentile - where are we in recent distribution?
    features['vol_percentile'] = rolling_percentile(baseline_vol, 240)
    
    # 2. Volatility regime - categorical but encoded numerically
    # Low: <33rd percentile, Medium: 33-66, High: >66
    features['vol_regime'] = np.digitize(
        features['vol_percentile'],
        bins=[0, 33, 66, 100]
    ) - 1
    
    # 3. Regime change indicator - did we just switch regimes?
    features['regime_change'] = np.abs(
        np.diff(features['vol_regime'], prepend=features['vol_regime'][0])
    )
    
    # 4. Volatility expansion - current vs recent average
    features['vol_expansion'] = (
        baseline_vol / rolling_mean(baseline_vol, 60)
    )
    
    # 5. Volatility persistence - autocorrelation of volatility
    features['vol_persistence'] = rolling_autocorr(baseline_vol, 60, lag=1)
    
    return features

def rolling_percentile(data, window):
    """Calculate rolling percentile rank"""
    result = np.zeros_like(data)
    
    for i in range(len(data)):
        start = max(0, i - window + 1)
        window_data = data[start:i+1]
        # Percentile rank of current value in window
        result[i] = (window_data < data[i]).sum() / len(window_data) * 100
    
    return result

def rolling_autocorr(data, window, lag):
    """Rolling autocorrelation at specific lag"""
    result = np.zeros_like(data)
    
    for i in range(lag, len(data)):
        start = max(0, i - window + 1)
        x = data[start:i+1]
        y = data[start-lag:i+1-lag] if start >= lag else data[0:i+1-lag]
        
        if len(x) > 1 and len(y) > 1:
            result[i] = np.corrcoef(x[:len(y)], y)[0, 1]
    
    return result

regime_features = volatility_regime_features(vol_features)
print(f"Built {len(regime_features)} regime features")
print(f"High vol regime: {(regime_features['vol_regime'] == 2).sum()} bars")

def create_feature_matrix(price_feats, vol_feats, regime_feats):
    """Combine all features into single matrix"""
    # Stack all features
    all_features = {**price_feats, **vol_feats, **regime_feats}
    
    # Convert to 2D array
    n_samples = len(gold_data)
    n_features = len(all_features)
    
    feature_matrix = np.zeros((n_samples, n_features))
    feature_names = []
    
    for i, (name, values) in enumerate(all_features.items()):
        feature_matrix[:, i] = values
        feature_names.append(name)
    
    return feature_matrix, feature_names

def normalize_features(X, method='robust'):
    """Normalize features - robust to outliers"""
    if method == 'robust':
        # Use median and IQR instead of mean and std
        median = np.median(X, axis=0)
        q75 = np.percentile(X, 75, axis=0)
        q25 = np.percentile(X, 25, axis=0)
        iqr = q75 - q25
        
        # Avoid division by zero
        iqr = np.where(iqr == 0, 1, iqr)
        
        X_normalized = (X - median) / iqr
    else:
        # Standard z-score
        mean = np.mean(X, axis=0)
        std = np.std(X, axis=0)
        std = np.where(std == 0, 1, std)
        X_normalized = (X - mean) / std
    
    return X_normalized

# Build final feature set
X, feature_names = create_feature_matrix(
    price_features, vol_features, regime_features
)

# Normalize
X_normalized = normalize_features(X, method='robust')

print(f"Feature matrix shape: {X_normalized.shape}")
print(f"Features: {', '.join(feature_names[:5])}...")
print(f"No NaN values: {not np.any(np.isnan(X_normalized))}")