rusty_backtest/features/

market_impact.rs

//! Market Impact Analysis
//!
//! Implements Kyle's Lambda and other market impact metrics for measuring
//! the price impact of trades and order flow.

use super::{TradeSide, TradeTick, decimal_to_f64_or_nan};
use smallvec::SmallVec;

/// Kyle's Lambda calculator for market impact coefficient
///
/// Kyle's Lambda (λ) measures the price impact per unit of net order flow.
/// Formula: λ = Δp / Δq
/// where Δp is the price change and Δq is the net volume imbalance
#[derive(Debug, Clone)]
pub struct KylesLambdaCalculator {
    /// Window size for rolling calculations
    window_size: usize,
    /// Price changes buffer (circular)
    price_changes: Vec<f64>,
    /// Volume imbalances buffer (circular)
    volume_imbalances: Vec<f64>,
    /// Current position in circular buffer
    buffer_position: usize,
    /// Number of samples collected
    samples_count: usize,
    /// Previous price for calculating changes
    previous_price: Option<f64>,
    /// Accumulated buy volume in current interval
    interval_buy_volume: f64,
    /// Accumulated sell volume in current interval
    interval_sell_volume: f64,
    /// Time of last interval calculation (nanoseconds)
    last_interval_time: u64,
    /// Interval duration in nanoseconds
    interval_duration_ns: u64,
}

impl KylesLambdaCalculator {
    /// Create a new Kyle's Lambda calculator
    ///
    /// # Arguments
    /// * `window_size` - Number of intervals to use for rolling calculation
    /// * `interval_duration_ms` - Duration of each interval in milliseconds
    #[must_use]
    pub fn new(window_size: usize, interval_duration_ms: u64) -> Self {
        Self {
            window_size,
            price_changes: vec![0.0; window_size],
            volume_imbalances: vec![0.0; window_size],
            buffer_position: 0,
            samples_count: 0,
            previous_price: None,
            interval_buy_volume: 0.0,
            interval_sell_volume: 0.0,
            last_interval_time: 0,
            interval_duration_ns: interval_duration_ms * 1_000_000,
        }
    }

    /// Update with a new trade
    pub fn update_trade(&mut self, trade: &TradeTick) -> Option<f64> {
        let timestamp = trade.timestamp_ns;
        let price = decimal_to_f64_or_nan(trade.price);
        let quantity = decimal_to_f64_or_nan(trade.quantity);

        // Initialize timing on first trade
        if self.last_interval_time == 0 {
            self.last_interval_time = timestamp;
            self.previous_price = Some(price);
        }

        // Accumulate volume based on trade side
        if trade.side == TradeSide::Buy {
            self.interval_buy_volume += quantity;
        } else {
            self.interval_sell_volume += quantity;
        }

        // Check if interval has elapsed
        if timestamp >= self.last_interval_time + self.interval_duration_ns
            && let Some(prev_price) = self.previous_price
        {
            // Calculate price change
            let price_change = price - prev_price;

            // Calculate net volume imbalance
            let volume_imbalance = self.interval_buy_volume - self.interval_sell_volume;

            // Store in circular buffer
            self.price_changes[self.buffer_position] = price_change;
            self.volume_imbalances[self.buffer_position] = volume_imbalance;

            // Update buffer position
            self.buffer_position = (self.buffer_position + 1) % self.window_size;
            self.samples_count = self.samples_count.saturating_add(1).min(self.window_size);

            // Reset for next interval
            self.interval_buy_volume = 0.0;
            self.interval_sell_volume = 0.0;
            self.last_interval_time = timestamp;
            self.previous_price = Some(price);

            // Calculate Kyle's Lambda if we have enough samples
            return self.calculate_lambda();
        }

        None
    }
    /// Alias for update_trade() to maintain backward compatibility
    pub fn add_trade(&mut self, trade: &TradeTick) {
        self.update_trade(trade);
    }

    /// Alias for calculate_lambda() to maintain backward compatibility
    pub fn calculate(&self) -> f64 {
        self.calculate_lambda().unwrap_or(0.0)
    }

