Unpredicting the Unpredictable A Chaos Theory-Based Trading Strategy in Backtrader

pip install backtrader yfinance pandas numpy matplotlib scipy

import backtrader as bt
import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform # For correlation dimension
from scipy.stats import linregress # For correlation dimension and Hurst exponent

# Set matplotlib style for better visualization
%matplotlib inline
plt.rcParams['figure.figsize'] = (10, 6)

# Download historical data for Ethereum (ETH-USD)
# Remember the instruction: yfinance download with auto_adjust=False and droplevel(axis=1, level=1).
print("Downloading ETH-USD data from 2021-01-01 to 2024-01-01...")
data = yf.download('ETH-USD', '2021-01-01', '2024-01-01', auto_adjust=False)
data.columns = data.columns.droplevel(1) # Drop the second level of multi-index columns
print("Data downloaded successfully.")
print(data.head()) # Display first few rows of the data

# Create a Backtrader data feed from the pandas DataFrame
data_feed = bt.feeds.PandasData(dataname=data)

class ChaosStrategy(bt.Strategy):
    params = (
        ('lookback', 100),            # Window for chaos analysis (needs to be sufficiently large)
        ('embedding_dim', 5),         # Embedding dimension for phase space reconstruction
        ('delay', 3),                 # Time delay for embedding
        ('lyap_threshold', 0.05),     # Lyapunov exponent threshold for "order"
        ('corr_dim_threshold', 2.0),  # Correlation dimension threshold for "order"
        ('trend_period', 30),         # Simple Moving Average period for trend filter
        ('stop_loss_pct', 0.02),      # 2% stop loss
    )
    
    def __init__(self):
        # Store returns history for chaos analysis
        self.returns_history = [] # Use returns, as chaotic dynamics are often more visible here
        
        # Calculate daily percentage returns
        self.returns = bt.indicators.PctChange(self.data.close, period=1)
        
        # Trend filter (Simple Moving Average on close price)
        self.trend_ma = bt.indicators.SMA(self.data.close, period=self.params.trend_period)
        
        # Variables to store chaos theory measures for the current bar
        self.lyapunov_exponent = 0.0
        self.correlation_dimension = 0.0
        self.hurst_exponent = 0.5 # Default for random walk
        
        # Track orders
        self.order = None
        self.stop_order = None

    def time_delay_embedding(self, series, m, tau):
        """
        Create time-delay embedding vectors for phase space reconstruction.
        Args:
            series (np.array): The time series data (e.g., returns).
            m (int): Embedding dimension.
            tau (int): Time delay.
        Returns:
            np.array: Embedded vectors.
        """
        N = len(series)
        # Ensure enough data points to form at least one embedded vector
        if N < m * tau:
            return np.array([])
            
        # Create embedded vectors
        embedded = np.zeros((N - (m-1)*tau, m))
        for i in range(m):
            embedded[:, i] = series[i*tau : N-(m-1-i)*tau]
        return embedded

    def calculate_lyapunov_exponent(self, embedded_vectors):
        """
        Calculates a simplified largest Lyapunov exponent using Wolf's algorithm logic.
        This is a computationally intensive and approximate implementation.
        """
        if len(embedded_vectors) < 20: # Need sufficient points
            return 0.0
            
        try:
            N = len(embedded_vectors)
            lyap_sum = 0
            count = 0
            
            # Iterate through trajectory points to find nearest neighbors and track divergence
            for i in range(N - 5): # Avoid index out of bounds for future_dist
                # Find the nearest neighbor to embedded_vectors[i]
                distances = np.linalg.norm(embedded_vectors[i] - embedded_vectors, axis=1)
                distances[i] = np.inf  # Exclude self from nearest neighbor search
                nearest_idx = np.argmin(distances)
                
                # Ensure nearest neighbor is also within bounds for future steps
                if nearest_idx < N - 5:
                    initial_dist = distances[nearest_idx]
                    if initial_dist > 0: # Avoid division by zero
                        # Track divergence over 'delay' steps (here hardcoded to 3 steps)
                        future_dist = np.linalg.norm(
                            embedded_vectors[i + 3] - embedded_vectors[nearest_idx + 3]
                        )
                        
                        if future_dist > 0: # Avoid log(0)
                            lyap_sum += np.log(future_dist / initial_dist)
                            count += 1
            
            # Average divergence rate over all considered pairs
            # Divided by the number of steps over which divergence was measured (3 here)
            return lyap_sum / (3 * count) if count > 0 else 0.0
            
        except Exception as e:
            # print(f"Error calculating Lyapunov: {e}")
            return 0.0

    def calculate_correlation_dimension(self, embedded_vectors):
        """
        Estimates the correlation dimension using a simplified Grassberger-Procaccia algorithm.
        This is a computationally intensive and approximate implementation.
        """
        if len(embedded_vectors) < 20: # Needs sufficient points
            return 0.0
            
        try:
            # Calculate all pairwise Euclidean distances between embedded vectors
            distances = pdist(embedded_vectors)
            
            if len(distances) == 0 or np.all(distances == 0):
                return 0.0
                
            # Define a range of radii (epsilon) values on a logarithmic scale
            # Ensure min_dist is greater than zero
            min_dist = np.min(distances[distances > 0])
            max_dist = np.max(distances)
            
            if min_dist >= max_dist: # Handle cases where all non-zero distances are the same
                return 0.0
                
            radii = np.logspace(np.log10(min_dist), np.log10(max_dist), 10) # 10 points
            correlation_integrals = []
            
            # Calculate the correlation integral C(r) for each radius r
            # C(r) is the proportion of pairs of points whose distance is less than r
            for r in radii:
                count = np.sum(distances < r)
                correlation_integrals.append(count / len(distances))
            
            # Estimate dimension from the slope of log(C(r)) vs log(r)
            # Take log of radii and integrals, add a small epsilon to integrals to avoid log(0)
            log_radii = np.log(radii[radii > 0]) # Ensure radii are > 0 for log
            log_integrals = np.log(np.array(correlation_integrals)[radii > 0] + 1e-10) # Add epsilon to avoid log(0)
            
            # Perform linear regression on the linear region of the log-log plot
            if len(log_radii) > 3: # Need enough points for regression
                slope, _, _, _, _ = linregress(log_radii, log_integrals)
                # The slope is an estimate of the correlation dimension
                return max(0.0, min(self.params.embedding_dim, slope)) # Bound within reasonable range
            
            return 0.0
            
        except Exception as e:
            # print(f"Error calculating Correlation Dimension: {e}")
            return 0.0

    def calculate_hurst_exponent(self, series):
        """
        Calculate Hurst exponent using the rescaled range (R/S) analysis.
        This is a simplified, single-window R/S calculation.
        """
        if len(series) < 20: # Needs at least 20 points for a reasonable estimate
            return 0.5 # Default to random walk if not enough data
            
        try:
            # 1. Calculate the mean-adjusted series
            mean_series = np.mean(series)
            centered_series = series - mean_series
            
            # 2. Calculate the cumulative sum (deviations from mean)
            cumsum = np.cumsum(centered_series)
            
            # 3. Calculate the range (R)
            R = np.max(cumsum) - np.min(cumsum)
            
            # 4. Calculate the standard deviation (S)
            S = np.std(series)
            
            if S == 0: # Avoid division by zero
                return 0.5
                
            # 5. Calculate the rescaled range (R/S)
            rs_ratio = R / S
            
            if rs_ratio <= 0: # Handle cases where R/S is non-positive (unlikely for real data but for safety)
                return 0.5
                
            # 6. Estimate Hurst exponent: H = log(R/S) / log(N)
            # This is a simplified direct calculation. More robust methods involve plotting log(R/S) vs log(N)
            # for multiple N and taking the slope.
            hurst = np.log(rs_ratio) / np.log(len(series))
            
            return max(0.0, min(1.0, hurst)) # Bound H between 0 and 1
            
        except Exception as e:
            # print(f"Error calculating Hurst: {e}")
            return 0.5

    def analyze_strange_attractor(self, series):
        """
        Combines all chaos metric calculations.
        """
        # Ensure series has enough data for embedding and subsequent calculations
        if len(series) < self.params.lookback:
            return 0.0, 0.0, 0.5 # Return defaults if not enough data
            
        # 1. Phase space reconstruction
        embedded = self.time_delay_embedding(
            series, self.params.embedding_dim, self.params.delay
        )
        
        if len(embedded) < 20: # Need enough embedded points for robust metric calculation
            return 0.0, 0.0, 0.5
            
        # 2. Calculate chaos measures
        lyap = self.calculate_lyapunov_exponent(embedded)
        corr_dim = self.calculate_correlation_dimension(embedded)
        hurst = self.calculate_hurst_exponent(series) # Hurst is usually on original series
        
        return lyap, corr_dim, hurst

    def notify_order(self, order):
        # Standard backtrader notify_order for managing stop-loss
        if order.status in [order.Completed]:
            if order.isbuy() and self.position.size > 0:
                stop_price = order.executed.price * (1 - self.params.stop_loss_pct)
                self.stop_order = self.sell(exectype=bt.Order.Stop, price=stop_price)
                # self.log(f'BUY EXECUTED, Price: {order.executed.price:.2f}, Size: {order.executed.size}, Stop Loss set at: {stop_price:.2f}')
            elif order.issell() and self.position.size < 0:
                stop_price = order.executed.price * (1 + self.params.stop_loss_pct)
                self.stop_order = self.buy(exectype=bt.Order.Stop, price=stop_price)
                # self.log(f'SELL EXECUTED (Short), Price: {order.executed.price:.2f}, Size: {order.executed.size}, Stop Loss set at: {stop_price:.2f}')
        
        if order.status in [order.Completed, order.Canceled, order.Rejected]:
            self.order = None
            if order == self.stop_order:
                self.stop_order = None
                
    def log(self, txt, dt=None):
        ''' Logging function for the strategy '''
        dt = dt or self.datas[0].datetime.date(0)
        # print(f'{dt.isoformat()}, {txt}') # Commented out for cleaner output during backtest
        
    def next(self):
        # Prevent new orders if one is already pending
        if self.order is not None:
            return
            
        # Store current returns in history list
        if not np.isnan(self.returns[0]):
            self.returns_history.append(self.returns[0])
        
        # Keep only the most recent 'lookback' * 2 history (arbitrary, just to ensure enough data for various calculations)
        # The actual `lookback` from params is used when slicing `recent_returns`.
        if len(self.returns_history) > self.params.lookback * 2:
            self.returns_history = self.returns_history[-self.params.lookback * 2:]
            
        # Skip if not enough data for the lookback period
        if len(self.returns_history) < self.params.lookback:
            return
            
        # Perform chaos analysis on the recent returns history
        recent_returns = np.array(self.returns_history[-self.params.lookback:])
        lyap, corr_dim, hurst = self.analyze_strange_attractor(recent_returns)
        
        # Store previous values and update current values for comparison in next iteration
        prev_lyap = self.lyapunov_exponent
        prev_corr_dim = self.correlation_dimension
        
        self.lyapunov_exponent = lyap
        self.correlation_dimension = corr_dim
        self.hurst_exponent = hurst
        
        # --- Trading logic based on Chaos Theory metrics ---
        # Strategy: Attempt to trade when the system transitions from chaotic to more ordered/predictable behavior.
        # This implies a shift to a regime where patterns might be more exploitable.
        
        # Signal 1: Transition from Chaos to Order (Lyapunov drops below threshold AND Correlation Dimension is low)
        # And align with the market trend using SMA.
        
        # Condition for potential order/predictability
        is_becoming_ordered = (lyap < self.params.lyap_threshold and prev_lyap >= self.params.lyap_threshold)
        is_simple_structure = (corr_dim < self.params.corr_dim_threshold)
        
        # Determine current market trend
        is_uptrend = self.data.close[0] > self.trend_ma[0]
        is_downtrend = self.data.close[0] < self.trend_ma[0]

        if is_becoming_ordered and is_simple_structure:
            if is_uptrend: # Enter long in an uptrend when predictability increases
                if self.position.size < 0:  # Close short first if exists
                    self.order = self.close()
                    if self.stop_order is not None: self.cancel(self.stop_order)
                elif not self.position:  # Go long if not in position
                    self.order = self.buy()
                    self.log(f'BUY SIGNAL (Chaos -> Order, Uptrend), Lyap: {lyap:.3f}, CD: {corr_dim:.3f}, Price: {self.data.close[0]:.2f}')
            elif is_downtrend: # Enter short in a downtrend when predictability increases
                if self.position.size > 0:  # Close long first if exists
                    self.order = self.close()
                    if self.stop_order is not None: self.cancel(self.stop_order)
                elif not self.position:  # Go short if not in position
                    self.order = self.sell()
                    self.log(f'SELL SIGNAL (Chaos -> Order, Downtrend), Lyap: {lyap:.3f}, CD: {corr_dim:.3f}, Price: {self.data.close[0]:.2f}')
        
        # Signal 2: Exit when system becomes chaotic again (predictability decreases)
        # This acts as a protective mechanism, exiting when the market becomes too unpredictable.
        elif lyap > self.params.lyap_threshold * 2 and self.position: # Use a higher threshold for exiting
            self.log(f'CLOSING POSITION (Chaos Resuming), Lyap: {lyap:.3f}, CD: {corr_dim:.3f}, Price: {self.data.close[0]:.2f}')
            if self.stop_order is not None:
                self.cancel(self.stop_order)
            self.order = self.close()
            
        # Signal 3 (Alternative/Complementary): Hurst Exponent-based signals
        # Trade when Hurst suggests clear trending (H > 0.6) or mean-reverting (H < 0.4) behavior.
        # This part could be used in conjunction with or as an alternative to Lyap/CorrDim.
        # For this tutorial, we will make it distinct to demonstrate its logic.
        # Note: If these signals are activated, they might override or conflict with Lyap/CorrDim signals.
        # You would typically choose one main set of entry criteria.
        
        # Uncomment and adjust if you want to also trade based on Hurst:
        # elif self.position.size == 0: # Only enter new positions if flat
        #     if (hurst > 0.6 and is_uptrend): # Strong persistence in uptrend
        #         self.order = self.buy()
        #         self.log(f'BUY SIGNAL (Hurst Trending), H: {hurst:.3f}, Price: {self.data.close[0]:.2f}')
        #     elif (hurst < 0.4 and is_downtrend): # Strong anti-persistence in downtrend (might imply reversal to trend)
        #         # This interpretation can be tricky; sometimes anti-persistence is traded as mean-reversion,
        #         # but here, in downtrend, could be interpreted as weakness for further shorting.
        #         # A better interpretation for Hurst < 0.5 + downtrend is a mean-reversion entry if price
        #         # temporarily bounces up, or a strong trending signal if it's "anti-persisting" a rally.
        #         # For simplicity, let's just use > 0.6 for long in uptrend.
        #         # If Hurst < 0.4 usually suggests mean-reversion, so it might signal a BUY in downtrend or SELL in uptrend.
        #         pass # Currently no explicit Hurst short signal

