feat(timeseries): add time-series integration helpers for financial analysis

- Implement 6 helper functions in src/timeseries_utils.rs:
  * prepare_for_hmm: Feature engineering for HMM regime detection
  * rolling_hurst_exponent: Mean-reversion detection (H < 0.5 = mean-reverting)
  * rolling_half_life: Mean-reversion speed for pairs trading
  * return_statistics: Risk metrics (mean, std, skew, kurt, sharpe)
  * create_lagged_features: ML feature matrix creation
  * rolling_correlation: Rolling correlation for pairs trading

- Add PyO3 bindings in src/timeseries_utils/python_bindings.rs:
  * All functions exposed with _py suffix
  * Proper signature decorators and error handling
  * Registered in lib.rs module system

- Update Python module exports:
  * python/optimizr/core.py: Import from _core
  * python/optimizr/__init__.py: Re-export all functions

- Create comprehensive example:
  * examples/timeseries_integration.py demonstrates all 6 functions
  * Includes integrated pairs trading workflow
  * Shows feature engineering for regime detection

- Technical details:
  * Fixed Array1<f64> type conversions for ndarray compatibility
  * Uses risk_metrics::hurst_exponent and estimate_half_life
  * Built successfully with maturin develop --release (40.93s)
  * All functions tested and working correctly

Part of Priority 3: Time-series integration helpers (Enhancement Strategy)
Addresses v0.3.0 roadmap: Bridge optimization with time-series analysis

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)

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<F: RustObjective>(
+    objective: &F,
+    population: &[Vec<f64>]
+) -> Vec<f64> {
+    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<DEResult>
+```
+
+**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<f64>,    // Successful F values
+    history_cr: Vec<f64>,   // 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<Vec<f64>> {
+    // Create features: returns, lagged returns, etc.
+}
+
+/// Rolling window risk metrics
+pub fn rolling_hurst_exponent(
+    returns: &[f64],
+    window_size: usize,
+) -> Vec<f64> {
+    // Compute Hurst exponent in rolling windows
+}
+
+/// Backtest parameter optimization
+pub fn optimize_strategy_params<F>(
+    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.

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)

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",
 ]

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:

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(())
 }

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<Vec<f64>> {
+    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<f64> {
+    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::<Vec<_>>();
+        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<f64> {
+    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::<f64>() / n;
+
+    // Standard deviation
+    let variance = returns.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / 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::<f64>() / n
+    } else {
+        0.0
+    };
+
+    // Excess kurtosis
+    let kurtosis = if std_dev > 1e-10 {
+        returns.iter().map(|r| ((r - mean) / std_dev).powi(4)).sum::<f64>() / 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<Vec<f64>> {
+    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<f64> {
+    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::<f64>() / window_size as f64;
+        let mean_y = y.iter().sum::<f64>() / 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
+        }
+    }
+}

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<f64>, lag_periods: Vec<usize>) -> Vec<Vec<f64>> {
+    prepare_for_hmm(&prices, &lag_periods)
+}
+
+#[pyfunction]
+#[pyo3(signature = (returns, window_size))]
+pub fn rolling_hurst_exponent_py(returns: Vec<f64>, window_size: usize) -> Vec<f64> {
+    rolling_hurst_exponent(&returns, window_size)
+}
+
+#[pyfunction]
+#[pyo3(signature = (prices, window_size))]
+pub fn rolling_half_life_py(prices: Vec<f64>, window_size: usize) -> Vec<f64> {
+    rolling_half_life(&prices, window_size)
+}
+
+#[pyfunction]
+#[pyo3(signature = (returns,))]
+pub fn return_statistics_py(returns: Vec<f64>) -> (f64, f64, f64, f64, f64) {
+    return_statistics(&returns)
+}
+
+#[pyfunction]
+#[pyo3(signature = (series, lags, include_original=true))]
+pub fn create_lagged_features_py(
+    series: Vec<f64>,
+    lags: Vec<usize>,
+    include_original: bool,
+) -> Vec<Vec<f64>> {
+    create_lagged_features(&series, &lags, include_original)
+}
+
+#[pyfunction]
+#[pyo3(signature = (series1, series2, window_size))]
+pub fn rolling_correlation_py(
+    series1: Vec<f64>,
+    series2: Vec<f64>,
+    window_size: usize,
+) -> Vec<f64> {
+    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(())
+}