Introduction: The Millisecond Battleground

High-Frequency Market Making Mechanics:
┌─────────────────────────────────────────────────────────────────────┐
|  Market Maker Quote Cycle (Microseconds):                         |
|                                                                  |
|  1. Read current LOB state (bid/ask volumes, depths).            |
|  2. Calculate optimal bid/ask prices (Avellaneda-Stoikov).       |
|  3. Update resting limit orders (cancel/replace).                |
|  4. Monitor fills and adjust inventory.                          |
|  5. Repeat 5,000+ times per second.                             |
|                                                                  |
|  Profit Model:                                                   |
|  Total Revenue = (Spread/2) × Executed_Volume                   |
|  - Inventory_Holding_Costs                                       |
|  - Adverse_Selection_Losses                                     |
└─────────────────────────────────────────────────────────────────────┘

Latency Arbitrage Example:
┌─────────────────────────────────────────────────────────────────────┐
|  Event: Large trade executes on CME (Chicago).                    |
|                                                                  |
|  Propagation Paths:                                              |
|  CME Chicago ──(4.2ms)──▶ NYSE New Jersey (Primary)             |
|  CME Chicago ──(4.2ms)──▶ NYSE New Jersey (Secondary)           |
|                                                                  |
|  If HFT Firm A is 0.5ms faster than HFT Firm B:                 |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Firm A sends order 0.5ms earlier, capturing the arbitrage │   |
|  │  spread before Firm B's order arrives.                     │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                                                                  |
|  Arbitrage Profit = (Price_CME - Price_NYSE) - Latency_Delta    |
|  The race is decided by physical proximity and hardware speed.  |
└─────────────────────────────────────────────────────────────────────┘

Co-Location Architecture:
┌─────────────────────────────────────────────────────────────────────┐
|  Exchange Data Center (e.g., Equinix NY4)                         |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Exchange Matching Engine (Core)                          │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                              │                                    |
|                    ┌─────────┴─────────┐                          |
|                    │   5-10 microsecond  │                          |
|                    │   fiber connection  │                          |
|                    └─────────┬─────────┘                          |
|                              │                                    |
|  ┌───────────────────────────┴───────────────────────────────┐   |
|  │  HFT Firm Server Rack (Co-located)                       │   |
|  │  ┌─────────────────────────────────────────────────────┐ │   |
|  │  │  FPGA/ASIC Hardware (Market Data → Decision → Order)│ │   |
|  │  └─────────────────────────────────────────────────────┘ │   |
|  └───────────────────────────────────────────────────────────┘   |
|                                                                  |
|  Distance from matching engine: < 100 meters                   |
|  Round-trip latency: < 20 microseconds                          |
└─────────────────────────────────────────────────────────────────────┘

Network Stack Latency Comparison:
┌─────────────────────────────────────────────────────────────────────┐
|  Standard OS Stack:                                               |
|  Network Wire → Kernel NIC Driver → OS Network Stack →          |
|  → Kernel→User Context Switch → Application → ...              |
|  Latency: ~50-100 microseconds                                  |
|                                                                  |
|  Kernel Bypass (DPDK):                                           |
|  Network Wire → NIC → Application Memory (Direct)               |
|  Latency: ~1-5 microseconds                                     |
|                                                                  |
|  FPGA Hardware:                                                  |
|  Network Wire → FPGA Logic (All processing in hardware)         |
|  Latency: < 1 microsecond (nanoseconds)                         |
└─────────────────────────────────────────────────────────────────────┘

OBI_t = (V_t^(bid) - V_t^(ask)) / (V_t^(bid) + V_t^(ask))

Where:
- V_t^(bid) = Total bid volume at the top k price levels.
- V_t^(ask) = Total ask volume at the top k price levels.

OBI Signal Interpretation:
┌─────────────────────────────────────────────────────────────────────┐
|  OBI Value        │  Interpretation                              |
|─────────────────────────────────────────────────────────────────────|
|  +0.8 to +1.0     │  Heavy buying pressure → Upward price tick   |
|                   │  probability > 80%                           |
|─────────────────────────────────────────────────────────────────────|
|  +0.2 to +0.8     │  Moderate buying pressure                    |
|─────────────────────────────────────────────────────────────────────|
|  -0.2 to +0.2     │  Neutral book; balanced pressure             |
|─────────────────────────────────────────────────────────────────────|
|  -0.8 to -0.2     │  Moderate selling pressure                   |
|─────────────────────────────────────────────────────────────────────|
|  -1.0 to -0.8     │  Heavy selling pressure → Downward price tick|
|                   │  probability > 80%                           |
└─────────────────────────────────────────────────────────────────────┘