# Create a Cerebro entity
cerebro = bt.Cerebro()

# Add the strategy
cerebro.addstrategy(ChaosStrategy)

# Add the data feed
cerebro.adddata(data_feed)

# Set the sizer: invest 95% of available cash on each trade
cerebro.addsizer(bt.sizers.PercentSizer, percents=95)

# Set starting cash
cerebro.broker.setcash(100000.0) # Start with $100,000

# Set commission (e.g., 0.1% per transaction)
cerebro.broker.setcommission(commission=0.001)

# --- Add Analyzers for comprehensive performance evaluation ---
cerebro.addanalyzer(bt.analyzers.SharpeRatio, _name='sharpe')
cerebro.addanalyzer(bt.analyzers.DrawDown, _name='drawdown')
cerebro.addanalyzer(bt.analyzers.Returns, _name='returns')
cerebro.addanalyzer(bt.analyzers.TradeAnalyzer, _name='tradeanalyzer')
cerebro.addanalyzer(bt.analyzers.SQN, _name='sqn') # System Quality Number

# Print starting portfolio value
print(f'Starting Portfolio Value: ${cerebro.broker.getvalue():,.2f}')

# Run the backtest
print("Running backtest...")
results = cerebro.run()
print("Backtest finished.")

# Print final portfolio value
final_value = cerebro.broker.getvalue()
print(f'Final Portfolio Value: ${final_value:,.2f}')
print(f'Return: {((final_value / 100000) - 1) * 100:.2f}%') # Calculate and print percentage return