    /// Calculate Kyle's Lambda using linear regression
    ///
    /// Performs OLS regression: Δp = λ * Δq + ε
    /// Returns the slope coefficient λ
    fn calculate_lambda(&self) -> Option<f64> {
        if self.samples_count < 10 {
            return None; // Need minimum samples for reliable estimate
        }

        let n = self.samples_count;

        // Calculate means
        let mut sum_x = 0.0;
        let mut sum_y = 0.0;

        for i in 0..n {
            sum_x += self.volume_imbalances[i];
            sum_y += self.price_changes[i];
        }

        let mean_x = sum_x / n as f64;
        let mean_y = sum_y / n as f64;

        // Calculate slope using least squares
        let mut numerator = 0.0;
        let mut denominator = 0.0;

        for i in 0..n {
            let x_diff = self.volume_imbalances[i] - mean_x;
            let y_diff = self.price_changes[i] - mean_y;
            numerator += x_diff * y_diff;
            denominator += x_diff * x_diff;
        }

        if denominator.abs() < 1e-10 {
            return None; // Avoid division by zero
        }

        Some(numerator / denominator)
    }

    /// Get current Kyle's Lambda estimate
    pub fn get_lambda(&self) -> Option<f64> {
        self.calculate_lambda()
    }

    /// Reset the calculator
    pub fn reset(&mut self) {
        self.price_changes.fill(0.0);
        self.volume_imbalances.fill(0.0);
        self.buffer_position = 0;
        self.samples_count = 0;
        self.previous_price = None;
        self.interval_buy_volume = 0.0;
        self.interval_sell_volume = 0.0;
        self.last_interval_time = 0;
    }
}

/// Advanced market impact metrics
#[derive(Debug, Clone)]
pub struct MarketImpactAnalyzer {
    /// Kyle's Lambda calculator
    kyles_lambda: KylesLambdaCalculator,
    /// Temporary impact decay estimator
    temp_impact_halflife: f64,
    /// Permanent impact ratio
    permanent_impact_ratio: f64,
    /// Recent impacts for decay modeling
    recent_impacts: SmallVec<[(u64, f64, f64); 32]>, // (timestamp, impact, volume)
}

impl MarketImpactAnalyzer {
    /// Create a new market impact analyzer
    #[must_use]
    pub fn new(window_size: usize, interval_duration_ms: u64) -> Self {
        Self {
            kyles_lambda: KylesLambdaCalculator::new(window_size, interval_duration_ms),
            temp_impact_halflife: 300_000.0, // 5 minutes in milliseconds
            permanent_impact_ratio: 0.3,     // 30% of impact is permanent
            recent_impacts: SmallVec::new(),
        }
    }

    /// Update with new trade and get current impact metrics
    pub fn update(&mut self, trade: &TradeTick) -> MarketImpactMetrics {
        // Update Kyle's Lambda
        let kyles_lambda = self.kyles_lambda.update_trade(trade);

        // Calculate instantaneous impact if we have lambda
        let instantaneous_impact = if let Some(lambda) = kyles_lambda {
            let signed_volume = if trade.side == TradeSide::Buy {
                decimal_to_f64_or_nan(trade.quantity)
            } else {
                -decimal_to_f64_or_nan(trade.quantity)
            };
            lambda * signed_volume
        } else {
            0.0
        };

        // Add to recent impacts
        if instantaneous_impact.abs() > 1e-10 {
            self.recent_impacts.push((
                trade.timestamp_ns,
                instantaneous_impact,
                decimal_to_f64_or_nan(trade.quantity),
            ));

            // Keep only recent impacts (last hour)
            let cutoff_time = trade.timestamp_ns.saturating_sub(3_600_000_000_000);
            self.recent_impacts.retain(|(ts, _, _)| *ts > cutoff_time);
        }

        // Calculate cumulative temporary impact with decay
        let cumulative_temp_impact = self.calculate_cumulative_impact(trade.timestamp_ns);

        MarketImpactMetrics {
            kyles_lambda,
            instantaneous_impact,
            cumulative_temp_impact,
            permanent_impact: instantaneous_impact * self.permanent_impact_ratio,
        }
    }