λ(t) = μ(t) + ∫₀ᵗ Σ α_ij · exp(-β_ij · (t - s)) · dN_j(s)

Where:
- λ(t) = Conditional intensity (arrival rate) at time t.
- μ(t) = Baseline arrival rate (exogenous).
- α_ij = Excitation parameter (how much event type j excites event type i).
- β_ij = Decay parameter (how quickly excitation fades over time).
- N_j(s) = Counting process for event type j up to time s.
- Integral represents self-excitation and cross-excitation feedback loops.

Hawkes Process Visualization:
┌─────────────────────────────────────────────────────────────────────┐
|  Event Arrival Intensity                                          |
|  ▲                                                                 |
|  │        ████████████                                             |
|  │       ██████████████     ████████                              |
|  │      ████████████████   ████████████  ████                     |
|  │     ██████████████████ ████████████████████                    |
|  │    ████████████████████████████████████████████                |
|  └─────────────────────────────────────────────────────▶ Time       |
|                                                                  |
|  • Large event triggers excitation (clustering).                 |
|  • Intensity decays exponentially (β parameter).                 |
|  • Multiple event types interact (α cross-excitation).           |
|                                                                  |
|  Application: Predicts order flow bursts to avoid adverse        |
|  selection from informed order clusters.                       |
└─────────────────────────────────────────────────────────────────────┘

r(s, q, t) = s - q · γ · σ² · (T - t)

Where:
- s = Current mid-price.
- q = Current inventory (positive = long, negative = short).
- γ = Risk aversion coefficient (inventory penalty).
- σ = Asset volatility.
- (T - t) = Remaining trading horizon.

Ask_Price = r + (Spread/2) + q · γ · σ² · (T - t)
Bid_Price  = r - (Spread/2) + q · γ · σ² · (T - t)

Avellaneda-Stoikov Quote Optimization:
┌─────────────────────────────────────────────────────────────────────┐
|  Inventory Position   │  Quote Adjustment                        |
|─────────────────────────────────────────────────────────────────────|
|  Long inventory       │  Lower bids & asks to reduce exposure    |
|  (q > 0)              │  (encourage sellers to hit bid)          |
|─────────────────────────────────────────────────────────────────────|
|  Short inventory      │  Raise bids & asks to reduce exposure    |
|  (q < 0)              │  (encourage buyers to hit ask)           |
|─────────────────────────────────────────────────────────────────────|
|  Neutral inventory    │  Symmetric quotes around mid-price       |
|  (q ≈ 0)              │  (maximize spread capture)              |
|─────────────────────────────────────────────────────────────────────|
|  High volatility      │  Wider spreads to compensate for         |
|  (σ high)             │  increased adverse selection risk        |
|─────────────────────────────────────────────────────────────────────|
|  Low volatility       │  Narrower spreads (competitive quotes)   |
|  (σ low)              │                                           |
└─────────────────────────────────────────────────────────────────────┘

import numpy as np
import pandas as pd
from typing import Dict, List, Tuple
from scipy.optimize import minimize

# -------------------- 1. ORDER BOOK IMBALANCE (OBI) CALCULATOR --------------------
class OrderBookImbalance:
    def __init__(self, depth_levels: int = 5):
        self.depth_levels = depth_levels
    
    def calculate_obi(self, bids: List[Tuple[float, float]], 
                      asks: List[Tuple[float, float]]) -> float:
        """
        Calculate Order Book Imbalance.
        bids: List of (price, volume) sorted descending.
        asks: List of (price, volume) sorted ascending.
        """
        bid_volume = sum([v for _, v in bids[:self.depth_levels]])
        ask_volume = sum([v for _, v in asks[:self.depth_levels]])
        
        if bid_volume + ask_volume == 0:
            return 0.0
        
        obi = (bid_volume - ask_volume) / (bid_volume + ask_volume)
        return obi

