Overall Statistics
Total Orders
56
Average Win
0.91%
Average Loss
-1.63%
Compounding Annual Return
-2.098%
Drawdown
11.400%
Expectancy
-0.221
Start Equity
100000
End Equity
89936.16
Net Profit
-10.064%
Sharpe Ratio
-1.481
Sortino Ratio
-0.413
Probabilistic Sharpe Ratio
0.003%
Loss Rate
50%
Win Rate
50%
Profit-Loss Ratio
0.56
Alpha
-0.056
Beta
0.001
Annual Standard Deviation
0.038
Annual Variance
0.001
Information Ratio
-0.773
Tracking Error
0.147
Treynor Ratio
-63.304
Total Fees
$104.11
Estimated Strategy Capacity
$50000000.00
Lowest Capacity Asset
BMY R735QTJ8XC9X
Portfolio Turnover
1.53%
Drawdown Recovery
9
# region imports
from AlgorithmImports import *
from stat_arb import PairStatArb
# endregion


class SpectralStatArb(QCAlgorithm):

    def initialize(self):
        self.set_start_date(self.end_date - timedelta(5 * 365))
        self.set_cash(100_000)
        self.settings.automatic_indicator_warm_up = True
        self.settings.seed_initial_prices = True
        self.universe_settings.resolution = Resolution.DAILY
        self.universe_settings.schedule.on(self.date_rules.month_start("SPY"))
        self._universe = self.add_universe(lambda fundamental: 
            [f.symbol for f in sorted([f for f in fundamental if f.has_fundamental_data and f.market_cap > 5e9],
            key=lambda f: f.dollar_volume)[-100:]]
        )
        self.laplacian_lookback = 252
        self.arb_lookback = 30
        self.entry_z_score = 2.0
        self.exit_z_score = 0.5
        self._model = None
        self._pair_spread = None
        self.set_warm_up(timedelta(45))

    def on_warmup_finished(self):
        # Configure training and prediction schedule.
        trading_time_rule = self.time_rules.at(8, 0)
        training_time_rule = self.time_rules.at(6, 0)
        self.train(self.date_rules.month_end(), training_time_rule, self._train)
        self.schedule.on(self.date_rules.every_day('SPY'), trading_time_rule, self._rebalance)
        # Rebalance today too.
        if self.live_mode:
            self._train()
            self._rebalance()
        else:
            self.train(self.date_rules.today, training_time_rule, self._train)
            self.schedule.on(self.date_rules.today, trading_time_rule, self._rebalance)

    def on_securities_changed(self, changes):
        for security in changes.added_securities:
            # Track returns for clustering.
            security.log_return = self.logr(security, 1)
            security.log_return.window.size = self.laplacian_lookback
            security.log_return.reset()
            history = self.history[TradeBar](security, self.laplacian_lookback + 1, Resolution.DAILY)
            for bar in history:
                security.log_return.update(bar.end_time, bar.close)
        for security in changes.removed_securities:
            # Drop the active pair if one of its legs leaves the universe.
            if self._pair_spread and security.symbol in self._pair_spread.symbols:
                self._drop_pair()
            self.deregister_indicator(security.log_return)

    def _train(self):
        # Initialize stat-arb trainer.
        model = PairStatArb(self.laplacian_lookback, self.arb_lookback)
        # Train pair selection and warm the selected spread indicator.
        pair_spread = model.train(self, self._universe)
        self._drop_pair()
        if pair_spread is None:
            self._model = None
            self.log("No cointegrated pairs found in the cluster this month.")
            return
        # Store model and active spread indicator.
        self._model = model
        self._pair_spread = pair_spread
        self.register_indicator(pair_spread.symbol_a, self._pair_spread, Resolution.DAILY)
        self.register_indicator(pair_spread.symbol_b, self._pair_spread, Resolution.DAILY)
        self.log(
            f"New Cointegrated Pair: {pair_spread.symbol_a.value} & {pair_spread.symbol_b.value} "
            f"(beta={pair_spread.hedge_ratio:.2f})"
        )

    def _rebalance(self):
        pair = self._pair_spread
        if pair is None or not pair.is_ready:
            return
        asset_a, asset_b = pair.symbols
        security_a = self.securities[asset_a]
        security_b = self.securities[asset_b]
        half_life = float(pair.half_life)
        if not np.isfinite(half_life) or half_life > 30:
            return
        dynamic_entry_z = self.entry_z_score if half_life < 10 else self.entry_z_score + 0.5
        z_score = float(pair.z_score)
        if not np.isfinite(z_score):
            return
        invested = security_a.invested or security_b.invested
        if not invested and abs(z_score) > dynamic_entry_z:
            # Short the spread when z-score is high, and long it when z-score is low.
            gross = 1 + abs(pair.hedge_ratio)
            if not np.isfinite(gross) or gross <= 0:
                return
            direction = np.sign(z_score)
            weight_a = -direction / gross
            weight_b = direction * pair.hedge_ratio / gross
            self.set_holdings([PortfolioTarget(asset_a, weight_a), PortfolioTarget(asset_b, weight_b)], False)
            self.log(f"ENTRY: {asset_a.value}/{asset_b.value} | HL: {half_life:.1f}d | Z: {z_score:.2f}")
        elif invested and abs(z_score) < self.exit_z_score:
            # Take profit once the spread reverts toward the mean.
            self.liquidate(asset_a)
            self.liquidate(asset_b)
            self.log("EXIT: Pair converged")

    def _drop_pair(self):
        # Close the pair positions and release the active spread indicator.
        if self._pair_spread:
            for symbol in self._pair_spread.symbols:
                self.liquidate(symbol)
            self.deregister_indicator(self._pair_spread)
        self._pair_spread = None
