diff --git a/Dockerfile b/Dockerfile index 2165c38..ea902f5 100644 --- a/Dockerfile +++ b/Dockerfile @@ -26,6 +26,11 @@ RUN apt-get update && apt-get install -y \ build-essential \ curl \ git \ + pkg-config \ + libssl-dev \ + libopenblas-dev \ + gfortran \ + patchelf \ && rm -rf /var/lib/apt/lists/* # Install Rust (needed for maturin) diff --git a/ENHANCEMENT_STRATEGY.md b/ENHANCEMENT_STRATEGY.md new file mode 100644 index 0000000..f7ada7c --- /dev/null +++ b/ENHANCEMENT_STRATEGY.md @@ -0,0 +1,407 @@ +# OptimizR Enhancement Strategy + +**Date**: January 2, 2025 +**Context**: Post-Polaroid Phase 4, exploring integration and improvements +**Based On**: v0.2.0 codebase review, roadmap analysis, synergy opportunities + +## Current State Analysis + +### ✅ What's Implemented (v0.2.0) + +1. **Core Algorithms**: + - Differential Evolution (5 strategies: rand1, best1, currenttobest1, rand2, best2) + - Hidden Markov Models (Baum-Welch, Viterbi) + - MCMC Sampling (Metropolis-Hastings, adaptive proposals) + - Grid Search + - Information Theory (mutual information, Shannon entropy) + +2. **Advanced Features (v0.2.0)**: + - Sparse Optimization (Sparse PCA, Box-Tao, Elastic Net) + - Optimal Control (HJB solver, regime switching, jump diffusion) + - Risk Metrics (Hurst exponent, half-life, bootstrap) + - Mathematical Toolkit (numerical differentiation, linear algebra, statistics) + +3. **Architecture**: + - Trait-based design (Optimizer, Sampler, InformationMeasure) + - Builder pattern for configuration + - Functional programming utilities (composition, memoization, pipes) + - Rayon dependency already present + - Feature flag infrastructure (`parallel` feature exists) + +### ⚠️ What's Missing/Incomplete + +1. **Parallelization BLOCKED**: + - Infrastructure exists (Rayon trait, ParallelExecutor trait in core.rs) + - DE has `parallel` parameter but **disabled** due to Python GIL + - Comment: "Python callbacks cannot be safely parallelized due to GIL" + - Grid search marked as "future: Expected 50-100x speedup" + +2. **Advanced DE Variants (Roadmap v0.3.0)**: + - JADE (jDE with archive) + - SHADE (Success-History based Adaptive DE) + - L-SHADE (with linear population reduction) + - Current: Only basic jDE adaptive control + +3. **Multi-Objective Optimization (Roadmap)**: + - NSGA-DE (Non-dominated Sorting) + - MODE (Multi-Objective DE) + - Pareto front computation + +4. **GPU Acceleration (Roadmap)**: + - CUDA kernels + - OpenCL support + - 10-100× additional speedup + +5. **Additional Algorithms (Roadmap)**: + - Particle Swarm Optimization (PSO) + - CMA-ES (Covariance Matrix Adaptation) + - Simulated Annealing + - Ant Colony Optimization + +## Synergy Opportunities: Polaroid + OptimizR + +### 1. Time-Series Feature Engineering for HMM +**Description**: Use Polaroid's time-series operations to create features for regime detection + +**Implementation**: +```python +# Polaroid: Fast feature creation +df = client.lag(['price'], periods=1) # Lagged prices +df = client.pct_change(['price'], periods=1) # Returns +df = client.diff(['price'], periods=1) # Price changes + +# OptimizR: Regime detection on features +returns = df['price_pct_change'].to_numpy() +hmm = HMM(n_states=3) # Bull, Bear, Sideways +hmm.fit(returns, n_iterations=100) +states = hmm.predict(returns) +``` + +**Value**: +- Polaroid provides fast feature engineering (50-200× faster for large datasets) +- OptimizR provides statistical inference (HMM regime detection) +- Combined: Real-time regime switching for trading strategies + +### 2. Risk Metrics on Time-Series Data +**Description**: Calculate advanced risk metrics using both systems + +**Implementation**: +```python +# Polaroid: Efficient return calculation +df = client.pct_change(['price'], periods=1) +returns = df['price_pct_change'].to_numpy() + +# OptimizR: Risk analysis +hurst = compute_hurst_exponent(returns) # Mean-reversion detection +half_life = estimate_half_life(returns) # Reversion time +risk_metrics = compute_risk_metrics(returns) # Comprehensive suite +``` + +**Value**: +- Fast preprocessing (Polaroid) + sophisticated analysis (OptimizR) +- Useful for pairs trading, mean-reversion strategies +- Real-time risk monitoring + +### 3. Optimal Control with Market Data +**Description**: Dynamic portfolio rebalancing with regime-dependent strategies + +**Implementation**: +```python +# Polaroid: Multi-asset feature creation +df = client.lag(['spy_price', 'vix'], periods=[1, 5, 20]) +df = client.pct_change(['spy_price'], periods=1) + +# OptimizR: Solve optimal control problem +# State: [price, volatility regime] +# Control: portfolio weights +value_fn = solve_hjb_regime_switching(...) +``` + +**Value**: +- Combines fast data processing with optimal control theory +- Regime-dependent strategies (bull vs bear market) +- Practical for HFT and algorithmic trading + +### 4. Parameter Optimization for Trading Strategies +**Description**: Use DE to optimize strategy parameters on time-series data + +**Implementation**: +```python +# Polaroid: Backtest execution (fast data ops) +def backtest_strategy(params): + df = client.lag(['price'], periods=int(params[0])) + # ... strategy logic ... + return -sharpe_ratio # Minimize negative Sharpe + +# OptimizR: Find optimal parameters +result = differential_evolution( + objective_fn=backtest_strategy, + bounds=[(1, 50), (0.01, 0.5)], # [lag_period, threshold] + maxiter=500, + strategy='rand1' +) +``` + +**Value**: +- Polaroid handles heavy data processing +- OptimizR finds optimal parameters +- 74-88× faster than SciPy DE + +## High-Priority Enhancements + +### Priority 1: Enable Parallelization for Pure-Rust Objectives + +**Problem**: `parallel` parameter exists but disabled due to Python GIL issues + +**Solution**: Create Rust-native objective function trait for GIL-free parallelization + +**Implementation Strategy**: +1. Add `RustObjectiveFn` trait separate from Python callbacks +2. Implement parallel evaluation for Rust-native functions +3. Keep Python callbacks sequential (GIL limitation) +4. Enable parallel grid search (no Python callbacks needed for grid) + +**Code Outline**: +```rust +// In src/core.rs or src/differential_evolution.rs + +/// Rust-native objective function (no Python, no GIL) +pub trait RustObjective: Send + Sync { + fn evaluate(&self, x: &[f64]) -> f64; +} + +/// Parallel evaluation for Rust objectives +#[cfg(feature = "parallel")] +fn evaluate_population_parallel( + objective: &F, + population: &[Vec] +) -> Vec { + use rayon::prelude::*; + population.par_iter() + .map(|individual| objective.evaluate(individual)) + .collect() +} + +// Python binding for benchmarking +#[pyfunction] +fn differential_evolution_rust( + objective_name: &str, // "sphere", "rosenbrock", "rastrigin" + bounds: Vec<(f64, f64)>, + parallel: bool, // Now actually works! + ... +) -> PyResult +``` + +**Benefits**: +- 10-100× speedup for built-in test functions (sphere, Rosenbrock, Rastrigin) +- Useful for benchmarking and testing +- Grid search can be parallelized (no callbacks) +- Foundation for future Rust-only mode + +**Effort**: Medium (1-2 hours) + +### Priority 2: Implement SHADE (Success-History Adaptive DE) + +**Problem**: Current adaptive DE uses basic jDE, SHADE is state-of-the-art + +**Solution**: Implement SHADE algorithm from Tanabe & Fukunaga (2013) + +**Key Features**: +- Historical memory of successful parameters (F, CR) +- Weighted random selection from memory +- Better than jDE on CEC benchmarks + +**Implementation Strategy**: +1. Add `SHADE` variant to `DEStrategy` enum +2. Create success history buffer (circular buffer of size H=10-100) +3. Update memory after each successful mutation +4. Sample (F, CR) from history using Cauchy/Normal distributions + +**Code Outline**: +```rust +pub enum DEAdaptive { + None, + JDE, // Current implementation + SHADE, // New: Success-history based + LSHADE, // Future: With linear population reduction +} + +struct SHADEMemory { + history_f: Vec, // Successful F values + history_cr: Vec, // Successful CR values + index: usize, // Circular buffer index + size: usize, // Memory size H +} + +impl SHADEMemory { + fn sample_f(&self) -> f64 { + // Cauchy distribution centered on random history entry + } + + fn sample_cr(&self) -> f64 { + // Normal distribution centered on random history entry + } + + fn update(&mut self, successful_f: f64, successful_cr: f64) { + // Add to circular buffer + } +} +``` + +**Benefits**: +- State-of-the-art adaptive control +- Better than jDE empirically +- Aligns with roadmap (v0.3.0) +- Minimal API changes + +**Effort**: Medium-High (2-4 hours with testing) + +### Priority 3: Time-Series Integration Helpers + +**Problem**: Using Polaroid + OptimizR requires manual glue code + +**Solution**: Create helper functions for common time-series + optimization patterns + +**Implementation Strategy**: +1. Add `timeseries_utils` module to OptimizR +2. Functions for common workflows +3. Optional Polaroid integration (via feature flag) + +**Code Outline**: +```rust +// New module: src/timeseries_utils.rs + +/// Prepare time-series data for HMM regime detection +pub fn prepare_for_hmm( + prices: &[f64], + lag_periods: &[usize], +) -> Vec> { + // Create features: returns, lagged returns, etc. +} + +/// Rolling window risk metrics +pub fn rolling_hurst_exponent( + returns: &[f64], + window_size: usize, +) -> Vec { + // Compute Hurst exponent in rolling windows +} + +/// Backtest parameter optimization +pub fn optimize_strategy_params( + objective_fn: F, + param_bounds: Vec<(f64, f64)>, + n_trials: usize, +) -> DEResult +where F: Fn(&[f64]) -> f64 +{ + // Wrapper around DE with sensible defaults +} +``` + +**Python Bindings**: +```python +from optimizr import timeseries_utils as tsu + +# Prepare features +features = tsu.prepare_for_hmm(prices, lag_periods=[1, 5, 20]) + +# Rolling risk metrics +rolling_hurst = tsu.rolling_hurst_exponent(returns, window_size=252) + +# Strategy optimization +def my_strategy(params): + # ... backtesting logic ... + return sharpe_ratio + +result = tsu.optimize_strategy_params( + my_strategy, + param_bounds=[(1, 50), (0.01, 0.5)], + n_trials=500 +) +``` + +**Benefits**: +- Reduces boilerplate for common use cases +- Makes integration obvious +- Encourages adoption +- Low effort, high value + +**Effort**: Low-Medium (1-2 hours) + +## Secondary Enhancements (Future Work) + +### 4. Multi-Objective Optimization (NSGA-DE) +- **Roadmap**: v0.3.0 +- **Use Case**: Portfolio optimization (maximize return, minimize risk) +- **Effort**: High (5-8 hours) + +### 5. GPU Acceleration +- **Roadmap**: v0.3.0 +- **Use Case**: Massive population sizes (10K-100K individuals) +- **Effort**: Very High (multi-day project) + +### 6. Additional Algorithms (PSO, CMA-ES, etc.) +- **Roadmap**: v0.3.0 +- **Use Case**: Algorithm portfolio for different problem types +- **Effort**: High per algorithm (3-5 hours each) + +## Recommended Implementation Order + +1. **Session 1 (Current)**: Time-Series Integration Helpers (1-2 hours) + - Low effort, immediate value + - Makes Polaroid + OptimizR integration obvious + - Creates examples for documentation + +2. **Session 2**: Enable Rust-Native Parallelization (1-2 hours) + - Unblocks major performance gain + - Grid search parallelization + - Foundation for future work + +3. **Session 3**: Implement SHADE (2-4 hours) + - State-of-the-art adaptive DE + - Aligns with roadmap + - Publishable improvement + +4. **Future**: Multi-objective, GPU, additional algorithms + - Larger projects + - Requires more research + +## Testing Strategy + +For each enhancement: +1. **Unit tests**: Algorithm correctness (sphere function, Rosenbrock) +2. **Benchmarks**: Performance comparison (before/after) +3. **Integration tests**: Polaroid + OptimizR workflows +4. **Documentation**: Usage examples, API docs + +## Git Commit Strategy (per MANDATORY rules) + +Each enhancement gets: +1. Feature branch: `feature/shade-algorithm` or `feature/rust-parallelization` +2. Implementation commits with tests +3. Benchmark results documented +4. Final commit: `feat(de): implement SHADE adaptive DE variant` +5. Push to origin +6. Log to historia/ + +## Success Metrics + +1. **Performance**: + - Rust parallelization: 10-100× speedup on multi-core + - SHADE: 10-20% better convergence than jDE on benchmarks + - Time-series helpers: Zero overhead (pure convenience) + +2. **Usability**: + - Integration examples in documentation + - Clear API documentation + - Python usage examples + +3. **Completeness**: + - All tests passing + - Benchmarks documented + - Changes committed to git + +--- + +**Next Action**: Implement Priority 3 (Time-Series Integration Helpers) as it's lowest effort with immediate value for demonstrating Polaroid + OptimizR synergy. diff --git a/examples/timeseries_integration.py b/examples/timeseries_integration.py new file mode 100644 index 0000000..8f4caef --- /dev/null +++ b/examples/timeseries_integration.py @@ -0,0 +1,256 @@ +""" +Time-Series Integration Helpers - Example Usage +=============================================== + +Demonstrates the 6 time-series utility functions for financial data analysis: +1. prepare_for_hmm_py: Feature engineering for regime detection +2. rolling_hurst_exponent_py: Mean-reversion detection +3. rolling_half_life_py: Pairs trading metrics +4. return_statistics_py: Risk analysis +5. create_lagged_features_py: ML feature creation +6. rolling_correlation_py: Correlation analysis + +These helpers bridge OptimizR's optimization capabilities with time-series analysis, +particularly useful for regime-switching models and pairs trading strategies. +""" + +import optimizr +import numpy as np + + +def example_prepare_for_hmm(): + """Example: Feature engineering for Hidden Markov Models""" + print("\n=== Example 1: prepare_for_hmm_py ===") + + # Simulate stock prices + prices = [100.0, 101.5, 99.8, 102.3, 103.7, 104.2, 103.1, 105.8, 107.2, 106.5] + + # Create feature matrix with 1 and 2-period lags + features = optimizr.prepare_for_hmm_py(prices, [1, 2]) + + print(f"Input: {len(prices)} price points") + print(f"Output: {len(features)} rows x {len(features[0])} columns") + print("\nFeature columns:") + print(" [0] Simple returns") + print(" [1] Log returns") + print(" [2] Volatility proxy (squared returns)") + print(" [3] Lagged returns (lag=1)") + print(" [4] Lagged returns (lag=2)") + print(f"\nFirst row: {[f'{x:.4f}' for x in features[0]]}") + print("\n💡 Use this with OptimizR's HMM for regime detection!") + + +def example_rolling_hurst(): + """Example: Detecting mean-reversion with Hurst exponent""" + print("\n=== Example 2: rolling_hurst_exponent_py ===") + + # Generate mean-reverting returns + np.random.seed(42) + returns = list(np.random.randn(20) * 0.02) + + # Compute rolling Hurst exponent + window = 10 + hurst_values = optimizr.rolling_hurst_exponent_py(returns, window) + + print(f"Returns: {len(returns)} observations") + print(f"Rolling Hurst (window={window}): {len(hurst_values)} values") + print(f"\nHurst values: {[f'{h:.3f}' for h in hurst_values[:5]]}...") + print("\nInterpretation:") + print(" H < 0.5: Mean-reverting (good for pairs trading)") + print(" H = 0.5: Random walk") + print(" H > 0.5: Trending") + + avg_hurst = np.mean(hurst_values) + if avg_hurst < 0.5: + print(f"\n📊 Average H = {avg_hurst:.3f} → Mean-reverting behavior detected!") + elif avg_hurst > 0.5: + print(f"\n📊 Average H = {avg_hurst:.3f} → Trending behavior detected!") + else: + print(f"\n📊 Average H = {avg_hurst:.3f} → Random walk behavior") + + +def example_rolling_half_life(): + """Example: Mean-reversion speed for pairs trading""" + print("\n=== Example 3: rolling_half_life_py ===") + + # Simulate spread between two cointegrated assets + np.random.seed(42) + spread = list(100 + np.cumsum(np.random.randn(30) * 0.5)) + + # Compute rolling half-life + window = 15 + half_lives = optimizr.rolling_half_life_py(spread, window) + + print(f"Spread: {len(spread)} observations") + print(f"Rolling half-life (window={window}): {len(half_lives)} values") + print(f"\nHalf-life values: {[f'{hl:.2f}' for hl in half_lives[:5]]}...") + print("\nInterpretation:") + print(" Lower half-life → Faster mean reversion") + print(" Higher half-life → Slower mean reversion") + + avg_hl = np.mean(half_lives) + print(f"\n📊 Average half-life: {avg_hl:.2f} periods") + print(f" → Spread reverts to mean in ~{avg_hl:.0f} periods on average") + + +def example_return_statistics(): + """Example: Comprehensive risk metrics""" + print("\n=== Example 4: return_statistics_py ===") + + # Sample returns from a trading strategy + returns = [0.02, -0.01, 0.015, 0.025, -0.005, 0.01, -0.02, 0.03, 0.005, -0.015] + + # Compute statistics + mean, std, skew, kurt, sharpe = optimizr.return_statistics_py(returns) + + print(f"Returns: {len(returns)} observations") + print("\nStatistics:") + print(f" Mean return: {mean:.4f} ({mean*100:.2f}%)") + print(f" Volatility (std): {std:.4f}") + print(f" Skewness: {skew:.4f} {'(left-tailed)' if skew < 0 else '(right-tailed)'}") + print(f" Kurtosis: {kurt:.4f} {'(fat tails)' if kurt > 0 else '(thin tails)'}") + print(f" Sharpe ratio: {sharpe:.4f}") + + print("\n📊 Risk Assessment:") + if sharpe > 2.0: + print(" ✅ Excellent risk-adjusted returns") + elif sharpe > 1.0: + print(" ✓ Good risk-adjusted returns") + else: + print(" ⚠️ Moderate risk-adjusted returns") + + +def example_lagged_features(): + """Example: Create features for ML models""" + print("\n=== Example 5: create_lagged_features_py ===") + + # Time series to predict + returns = [0.01, 0.02, -0.01, 0.015, 0.005, -0.005, 0.025, 0.01, -0.01, 0.02] + + # Create lagged feature matrix + lags = [1, 2, 3] + features = optimizr.create_lagged_features_py(returns, lags, include_original=True) + + print(f"Original series: {len(returns)} observations") + print(f"Lagged features: {len(features)} rows x {len(features[0])} columns") + print("\nFeature columns:") + print(f" [0] Original value (t)") + print(f" [1] Lag-1 (t-1)") + print(f" [2] Lag-2 (t-2)") + print(f" [3] Lag-3 (t-3)") + print(f"\nFirst row: {[f'{x:.4f}' for x in features[0]]}") + print("\n💡 Use this for ML prediction models (LSTM, Random Forest, etc.)") + + +def example_rolling_correlation(): + """Example: Pairs trading correlation analysis""" + print("\n=== Example 6: rolling_correlation_py ===") + + # Two potentially cointegrated assets + np.random.seed(42) + asset1_returns = list(np.random.randn(25) * 0.02) + asset2_returns = list(np.random.randn(25) * 0.02 + np.array(asset1_returns) * 0.6) + + # Compute rolling correlation + window = 10 + correlations = optimizr.rolling_correlation_py(asset1_returns, asset2_returns, window) + + print(f"Asset 1 returns: {len(asset1_returns)} observations") + print(f"Asset 2 returns: {len(asset2_returns)} observations") + print(f"Rolling correlation (window={window}): {len(correlations)} values") + print(f"\nCorrelation values: {[f'{c:.3f}' for c in correlations[:5]]}...") + + avg_corr = np.mean(correlations) + print(f"\n📊 Average correlation: {avg_corr:.3f}") + if avg_corr > 0.7: + print(" → Strong positive correlation (good for pairs trading)") + elif avg_corr > 0.3: + print(" → Moderate correlation") + else: + print(" → Weak correlation (not ideal for pairs trading)") + + +def example_integrated_workflow(): + """Example: Complete pairs trading analysis workflow""" + print("\n" + "="*70) + print("=== Integrated Workflow: Pairs Trading Analysis ===") + print("="*70) + + # Generate synthetic pair of assets + np.random.seed(42) + n = 50 + asset1 = list(100 + np.cumsum(np.random.randn(n) * 0.5)) + asset2 = list(100 + np.cumsum(np.random.randn(n) * 0.5 + + (np.array(asset1) - 100) * 0.6)) + + # Compute spread + spread = [a1 - a2 for a1, a2 in zip(asset1, asset2)] + + # Step 1: Check mean-reversion with Hurst exponent + print("\n1. Mean-reversion check (Hurst exponent):") + hurst_values = optimizr.rolling_hurst_exponent_py(spread, 20) + avg_hurst = np.mean(hurst_values) + print(f" Average Hurst: {avg_hurst:.3f}") + mean_reverting = avg_hurst < 0.5 + print(f" Mean-reverting: {'✅ Yes' if mean_reverting else '❌ No'}") + + # Step 2: Estimate reversion speed + print("\n2. Mean-reversion speed (half-life):") + half_lives = optimizr.rolling_half_life_py(spread, 20) + avg_hl = np.mean(half_lives) + print(f" Average half-life: {avg_hl:.2f} periods") + print(f" Reversion time: ~{avg_hl:.0f} periods") + + # Step 3: Check correlation stability + print("\n3. Correlation stability:") + returns1 = [(asset1[i] - asset1[i-1])/asset1[i-1] for i in range(1, len(asset1))] + returns2 = [(asset2[i] - asset2[i-1])/asset2[i-1] for i in range(1, len(asset2))] + correlations = optimizr.rolling_correlation_py(returns1, returns2, 15) + avg_corr = np.mean(correlations) + print(f" Average correlation: {avg_corr:.3f}") + print(f" Correlation stability: {'✅ High' if avg_corr > 0.7 else '⚠️ Moderate' if avg_corr > 0.3 else '❌ Low'}") + + # Step 4: Risk metrics for spread returns + print("\n4. Spread risk metrics:") + spread_returns = [(spread[i] - spread[i-1]) for i in range(1, len(spread))] + mean, std, skew, kurt, sharpe = optimizr.return_statistics_py(spread_returns) + print(f" Mean: {mean:.4f}, Volatility: {std:.4f}") + print(f" Sharpe: {sharpe:.3f}") + + # Final recommendation + print("\n" + "="*70) + print("📊 Trading Recommendation:") + if mean_reverting and avg_hl < 20 and avg_corr > 0.5: + print("✅ STRONG PAIR: Good candidate for pairs trading") + print(f" - Fast mean reversion ({avg_hl:.1f} periods)") + print(f" - Stable correlation ({avg_corr:.2f})") + print(f" - Predictable behavior (H={avg_hurst:.2f})") + elif mean_reverting and avg_corr > 0.3: + print("⚠️ MODERATE PAIR: Consider with caution") + print(f" - Mean reversion detected") + print(f" - Moderate correlation ({avg_corr:.2f})") + else: + print("❌ WEAK PAIR: Not recommended for pairs trading") + print(f" - Low correlation or trending behavior") + print("="*70) + + +if __name__ == "__main__": + print("=" * 70) + print("OptimizR Time-Series Integration Helpers") + print("=" * 70) + + # Run all examples + example_prepare_for_hmm() + example_rolling_hurst() + example_rolling_half_life() + example_return_statistics() + example_lagged_features() + example_rolling_correlation() + + # Integrated workflow + example_integrated_workflow() + + print("\n" + "=" * 70) + print("✅ All examples completed successfully!") + print("=" * 70) diff --git a/python/optimizr/__init__.py b/python/optimizr/__init__.py index 5109da6..495c66d 100644 --- a/python/optimizr/__init__.py +++ b/python/optimizr/__init__.py @@ -23,6 +23,13 @@ from optimizr.core import ( compute_risk_metrics_py, estimate_half_life_py, bootstrap_returns_py, + # Time-series utilities + prepare_for_hmm_py, + rolling_hurst_exponent_py, + rolling_half_life_py, + return_statistics_py, + create_lagged_features_py, + rolling_correlation_py, ) # Try to import maths_toolkit from Rust backend @@ -47,5 +54,12 @@ __all__ = [ "compute_risk_metrics_py", "estimate_half_life_py", "bootstrap_returns_py", + # Time-series utilities + "prepare_for_hmm_py", + "rolling_hurst_exponent_py", + "rolling_half_life_py", + "return_statistics_py", + "create_lagged_features_py", + "rolling_correlation_py", "maths_toolkit", ] diff --git a/python/optimizr/core.py b/python/optimizr/core.py index 0a3c683..f462f23 100644 --- a/python/optimizr/core.py +++ b/python/optimizr/core.py @@ -21,6 +21,13 @@ try: compute_risk_metrics_py, estimate_half_life_py, bootstrap_returns_py, + # Time-series utilities + prepare_for_hmm_py, + rolling_hurst_exponent_py, + rolling_half_life_py, + return_statistics_py, + create_lagged_features_py, + rolling_correlation_py, ) RUST_AVAILABLE = True except ImportError: diff --git a/src/lib.rs b/src/lib.rs index b1942d5..4da1404 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -33,6 +33,7 @@ use pyo3::types::PyModule; pub mod core; pub mod functional; pub mod maths_toolkit; // Mathematical utilities +pub mod timeseries_utils; // Time-series integration helpers // Modular structure (trait-based, generic) pub mod de; @@ -98,5 +99,8 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(risk_metrics::estimate_half_life_py, m)?)?; m.add_function(wrap_pyfunction!(risk_metrics::bootstrap_returns_py, m)?)?; + // Time-series utility functions + timeseries_utils::python_bindings::register_python_functions(m)?; + Ok(()) } diff --git a/src/timeseries_utils.rs b/src/timeseries_utils.rs new file mode 100644 index 0000000..9012cca --- /dev/null +++ b/src/timeseries_utils.rs @@ -0,0 +1,433 @@ +//! Time-Series Utilities for Finance and Trading +//! +//! Helper functions for common workflows combining time-series preprocessing +//! with optimization and statistical inference. Designed to work seamlessly +//! with Polaroid time-series operations and OptimizR algorithms. +//! +//! # Use Cases +//! +//! - Regime detection with HMM +//! - Rolling risk metrics (Hurst exponent, half-life) +//! - Trading strategy parameter optimization +//! - Feature engineering for financial data +//! +//! # Examples +//! +//! ```rust +//! use optimizr::timeseries_utils::*; +//! +//! // Prepare features for HMM regime detection +//! let prices = vec![100.0, 101.0, 99.5, 102.0]; +//! let features = prepare_for_hmm(&prices, &[1, 5]); +//! +//! // Rolling Hurst exponent +//! let returns = vec![0.01, -0.02, 0.015, 0.005]; +//! let rolling_hurst = rolling_hurst_exponent(&returns, 20); +//! ``` + +use crate::risk_metrics::{hurst_exponent, estimate_half_life}; +use ndarray::Array1; + +#[cfg(feature = "python-bindings")] +pub mod python_bindings; + +/// Prepare time-series price data for HMM regime detection +/// +/// Creates features from price series including: +/// - Returns (percent change) +/// - Lagged returns +/// - Log returns +/// - Volatility proxy (absolute returns) +/// +/// # Arguments +/// +/// * `prices` - Raw price series (e.g., stock prices, crypto prices) +/// * `lag_periods` - Lags to compute for returns (e.g., [1, 5, 20] for daily, weekly, monthly) +/// +/// # Returns +/// +/// Matrix where each row is a feature vector suitable for HMM training. +/// First row will have NaN values due to lagging. +/// +/// # Example +/// +/// ``` +/// use optimizr::timeseries_utils::prepare_for_hmm; +/// +/// let prices = vec![100.0, 101.0, 99.5, 102.0, 103.5]; +/// let features = prepare_for_hmm(&prices, &[1, 2]); +/// +/// // Features: [return_t, return_t-1, return_t-2, log_return, abs_return] +/// assert_eq!(features[0].len(), 5); +/// ``` +pub fn prepare_for_hmm(prices: &[f64], lag_periods: &[usize]) -> Vec> { + let n = prices.len(); + if n < 2 { + return vec![]; + } + + // Calculate returns + let mut returns = Vec::with_capacity(n - 1); + let mut log_returns = Vec::with_capacity(n - 1); + let mut abs_returns = Vec::with_capacity(n - 1); + + for i in 1..n { + let ret = (prices[i] - prices[i - 1]) / prices[i - 1]; + returns.push(ret); + log_returns.push((prices[i] / prices[i - 1]).ln()); + abs_returns.push(ret.abs()); + } + + // Create feature matrix + let max_lag = *lag_periods.iter().max().unwrap_or(&0); + let mut features = Vec::new(); + + for t in max_lag..returns.len() { + let mut feature_vec = Vec::new(); + + // Current return + feature_vec.push(returns[t]); + + // Lagged returns + for &lag in lag_periods { + if t >= lag { + feature_vec.push(returns[t - lag]); + } else { + feature_vec.push(f64::NAN); + } + } + + // Log return and volatility proxy + feature_vec.push(log_returns[t]); + feature_vec.push(abs_returns[t]); + + features.push(feature_vec); + } + + features +} + +/// Compute Hurst exponent in rolling windows +/// +/// The Hurst exponent H indicates: +/// - H < 0.5: Mean-reverting series +/// - H = 0.5: Random walk (no memory) +/// - H > 0.5: Trending series (momentum) +/// +/// # Arguments +/// +/// * `returns` - Return series (percent changes or log returns) +/// * `window_size` - Rolling window size (e.g., 252 for 1 year of daily data) +/// +/// # Returns +/// +/// Vector of Hurst exponents, one per window. Length = returns.len() - window_size + 1 +/// +/// # Example +/// +/// ``` +/// use optimizr::timeseries_utils::rolling_hurst_exponent; +/// +/// let returns = vec![0.01; 300]; // Synthetic data +/// let rolling_h = rolling_hurst_exponent(&returns, 252); +/// +/// assert_eq!(rolling_h.len(), 300 - 252 + 1); +/// ``` +pub fn rolling_hurst_exponent(returns: &[f64], window_size: usize) -> Vec { + let n = returns.len(); + if n < window_size { + return vec![]; + } + + let mut rolling_h = Vec::with_capacity(n - window_size + 1); + + for i in 0..=(n - window_size) { + let window = &returns[i..i + window_size]; + let window_arr = Array1::from_vec(window.to_vec()); + // Use standard window sizes for R/S analysis + let window_sizes = vec![8, 16, 32].into_iter().filter(|&w| w <= window_size / 4).collect::>(); + let h = if window_sizes.is_empty() { + 0.5 // Default to random walk if window too small + } else { + hurst_exponent(&window_arr, &window_sizes) + .map(|result| result.hurst_exponent) + .unwrap_or(0.5) + }; + rolling_h.push(h); + } + + rolling_h +} + +/// Compute half-life of mean reversion in rolling windows +/// +/// Half-life is the expected time for a price series to revert halfway +/// to its mean. Useful for pairs trading and mean-reversion strategies. +/// +/// # Arguments +/// +/// * `prices` - Price series (not returns) +/// * `window_size` - Rolling window size +/// +/// # Returns +/// +/// Vector of half-life estimates in same time units as data +/// (e.g., days if daily prices) +/// +/// # Example +/// +/// ``` +/// use optimizr::timeseries_utils::rolling_half_life; +/// +/// let prices = vec![100.0; 300]; // Synthetic data +/// let rolling_hl = rolling_half_life(&prices, 100); +/// +/// assert_eq!(rolling_hl.len(), 300 - 100 + 1); +/// ``` +pub fn rolling_half_life(prices: &[f64], window_size: usize) -> Vec { + let n = prices.len(); + if n < window_size { + return vec![]; + } + + let mut rolling_hl = Vec::with_capacity(n - window_size + 1); + + for i in 0..=(n - window_size) { + let window = &prices[i..i + window_size]; + let window_arr = Array1::from_vec(window.to_vec()); + let hl = estimate_half_life(&window_arr).unwrap_or(f64::INFINITY); + rolling_hl.push(hl); + } + + rolling_hl +} + +/// Calculate basic time-series statistics +/// +/// Returns summary statistics useful for risk analysis and diagnostics. +/// +/// # Arguments +/// +/// * `returns` - Return series +/// +/// # Returns +/// +/// Tuple of (mean, std_dev, skewness, kurtosis, sharpe_ratio) +/// +/// # Example +/// +/// ``` +/// use optimizr::timeseries_utils::return_statistics; +/// +/// let returns = vec![0.01, -0.02, 0.015, 0.005, -0.01]; +/// let (mean, std, skew, kurt, sharpe) = return_statistics(&returns); +/// ``` +pub fn return_statistics(returns: &[f64]) -> (f64, f64, f64, f64, f64) { + let n = returns.len() as f64; + if returns.is_empty() { + return (0.0, 0.0, 0.0, 0.0, 0.0); + } + + // Mean + let mean = returns.iter().sum::() / n; + + // Standard deviation + let variance = returns.iter().map(|r| (r - mean).powi(2)).sum::() / n; + let std_dev = variance.sqrt(); + + // Skewness + let skewness = if std_dev > 1e-10 { + returns.iter().map(|r| ((r - mean) / std_dev).powi(3)).sum::() / n + } else { + 0.0 + }; + + // Excess kurtosis + let kurtosis = if std_dev > 1e-10 { + returns.iter().map(|r| ((r - mean) / std_dev).powi(4)).sum::() / n - 3.0 + } else { + 0.0 + }; + + // Sharpe ratio (assuming 252 trading days, risk-free rate = 0) + let sharpe_ratio = if std_dev > 1e-10 { + (mean * 252.0_f64.sqrt()) / std_dev + } else { + 0.0 + }; + + (mean, std_dev, skewness, kurtosis, sharpe_ratio) +} + +/// Create lagged features for machine learning +/// +/// Useful for creating supervised learning datasets from time series. +/// +/// # Arguments +/// +/// * `series` - Input time series (prices, returns, etc.) +/// * `lags` - Vector of lag periods (e.g., [1, 2, 5, 10]) +/// * `include_original` - Whether to include t=0 (current value) +/// +/// # Returns +/// +/// Matrix where each row is [t, t-lag1, t-lag2, ...] if include_original=true +/// or [t-lag1, t-lag2, ...] if false +/// +/// # Example +/// +/// ``` +/// use optimizr::timeseries_utils::create_lagged_features; +/// +/// let prices = vec![100.0, 101.0, 99.5, 102.0, 103.5, 104.0]; +/// let features = create_lagged_features(&prices, &[1, 2], true); +/// +/// // Each row: [price_t, price_t-1, price_t-2] +/// ``` +pub fn create_lagged_features(series: &[f64], lags: &[usize], include_original: bool) -> Vec> { + let n = series.len(); + let max_lag = *lags.iter().max().unwrap_or(&0); + + if n <= max_lag { + return vec![]; + } + + let mut features = Vec::new(); + + for t in max_lag..n { + let mut feature_vec = Vec::new(); + + if include_original { + feature_vec.push(series[t]); + } + + for &lag in lags { + feature_vec.push(series[t - lag]); + } + + features.push(feature_vec); + } + + features +} + +/// Compute rolling correlation between two series +/// +/// Useful for pairs trading and correlation analysis. +/// +/// # Arguments +/// +/// * `series1` - First time series +/// * `series2` - Second time series (must be same length as series1) +/// * `window_size` - Rolling window size +/// +/// # Returns +/// +/// Vector of correlation coefficients +/// +/// # Example +/// +/// ``` +/// use optimizr::timeseries_utils::rolling_correlation; +/// +/// let spy = vec![100.0, 101.0, 99.5, 102.0, 103.5]; +/// let qqq = vec![200.0, 202.0, 199.0, 204.0, 207.0]; +/// let corr = rolling_correlation(&spy, &qqq, 3); +/// ``` +pub fn rolling_correlation(series1: &[f64], series2: &[f64], window_size: usize) -> Vec { + let n = series1.len(); + if n != series2.len() || n < window_size { + return vec![]; + } + + let mut rolling_corr = Vec::with_capacity(n - window_size + 1); + + for i in 0..=(n - window_size) { + let x = &series1[i..i + window_size]; + let y = &series2[i..i + window_size]; + + // Compute correlation + let mean_x = x.iter().sum::() / window_size as f64; + let mean_y = y.iter().sum::() / window_size as f64; + + let mut cov = 0.0; + let mut var_x = 0.0; + let mut var_y = 0.0; + + for j in 0..window_size { + let dx = x[j] - mean_x; + let dy = y[j] - mean_y; + cov += dx * dy; + var_x += dx * dx; + var_y += dy * dy; + } + + let corr = if var_x > 1e-10 && var_y > 1e-10 { + cov / (var_x.sqrt() * var_y.sqrt()) + } else { + 0.0 + }; + + rolling_corr.push(corr); + } + + rolling_corr +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_prepare_for_hmm() { + let prices = vec![100.0, 101.0, 99.5, 102.0, 103.5, 104.0]; + let features = prepare_for_hmm(&prices, &[1, 2]); + + assert!(!features.is_empty()); + // Each feature vector: [return, lag1, lag2, log_return, abs_return] + assert_eq!(features[0].len(), 5); + } + + #[test] + fn test_rolling_hurst_exponent() { + let returns = vec![0.01; 50]; + let rolling_h = rolling_hurst_exponent(&returns, 20); + + assert_eq!(rolling_h.len(), 50 - 20 + 1); + // Constant series should have H close to 0.5 + for h in rolling_h { + assert!((h - 0.5).abs() < 0.3); // Allow some numerical variation + } + } + + #[test] + fn test_return_statistics() { + let returns = vec![0.01, -0.02, 0.015, 0.005, -0.01]; + let (mean, std, skew, kurt, sharpe) = return_statistics(&returns); + + assert!((mean).abs() < 0.1); // Small mean + assert!(std > 0.0); // Non-zero volatility + // Skewness and kurtosis can vary widely + } + + #[test] + fn test_create_lagged_features() { + let series = vec![1.0, 2.0, 3.0, 4.0, 5.0]; + let features = create_lagged_features(&series, &[1, 2], true); + + assert_eq!(features.len(), 3); // 5 - 2 = 3 valid rows + assert_eq!(features[0], vec![3.0, 2.0, 1.0]); // t=2: [val_t, val_t-1, val_t-2] + } + + #[test] + fn test_rolling_correlation() { + let x = vec![1.0, 2.0, 3.0, 4.0, 5.0]; + let y = vec![2.0, 4.0, 6.0, 8.0, 10.0]; // Perfect correlation + let corr = rolling_correlation(&x, &y, 3); + + assert_eq!(corr.len(), 3); + for c in corr { + assert!((c - 1.0).abs() < 1e-6); // Perfect positive correlation + } + } +} diff --git a/src/timeseries_utils/python_bindings.rs b/src/timeseries_utils/python_bindings.rs new file mode 100644 index 0000000..274f0a8 --- /dev/null +++ b/src/timeseries_utils/python_bindings.rs @@ -0,0 +1,62 @@ +//! Python bindings for time-series utilities + +use pyo3::prelude::*; + +use super::{ + create_lagged_features, prepare_for_hmm, return_statistics, rolling_correlation, + rolling_half_life, rolling_hurst_exponent, +}; + +#[pyfunction] +#[pyo3(signature = (prices, lag_periods))] +pub fn prepare_for_hmm_py(prices: Vec, lag_periods: Vec) -> Vec> { + prepare_for_hmm(&prices, &lag_periods) +} + +#[pyfunction] +#[pyo3(signature = (returns, window_size))] +pub fn rolling_hurst_exponent_py(returns: Vec, window_size: usize) -> Vec { + rolling_hurst_exponent(&returns, window_size) +} + +#[pyfunction] +#[pyo3(signature = (prices, window_size))] +pub fn rolling_half_life_py(prices: Vec, window_size: usize) -> Vec { + rolling_half_life(&prices, window_size) +} + +#[pyfunction] +#[pyo3(signature = (returns,))] +pub fn return_statistics_py(returns: Vec) -> (f64, f64, f64, f64, f64) { + return_statistics(&returns) +} + +#[pyfunction] +#[pyo3(signature = (series, lags, include_original=true))] +pub fn create_lagged_features_py( + series: Vec, + lags: Vec, + include_original: bool, +) -> Vec> { + create_lagged_features(&series, &lags, include_original) +} + +#[pyfunction] +#[pyo3(signature = (series1, series2, window_size))] +pub fn rolling_correlation_py( + series1: Vec, + series2: Vec, + window_size: usize, +) -> Vec { + rolling_correlation(&series1, &series2, window_size) +} + +pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(prepare_for_hmm_py, m)?)?; + m.add_function(wrap_pyfunction!(rolling_hurst_exponent_py, m)?)?; + m.add_function(wrap_pyfunction!(rolling_half_life_py, m)?)?; + m.add_function(wrap_pyfunction!(return_statistics_py, m)?)?; + m.add_function(wrap_pyfunction!(create_lagged_features_py, m)?)?; + m.add_function(wrap_pyfunction!(rolling_correlation_py, m)?)?; + Ok(()) +}