# -------------------- 2. HAWKES PROCESS SIMULATOR --------------------
class HawkesProcess:
    def __init__(self, mu: float, alpha: float, beta: float):
        """
        mu: Baseline intensity.
        alpha: Self-excitation parameter (0 < alpha < beta).
        beta: Decay parameter.
        """
        self.mu = mu
        self.alpha = alpha
        self.beta = beta
        self.event_times = []
        self.intensity = mu
    
    def simulate(self, T: float, max_events: int = 1000) -> List[float]:
        """
        Simulate a Hawkes process up to time T using Ogata's thinning algorithm.
        Returns: List of event times.
        """
        t = 0
        events = []
        self.event_times = []
        
        while t < T and len(events) < max_events:
            # Current intensity
            current_intensity = self.get_intensity(t)
            
            # Sample inter-arrival time from exponential distribution
            u = np.random.random()
            dt = -np.log(u) / current_intensity if current_intensity > 0 else T - t
            t += dt
            
            if t > T:
                break
            
            # Thinning step: accept event with probability intensity / max_intensity
            current_intensity = self.get_intensity(t)
            max_intensity = self.mu + self.alpha * len(events)
            
            if current_intensity > 0 and np.random.random() < current_intensity / max_intensity:
                events.append(t)
                self.event_times.append(t)
        
        return events
    
    def get_intensity(self, t: float) -> float:
        """Calculate conditional intensity at time t."""
        intensity = self.mu
        for tau in self.event_times:
            if tau < t:
                intensity += self.alpha * np.exp(-self.beta * (t - tau))
        return intensity

# -------------------- 3. AVELLANEDA-STOIKOV MARKET MAKER --------------------
class AvellanedaStoikovMarketMaker:
    def __init__(self, gamma: float = 0.1, sigma: float = 0.02, 
                 T: float = 3600, spread_mult: float = 1.5):
        """
        gamma: Risk aversion coefficient.
        sigma: Volatility.
        T: Remaining trading horizon (in seconds).
        spread_mult: Multiplier for optimal spread.
        """
        self.gamma = gamma
        self.sigma = sigma
        self.T = T
        self.spread_mult = spread_mult
        self.inventory = 0  # Current position
        self.mid_price = 100.0
    
    def reservation_price(self, t: float) -> float:
        """Calculate reservation price at time t."""
        remaining = max(0, self.T - t)
        return self.mid_price - self.inventory * self.gamma * self.sigma**2 * remaining
    
    def optimal_spread(self, t: float) -> float:
        """Calculate optimal half-spread."""
        remaining = max(0, self.T - t)
        # Avellaneda-Stoikov spread formula (simplified)
        base_spread = self.gamma * self.sigma**2 * remaining
        # Add inventory adjustment
        inventory_adj = abs(self.inventory) * self.gamma * self.sigma * np.sqrt(remaining)
        half_spread = self.spread_mult * (base_spread + inventory_adj)
        return half_spread
    
    def quote(self, t: float) -> Dict[str, float]:
        """Generate bid and ask quotes."""
        r = self.reservation_price(t)
        half_spread = self.optimal_spread(t)
        
        bid_price = r - half_spread
        ask_price = r + half_spread
        
        return {
            'bid': bid_price,
            'ask': ask_price,
            'reservation_price': r,
            'half_spread': half_spread,
            'inventory': self.inventory
        }
    
    def update_inventory(self, trade_side: str, quantity: float):
        """Update inventory after a trade execution."""
        if trade_side == 'buy':
            self.inventory += quantity
        elif trade_side == 'sell':
            self.inventory -= quantity
    
    def pnl(self, trade_price: float, trade_side: str, quantity: float) -> float:
        """Calculate P&L from a trade."""
        if trade_side == 'buy':
            return -trade_price * quantity
        else:
            return trade_price * quantity