from AlgorithmImports import *


class PairSpread(PythonIndicator):

    def __init__(self, symbol_a, symbol_b, hedge_ratio: float, period: int, name: str = "PairSpread"):
        super().__init__()
        self.name = name
        self.symbol_a = symbol_a
        self.symbol_b = symbol_b
        self.symbols = (symbol_a, symbol_b)
        self.hedge_ratio = float(hedge_ratio)
        self.value = 0.0
        self.spread = 0.0
        self.z_score = 0.0
        self.half_life = np.inf
        # Store spread values.
        self._window = RollingWindow[float](period)
        # Cache the latest price and timestamp for each leg.
        self._prices = {}
        self._last_time = None

    def update(self, input_: BaseData):
        price = float(input_.value)
        if input_.symbol not in self.symbols or price <= 0:
            return self.is_ready
        self._prices[input_.symbol] = (input_.end_time, price)
        if self.symbol_a not in self._prices or self.symbol_b not in self._prices:
            return self.is_ready
        time_a, price_a = self._prices[self.symbol_a]
        time_b, price_b = self._prices[self.symbol_b]
        if time_a != time_b or time_a == self._last_time:
            return self.is_ready
        self._last_time = time_a
        # Calculate the log-price spread using the selected hedge ratio.
        self.value = float(np.log(price_a) - self.hedge_ratio * np.log(price_b))
        self.spread = self.value
        # Reset signal values until the spread window can recompute them.
        self.z_score = 0.0
        self.half_life = np.inf
        self._window.add(self.value)
        # Recompute z-score and half-life from the complete trailing spread window.
        if self._window.is_ready:
            spreads = np.array(list(self._window)[::-1], dtype=float)
            sd = np.std(spreads)
            if sd > 0 and np.isfinite(sd):
                self.z_score = float((self.value - np.mean(spreads)) / sd)
            # Estimate mean-reversion speed with a lagged spread regression.
            lagged = spreads[:-1]
            if len(lagged) > 1 and not np.all(lagged == lagged[0]):
                slope = np.polyfit(lagged, spreads[1:] - lagged, 1)[0]
                if np.isfinite(slope) and slope < 0:
                    self.half_life = float(-np.log(2) / slope)
        return self.is_ready

    @property
    def is_ready(self) -> bool:
        return self._window.is_ready