# --- Get and print analysis results ---
strat = results[0] # Access the strategy instance from the results

print("\n--- Strategy Performance Metrics ---")

# 1. Returns Analysis
returns_analysis = strat.analyzers.returns.get_analysis()
total_return = returns_analysis.get('rtot', 'N/A') * 100
annual_return = returns_analysis.get('rnorm100', 'N/A')
print(f"Total Return: {total_return:.2f}%")
print(f"Annualized Return: {annual_return:.2f}%")

# 2. Sharpe Ratio (Risk-adjusted return)
sharpe_ratio = strat.analyzers.sharpe.get_analysis()
print(f"Sharpe Ratio: {sharpe_ratio.get('sharperatio', 'N/A'):.2f}")

# 3. Drawdown Analysis (Measure of risk)
drawdown_analysis = strat.analyzers.drawdown.get_analysis()
max_drawdown = drawdown_analysis.get('maxdrawdown', 'N/A')
print(f"Max Drawdown: {max_drawdown:.2f}%")
print(f"Longest Drawdown Duration: {drawdown_analysis.get('maxdrawdownperiod', 'N/A')} bars")

# 4. Trade Analysis (Details about trades)
trade_analysis = strat.analyzers.tradeanalyzer.get_analysis()
total_trades = trade_analysis.get('total', {}).get('total', 0)
won_trades = trade_analysis.get('won', {}).get('total', 0)
lost_trades = trade_analysis.get('lost', {}).get('total', 0)
win_rate = (won_trades / total_trades) * 100 if total_trades > 0 else 0
print(f"Total Trades: {total_trades}")
print(f"Winning Trades: {won_trades} ({win_rate:.2f}%)")
print(f"Losing Trades: {lost_trades} ({100-win_rate:.2f}%)")
print(f"Average Win (PnL): {trade_analysis.get('won',{}).get('pnl',{}).get('average', 'N/A'):.2f}")
print(f"Average Loss (PnL): {trade_analysis.get('lost',{}).get('pnl',{}).get('average', 'N/A'):.2f}")
print(f"Ratio Avg Win/Avg Loss: {abs(trade_analysis.get('won',{}).get('pnl',{}).get('average', 0) / trade_analysis.get('lost',{}).get('pnl',{}).get('average', 1)):.2f}")