    /// Calculate cumulative temporary impact with exponential decay
    fn calculate_cumulative_impact(&self, current_time: u64) -> f64 {
        let mut total_impact = 0.0;
        let decay_rate = 0.693 / self.temp_impact_halflife; // ln(2) / halflife

        for &(timestamp, impact, _) in &self.recent_impacts {
            let time_elapsed = (current_time - timestamp) as f64 / 1_000_000.0; // to milliseconds
            let decay_factor = (-decay_rate * time_elapsed).exp();
            total_impact += impact * decay_factor * (1.0 - self.permanent_impact_ratio);
        }

        total_impact
    }

    /// Reset all calculators
    pub fn reset(&mut self) {
        self.kyles_lambda.reset();
        self.recent_impacts.clear();
    }
}

/// Market impact metrics
#[derive(Debug, Clone, Copy)]
pub struct MarketImpactMetrics {
    /// Kyle's Lambda coefficient
    pub kyles_lambda: Option<f64>,
    /// Instantaneous price impact of the last trade
    pub instantaneous_impact: f64,
    /// Cumulative temporary impact (with decay)
    pub cumulative_temp_impact: f64,
    /// Permanent impact component
    pub permanent_impact: f64,
}

/// SIMD-optimized batch Kyle's Lambda calculation using safe operations
pub mod simd {
    use wide::f64x4;

    /// Calculate Kyle's Lambda for multiple symbols in parallel using safe SIMD
    ///
    /// This function processes 4 symbols at a time using portable `wide` operations.
    /// Completely safe, NaN-aware, and works on all platforms.
    ///
    /// # Arguments
    /// * `price_changes` - Slice of price change arrays, one per symbol
    /// * `volume_imbalances` - Slice of volume imbalance arrays, one per symbol
    ///
    /// # Returns
    /// Vector of Kyle's Lambda coefficients, one per symbol
    pub fn calculate_kyles_lambda_batch(
        price_changes: &[&[f64]],
        volume_imbalances: &[&[f64]],
    ) -> Vec<f64> {
        let n_symbols = price_changes.len();
        let mut results = vec![0.0; n_symbols];

        // Process 4 symbols at a time using safe SIMD
        let chunks = n_symbols / 4;

        for chunk in 0..chunks {
            let base_idx = chunk * 4;

            // Ensure all 4 symbols have data
            if base_idx + 3 >= n_symbols {
                break;
            }

            // Get sample count (assume all symbols have same length)
            let n_samples = price_changes[base_idx].len();
            if n_samples == 0 {
                continue;
            }

            // Calculate means for 4 symbols simultaneously
            let mut sum_x = f64x4::ZERO;
            let mut sum_y = f64x4::ZERO;

            for i in 0..n_samples {
                // Load volume imbalances for 4 symbols (NaN-safe)
                let x_vec = f64x4::from([
                    volume_imbalances[base_idx][i],
                    volume_imbalances[base_idx + 1][i],
                    volume_imbalances[base_idx + 2][i],
                    volume_imbalances[base_idx + 3][i],
                ]);

                // Load price changes for 4 symbols (NaN-safe)
                let y_vec = f64x4::from([
                    price_changes[base_idx][i],
                    price_changes[base_idx + 1][i],
                    price_changes[base_idx + 2][i],
                    price_changes[base_idx + 3][i],
                ]);

                // Accumulate sums (NaN-safe addition)
                sum_x += x_vec;
                sum_y += y_vec;
            }

            // Calculate means
            let n_vec = f64x4::splat(n_samples as f64);
            let mean_x = sum_x / n_vec;
            let mean_y = sum_y / n_vec;

            // Calculate regression coefficients
            let mut numerator = f64x4::ZERO;
            let mut denominator = f64x4::ZERO;

            for i in 0..n_samples {
                // Load data for 4 symbols (NaN-safe)
                let x_vec = f64x4::from([
                    volume_imbalances[base_idx][i],
                    volume_imbalances[base_idx + 1][i],
                    volume_imbalances[base_idx + 2][i],
                    volume_imbalances[base_idx + 3][i],
                ]);

                let y_vec = f64x4::from([
                    price_changes[base_idx][i],
                    price_changes[base_idx + 1][i],
                    price_changes[base_idx + 2][i],
                    price_changes[base_idx + 3][i],
                ]);

                // Calculate differences from means (NaN-safe)
                let x_diff = x_vec - mean_x;
                let y_diff = y_vec - mean_y;

                // Update numerator and denominator (NaN-safe)
                let xy_prod = x_diff * y_diff;
                let xx_prod = x_diff * x_diff;

                numerator += xy_prod;
                denominator += xx_prod;
            }

            // Calculate lambda = numerator / denominator (NaN-safe division)
            let lambda = numerator / denominator;

            // Extract results safely
            let lambda_array = lambda.as_array_ref();
            results[base_idx] = lambda_array[0];
            results[base_idx + 1] = lambda_array[1];
            results[base_idx + 2] = lambda_array[2];
            results[base_idx + 3] = lambda_array[3];

            // Handle NaN results (division by zero or invalid inputs)
            for j in 0..4 {
                if results[base_idx + j].is_nan() || results[base_idx + j].is_infinite() {
                    results[base_idx + j] = 0.0;
                }
            }
        }

        // Handle remaining symbols (scalar fallback)
        for i in (chunks * 4)..n_symbols {
            results[i] = calculate_single_lambda(price_changes[i], volume_imbalances[i]);
        }

        results
    }