from AlgorithmImports import *
from pair_spread import PairSpread
from statsmodels.tsa.stattools import adfuller


class PairStatArb:

    # Initialize the strategy parameters for graph clustering and statistical arbitrage lookbacks.
    def __init__(self, laplacian_lookback: int, arb_lookback: int):
        self._laplacian_lookback = laplacian_lookback
        self._arb_lookback = arb_lookback

    def train(self, algorithm, universe):
        ready_securities = sorted([algorithm.securities[s] for s in universe.selected if algorithm.securities[s].log_return.window.is_ready], key=lambda s: s.symbol.value)
        # Ensure there are at least two valid assets available to form pairs.
        if len(ready_securities) < 2:
            return None
        # Build the log-return matrix from each indicator's own window.
        df_returns = pd.DataFrame({security.symbol.value: [float(point.value) for point in security.log_return.window] for security in ready_securities})
        # Build the weighted correlation graph and its Laplacian.
        weights = np.abs(df_returns.corr().fillna(0).values)
        np.fill_diagonal(weights, 0)
        laplacian = np.diag(np.sum(weights, axis=1)) - weights
        # Compute the eigenvalues and eigenvectors to find the Fiedler vector for spectral clustering.
        _, eigenvectors = np.linalg.eigh(laplacian)
        fiedler_vector = eigenvectors[:, 1]
        # Standardize the Fiedler vector orientation based on its maximum absolute value.
        if fiedler_vector[np.argmax(np.abs(fiedler_vector))] < 0:
            fiedler_vector = -fiedler_vector
        # Select assets belonging to the positive partition of the Fiedler vector.
        cluster_pos = [ready_securities[i].symbol for i in range(len(ready_securities)) if fiedler_vector[i] > 0]
        if len(cluster_pos) < 2:
            return None
        # Fetch the historical close prices for the identified cluster.
        history = algorithm.history(cluster_pos, self._laplacian_lookback, Resolution.DAILY)
        if history.empty:
            return None
        prices = history["close"].unstack(level=0)
        # Iterate over all asset combinations within the cluster to locate a cointegrated pair.
        for i in range(len(cluster_pos)):
            for j in range(i + 1, len(cluster_pos)):
                pair_spread = self._try_create_pair(algorithm, prices, cluster_pos[i], cluster_pos[j])
                if pair_spread:
                    return pair_spread
        return None

    def _try_create_pair(self, algorithm, prices, asset_a, asset_b):
        aligned = pd.concat([np.log(prices[asset_a]), np.log(prices[asset_b])], axis=1).dropna()
        if len(aligned) < self._arb_lookback:
            return None
        y = aligned.iloc[:, 0]
        x = aligned.iloc[:, 1]
        if np.all(x == x.iloc[0]):
            return None
        # Determine the hedge ratio using an ordinary least squares linear regression.
        hedge_ratio = float(np.polyfit(x, y, 1)[0])
        # Construct the spread series using the computed hedge ratio.
        spread = y - hedge_ratio * x
        # Verify the stationarity of the spread using the Augmented Dickey-Fuller test.
        if np.all(spread == spread.iloc[0]) or adfuller(spread)[1] >= 0.05:
            return None
        pair_spread = PairSpread(asset_a, asset_b, hedge_ratio, self._arb_lookback)
        # Fetch historical bar data to warm up the indicator windows of the pair.
        pair_history = algorithm.history[TradeBar]([asset_a, asset_b], self._arb_lookback + 5, Resolution.DAILY)
        for bars in pair_history:
            if asset_a in bars and asset_b in bars:
                pair_spread.update(bars[asset_a])
                pair_spread.update(bars[asset_b])
        return pair_spread