# -------------------- 4. SIMULATION --------------------
def run_hft_simulation():
    print("=" * 70)
    print("HIGH-FREQUENCY TRADING SIMULATION")
    print("=" * 70)
    
    # ---------- OBI Calculation ----------
    print("\n[1] ORDER BOOK IMBALANCE (OBI)")
    obi_calc = OrderBookImbalance(depth_levels=5)
    
    # Sample LOB data
    bids = [(100.05, 500), (100.00, 800), (99.95, 1200), (99.90, 900), (99.85, 700)]
    asks = [(100.10, 400), (100.15, 600), (100.20, 1000), (100.25, 800), (100.30, 500)]
    
    obi = obi_calc.calculate_obi(bids, asks)
    print(f"Bid Volume (top 5): {sum([v for _, v in bids[:5]]):.0f}")
    print(f"Ask Volume (top 5): {sum([v for _, v in asks[:5]]):.0f}")
    print(f"OBI: {obi:.4f}")
    print(f"Signal: {'🟢 BUYING PRESSURE' if obi > 0.3 else '🔴 SELLING PRESSURE' if obi < -0.3 else '⚪ NEUTRAL'}")
    
    # ---------- Hawkes Process ----------
    print("\n[2] HAWKES PROCESS ORDER ARRIVAL")
    hawkes = HawkesProcess(mu=0.5, alpha=0.8, beta=2.0)
    events = hawkes.simulate(T=100, max_events=50)
    
    print(f"Generated {len(events)} events in 100 seconds")
    print(f"Event times (first 10): {[f'{t:.3f}' for t in events[:10]]}")
    
    # Calculate clustering effect
    if len(events) > 1:
        interarrivals = np.diff(events)
        print(f"Mean inter-arrival time: {np.mean(interarrivals):.4f}s")
        print(f"Min inter-arrival time: {np.min(interarrivals):.4f}s")
        print(f"Max inter-arrival time: {np.max(interarrivals):.4f}s")
        print("Clustering detected: Yes (bursts of events)" if np.std(interarrivals) > np.mean(interarrivals) else "No clustering")
    
    # ---------- Avellaneda-Stoikov Market Making ----------
    print("\n[3] AVELLANEDA-STOIKOV MARKET MAKER")
    
    mm = AvellanedaStoikovMarketMaker(
        gamma=0.1, 
        sigma=0.02, 
        T=3600, 
        spread_mult=1.5
    )
    
    # Simulate dynamic quoting over 10 minutes
    print("\nDynamic Quotes (every 60 seconds):")
    print("-" * 50)
    print("Time  │  Mid  │  Bid  │  Ask  │  Inv  │  Reservation")
    print("-" * 50)
    
    for t in range(0, 601, 60):
        # Random walk for mid-price
        mm.mid_price += np.random.randn() * 0.05
        # Simulate random trades affecting inventory
        if np.random.random() < 0.3:
            side = np.random.choice(['buy', 'sell'])
            qty = np.random.uniform(1, 5)
            mm.update_inventory(side, qty)
        
        quotes = mm.quote(t)
        print(f"{t:4d}s │ {mm.mid_price:6.2f} │ {quotes['bid']:6.2f} │ {quotes['ask']:6.2f} │ "
              f"{quotes['inventory']:+6.2f} │ {quotes['reservation_price']:6.2f}")
    
    # Calculate P&L from simulated trades
    print("\nSimulated Trading P&L:")
    trades = [
        ('buy', 100.20, 10),
        ('sell', 100.30, 5),
        ('buy', 100.15, 8),
        ('sell', 100.25, 12),
        ('buy', 100.10, 6),
    ]
    
    total_pnl = 0
    for side, price, qty in trades:
        pnl = mm.pnl(price, side, qty)
        total_pnl += pnl
        print(f"{side.upper():4s} {qty:3.0f} shares @ ${price:.2f} → P&L: ${pnl:+.2f}")
    
    print(f"\nTotal P&L: ${total_pnl:+.2f}")
    
    # ---------- Market Making Profitability ----------
    print("\n[4] MARKET MAKING PROFITABILITY ANALYSIS")
    mm_capture = AvellanedaStoikovMarketMaker(gamma=0.08, sigma=0.015, T=3600, spread_mult=1.2)
    
    # Simulate 1000 quote cycles
    profits = []
    for i in range(1000):
        mm_capture.mid_price = 100 + np.random.randn() * 0.1
        t = np.random.uniform(0, 3600)
        quotes = mm_capture.quote(t)
        
        # Random execution at bid or ask
        if np.random.random() < 0.1:  # 10% fill rate
            if np.random.random() < 0.5:
                # Sell at bid (we provide liquidity)
                profits.append(quotes['bid'] - mm_capture.mid_price)
            else:
                # Buy at ask (we provide liquidity)
                profits.append(mm_capture.mid_price - quotes['ask'])
    
    if profits:
        print(f"Total trades: {len(profits)}")
        print(f"Average profit per trade: ${np.mean(profits):.4f}")
        print(f"Total profit: ${np.sum(profits):.2f}")
        print(f"Win rate: {(np.array(profits) > 0).mean()*100:.1f}%")

# -------------------- 5. EXECUTION --------------------
if __name__ == "__main__":
    run_hft_simulation()