    /// Calculate Kyle's Lambda for a single symbol using scalar operations
    ///
    /// This is the fallback function for individual symbols or when SIMD processing
    /// is not beneficial. Handles NaN values gracefully.
    fn calculate_single_lambda(price_changes: &[f64], volume_imbalances: &[f64]) -> f64 {
        if price_changes.is_empty() || volume_imbalances.is_empty() {
            return 0.0;
        }

        if price_changes.len() != volume_imbalances.len() {
            return 0.0;
        }

        let _n = price_changes.len() as f64;

        // Calculate means with NaN handling
        let sum_x: f64 = volume_imbalances.iter().filter(|x| x.is_finite()).sum();
        let sum_y: f64 = price_changes.iter().filter(|y| y.is_finite()).sum();
        let valid_count = price_changes
            .iter()
            .zip(volume_imbalances.iter())
            .filter(|(y, x)| y.is_finite() && x.is_finite())
            .count() as f64;

        if valid_count < 2.0 {
            return 0.0; // Need at least 2 valid samples
        }

        let mean_x = sum_x / valid_count;
        let mean_y = sum_y / valid_count;

        let mut numerator = 0.0;
        let mut denominator = 0.0;

        // Calculate regression coefficients with NaN handling
        for i in 0..price_changes.len() {
            let x = volume_imbalances[i];
            let y = price_changes[i];

            if x.is_finite() && y.is_finite() {
                let x_diff = x - mean_x;
                let y_diff = y - mean_y;
                numerator += x_diff * y_diff;
                denominator += x_diff * x_diff;
            }
        }

        if denominator.abs() < 1e-10 || !denominator.is_finite() {
            0.0
        } else {
            let result = numerator / denominator;
            if result.is_finite() { result } else { 0.0 }
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use rust_decimal::Decimal;
    use rust_decimal_macros::dec;

    #[test]
    fn test_kyles_lambda_basic() {
        let mut calculator = KylesLambdaCalculator::new(20, 100); // 100ms intervals for faster testing

        // Simulate trades with clear price impact pattern
        let base_time = 1_000_000_000u64;
        let mut last_lambda = None;

        // Generate 15 intervals to ensure we have enough samples (need 10 minimum)
        for interval in 0..15 {
            let interval_time = base_time + (interval as u64) * 200_000_000; // 200ms between intervals
            let price_base = dec!(100.0) + Decimal::from(interval) / Decimal::from(10); // Price increases gradually

            // Alternate between net buy and net sell intervals
            let (side, volume) = if interval % 2 == 0 {
                (TradeSide::Buy, dec!(100.0)) // Net buy volume
            } else {
                (TradeSide::Sell, dec!(50.0)) // Net sell volume
            };

            // Create a few trades within each interval
            for j in 0..3 {
                let trade = TradeTick {
                    symbol: "TEST".into(),
                    timestamp_ns: interval_time + j * 20_000_000, // 20ms apart
                    side,
                    price: price_base,
                    quantity: volume / Decimal::from(3),
                };
                let lambda = calculator.update_trade(&trade);
                if lambda.is_some() {
                    last_lambda = lambda;
                }
            }
        }

        // After 15 intervals, we should have a lambda estimate
        assert!(
            last_lambda.is_some(),
            "Expected Kyle's Lambda to be calculated after 15 intervals"
        );
    }

    #[test]
    fn test_market_impact_analyzer() {
        let mut analyzer = MarketImpactAnalyzer::new(20, 1000);

        let trade = TradeTick {
            symbol: "TEST".into(),
            timestamp_ns: 1_000_000_000,
            side: TradeSide::Buy,
            price: dec!(100.0),
            quantity: dec!(50.0),
        };

        let metrics = analyzer.update(&trade);

        // Initially no lambda estimate
        assert!(metrics.kyles_lambda.is_none());
        assert_eq!(metrics.instantaneous_impact, 0.0);
    }

    #[test]
    fn test_simd_batch_calculation() {
        use crate::features::market_impact::simd::calculate_kyles_lambda_batch;

        // Create test data for 4 symbols
        let price_changes = [
            vec![0.1, 0.2, -0.1, 0.3, -0.2],
            vec![0.15, 0.25, -0.05, 0.35, -0.15],
            vec![0.05, 0.15, -0.15, 0.25, -0.25],
            vec![0.2, 0.3, -0.2, 0.4, -0.1],
        ];

        let volume_imbalances = [
            vec![100.0, 200.0, -100.0, 300.0, -200.0],
            vec![150.0, 250.0, -50.0, 350.0, -150.0],
            vec![50.0, 150.0, -150.0, 250.0, -250.0],
            vec![200.0, 300.0, -200.0, 400.0, -100.0],
        ];

        let price_refs: Vec<&[f64]> = price_changes.iter().map(|v| v.as_slice()).collect();
        let volume_refs: Vec<&[f64]> = volume_imbalances.iter().map(|v| v.as_slice()).collect();

        // Now safe - no unsafe block needed!
        let results = calculate_kyles_lambda_batch(&price_refs, &volume_refs);

        assert_eq!(results.len(), 4);
        // All results should be positive given the positive correlation in test data
        for lambda in &results {
            assert!(*lambda > 0.0, "Lambda should be positive, got: {lambda}");
        }
    }

    #[test]
    fn test_simd_nan_safety() {
        use crate::features::market_impact::simd::calculate_kyles_lambda_batch;

        // Test with NaN values to ensure safe handling
        let price_changes = [
            vec![0.1, f64::NAN, -0.1, 0.3, -0.2],
            vec![0.15, 0.25, -0.05, f64::INFINITY, -0.15],
            vec![f64::NEG_INFINITY, 0.15, -0.15, 0.25, -0.25],
            vec![0.2, 0.3, -0.2, 0.4, f64::NAN],
        ];

        let volume_imbalances = [
            vec![100.0, 200.0, f64::NAN, 300.0, -200.0],
            vec![150.0, f64::INFINITY, -50.0, 350.0, -150.0],
            vec![50.0, 150.0, -150.0, f64::NEG_INFINITY, -250.0],
            vec![f64::NAN, 300.0, -200.0, 400.0, -100.0],
        ];

        let price_refs: Vec<&[f64]> = price_changes.iter().map(|v| v.as_slice()).collect();
        let volume_refs: Vec<&[f64]> = volume_imbalances.iter().map(|v| v.as_slice()).collect();

        // Should handle NaN/Infinity gracefully without panicking
        let results = calculate_kyles_lambda_batch(&price_refs, &volume_refs);

        assert_eq!(results.len(), 4);
        // All results should be finite (NaN/Inf converted to 0.0)
        for (i, lambda) in results.iter().enumerate() {
            assert!(
                lambda.is_finite(),
                "Lambda[{i}] should be finite, got: {lambda}"
            );
        }
    }
}