# 5. System Quality Number (SQN) - Dr. Van Tharp's measure of system quality
sqn_analysis = strat.analyzers.sqn.get_analysis()
print(f"System Quality Number (SQN): {sqn_analysis.get('sqn', 'N/A'):.2f}")


# --- Plot the results ---
print("\nPlotting results...")
# Adjust matplotlib plotting parameters to prevent warnings with large datasets
plt.rcParams['figure.max_open_warning'] = 0
plt.rcParams['agg.path.chunksize'] = 10000 # Helps with performance for large plots

try:
    # iplot=False for static plot, style='candlestick' for candlestick chart
    # plotreturn=True to show the equity curve in a separate subplot
    # Volume=False to remove the volume subplot as it might not be directly relevant for Chaos visualization
    fig = cerebro.plot(iplot=False, style='candlestick',
                 barup=dict(fill=False, lw=1.0, ls='-', color='green'), # Customize bullish candles
                 bardown=dict(fill=False, lw=1.0, ls='-', color='red'),  # Customize bearish candles
                 plotreturn=True, # Show equity curve
                 numfigs=1, # Ensure only one figure is generated
                 volume=False # Exclude volume plot, as it can clutter for this strategy
                )[0][0] # Access the figure object to save/show
    
    plt.show() # Display the plot
except Exception as e:
    print(f"Plotting error: {e}")
    print("Strategy completed successfully, but plotting was skipped due to an error.")