# QUANTCONNECT - WORK IN PROGRESS:
# Strategy Idea - see research jupyter notebook for implementation:
#
# 1. Select a High Market Cap Stock or ETF:
# These assets are liquid and exhibit relatively stable price movements
#
# 2. Analyze Daily Percentage Changes:
# Compute the daily percentage changes from historical closing prices.
# Plot the distribution of these changes, which should approximate a normal distribution 
# for large-cap stocks or ETFs.
#
# 3. Calculate Key Metrics:
# Compute the mean, E[X] and standard deviation σ of the daily percentage changes.
# Example from the simulation on $AAPl: E[X]=0.000877 (about 0.0877%), and σ=0.0196 (about 1.96%).
#
# 4. Set Up Butterfly Spreads:
# Place a 2x butterfly spread at E[X]. 
# Place a 1x butterfly spread at E[X]−σ (left standard deviation).
# Place a 1x butterfly spread at E[X]+σ (right standard deviation).
#
# **My Hypothesis**
# Defined Risk/Reward: Butterfly spreads have limited risk and reward, 
# making them suitable for range-bound strategies.
#
# Utilization of Statistical Metrics: Using historical data to determine 
# expected values and deviations adds a quantitative edge.
#
# High Return Potential: Butterfly spreads can yield high returns if the stock price 
# ends near the strike price at expiration.
# sources: 
# OptionsStrat: 15D till expiration Long Call Butterfly on $SPY: 
# Debit & Max Loss: $77.50 (-100%)
# Max Profit: $422.50 (+544%)
# https://optionstrat.com/build/long-call-butterfly/SPY/.SPY250110C596x2,.SPY250110C601x-4,.SPY250110C606x2
#
# Since butterfly spreads are relatively cheap to open and have a very high max profit, 
# we can open several at our determined values (E[x] & std()'s) where assuming that one
# will profit, it will cover the losses of the other 2
#
# Alternatively, this strategy can be adjusted for iron condors

from AlgorithmImports import *
import numpy as np

class ButterflySpreadStrategy(QCAlgorithm): 
    def initialize(self):
        self.set_start_date(2022, 1, 1)
        self.set_end_date(2022, 10, 1)
        self.set_cash(500000)

        option = self.add_option("SPY", Resolution.Minute)
        self.symbol = option.symbol
        option.set_filter(self.universe_filter)

        self.mean_daily_return = None
        self.std_dev_daily_return = None
        self.positions = {}
        self.position_expiries = {}
        self.position_quantities = {}
        
        self.initialize_statistics()

    def initialize_statistics(self):
        history = self.History(self.symbol.Underlying, 252, Resolution.Daily)
        if not history.empty:
            daily_returns = history['close'].pct_change().dropna()
            self.mean_daily_return = float(np.mean(daily_returns))
            self.std_dev_daily_return = float(np.std(daily_returns))

    def universe_filter(self, universe):
        return (universe
                .include_weeklys()
                .strikes(-5, 5)
                .expiration(timedelta(15), timedelta(45)))

    def calculate_equidistant_strikes(self, atm_price, available_strikes):
        sorted_strikes = sorted(available_strikes)
        atm_index = min(range(len(sorted_strikes)), key=lambda i: abs(sorted_strikes[i] - atm_price))
        
        strike_spacing = 1.0
        for i in range(1, len(sorted_strikes)):
            if atm_index - i >= 0 and atm_index + i < len(sorted_strikes):
                lower_strike = sorted_strikes[atm_index - i]
                middle_strike = sorted_strikes[atm_index]
                upper_strike = sorted_strikes[atm_index + i]
                
                if (abs(middle_strike - lower_strike - strike_spacing) < 0.01 and 
                    abs(upper_strike - middle_strike - strike_spacing) < 0.01):
                    return lower_strike, middle_strike, upper_strike
        
        return None, None, None

    def manage_positions(self, data):
        current_date = self.Time
        positions_to_remove = []
        
        for strategy_id, expiry in self.position_expiries.items():
            days_to_expiry = (expiry - current_date).days
            if days_to_expiry <= 5:
                if strategy_id in self.positions:
                    # Close the position by selling the strategy
                    strategy = self.positions[strategy_id]
                    quantity = self.position_quantities[strategy_id]
                    self.sell(strategy, quantity)
                    positions_to_remove.append(strategy_id)
        
        for strategy_id in positions_to_remove:
            del self.positions[strategy_id]
            del self.position_expiries[strategy_id]
            del self.position_quantities[strategy_id]

    def on_data(self, data):
        self.manage_positions(data)
        
        if self.portfolio.invested or self.mean_daily_return is None:
            return

        chain = data.option_chains.get(self.symbol, None)
        if not chain or not chain.underlying:
            return

        expiry = min([contract.expiry for contract in chain])
        calls = [contract for contract in chain 
                if contract.right == OptionRight.Call 
                and contract.expiry == expiry
                and contract.volume > 10]

        if len(calls) == 0:
            return

        atm_price = float(chain.underlying.price)
        available_strikes = sorted(set([call.strike for call in calls]))

        target_prices = [atm_price]  # Focus on ATM butterflies only

        for target_price in target_prices:
            lower_strike, middle_strike, upper_strike = self.calculate_equidistant_strikes(target_price, available_strikes)
            
            if all([lower_strike, middle_strike, upper_strike]):
                strategy_id = f"butterfly_{lower_strike}_{middle_strike}_{upper_strike}"
                
                if strategy_id not in self.positions:
                    option_strategy = OptionStrategies.call_butterfly(
                        self.symbol, upper_strike, middle_strike, lower_strike, expiry
                    )
                    quantity = 1
                    self.positions[strategy_id] = option_strategy
                    self.position_expiries[strategy_id] = expiry
                    self.position_quantities[strategy_id] = quantity
                    self.buy(option_strategy, quantity)

    def on_order_event(self, order_event):
        if order_event.status == OrderStatus.Filled:
            self.log(f"Order filled - {order_event.symbol}: {order_event.fill_quantity} @ {order_event.fill_price}")


# Since Notebook is not shareable as a public backtest:
# Here is the code to Research:
#
# from QuantConnect import *
# from QuantConnect.Research import *
# import numpy as np
# import pandas as pd
# import matplotlib.pyplot as plt
# import seaborn as sns
# from scipy import stats

# # Initialize QuantBook
# qb = QuantBook()

# # Add SPY data
# spy = qb.AddEquity("SPY")
# symbol = spy.Symbol

# # Get historical data (3 years of daily data)
# history = qb.History(symbol, 252*3, Resolution.Daily)
# if history.empty:
#     raise ValueError("No historical data retrieved")

# # Calculate daily returns
# daily_returns = history['close'].pct_change().dropna()

# # Calculate mean and standard deviation
# mean_return = float(np.mean(daily_returns))
# std_dev = float(np.std(daily_returns))

# # Create the plots
# plt.figure(figsize=(15, 10))

# # Plot 1: Historical Prices
# plt.subplot(2, 2, 1)
# history['close'].plot(title='SPY Historical Prices')
# plt.grid(True)

# # Plot 2: Daily Returns Distribution
# plt.subplot(2, 2, 2)
# sns.histplot(daily_returns, stat='density', kde=True)
# plt.axvline(mean_return, color='r', linestyle='--', label='Mean')
# plt.axvline(mean_return + std_dev, color='g', linestyle='--', label='+1 Std Dev')
# plt.axvline(mean_return - std_dev, color='g', linestyle='--', label='-1 Std Dev')
# plt.title('Distribution of Daily Returns')
# plt.legend()
# plt.grid(True)

# # Plot 3: QQ Plot to verify normality
# plt.subplot(2, 2, 3)
# stats.probplot(daily_returns, dist="norm", plot=plt)
# plt.title('Q-Q Plot of Daily Returns')

# # Calculate butterfly spread strike prices
# current_price = float(history['close'].iloc[-1])
# expected_price = current_price * (1 + mean_return)
# left_wing = current_price * (1 + mean_return - std_dev)
# right_wing = current_price * (1 + mean_return + std_dev)

# # Print strategy details
# print("\nButterfly Spread Strategy Parameters:")
# print(f"Current Price: ${current_price:.2f}")
# print(f"Mean Daily Return: {mean_return:.4%}")
# print(f"Standard Deviation: {std_dev:.4%}")
# print("\nButterfly Spread Strike Prices:")
# print(f"Left Wing (1x): ${left_wing:.2f}")
# print(f"Body (2x): ${expected_price:.2f}")
# print(f"Right Wing (1x): ${right_wing:.2f}")

# plt.tight_layout()
# plt.show()

# # Optional: Get option chain data for verification
# option_chain = qb.GetOptionHistory(symbol, datetime.now())
# strikes = option_chain.GetStrikes()
# if len(strikes) > 0:
#     print("\nAvailable Option Strikes near targets:")
#     strikes = sorted(strikes)
#     print(f"Left Wing Target: {min(strikes, key=lambda x: abs(x - left_wing))}")
#     print(f"Body Target: {min(strikes, key=lambda x: abs(x - expected_price))}")
#     print(f"Right Wing Target: {min(strikes, key=lambda x: abs(x - right_wing))}")






# And here is the possible improvement idea using monte carlo to get a ideal sortino ratio:


# Add after the existing imports
# from datetime import datetime

# # Add these functions after the imports and before the main code
# def simulate_paths(current_price, mean_return, std_dev, days, num_simulations=10000):
#     daily_returns = np.random.normal(mean_return, std_dev, (num_simulations, days))
#     paths = current_price * np.exp(np.cumsum(daily_returns, axis=1))
#     return paths

# def calculate_sortino_ratio(returns, risk_free_rate=0.02):
#     excess_returns = returns - risk_free_rate
#     downside_returns = np.where(returns < 0, returns, 0)
#     downside_std = np.std(downside_returns)
#     if downside_std == 0:
#         return 0
#     return (np.mean(excess_returns)) / downside_std

# def calculate_butterfly_payoff(paths, lower, center, upper):
#     payoffs = np.maximum(paths[:, -1] - lower, 0) - 2 * np.maximum(paths[:, -1] - center, 0) + np.maximum(paths[:, -1] - upper, 0)
#     return payoffs

# def optimize_butterfly_strikes(paths, current_price, mean_return, std_dev):
#     strike_ranges = {
#         'center': np.arange(0.98, 1.02, 0.001) * current_price,
#         'width': np.arange(0.01, 0.05, 0.001) * current_price
#     }
    
#     best_sortino = -np.inf
#     optimal_strikes = None
    
#     for center in strike_ranges['center']:
#         for width in strike_ranges['width']:
#             lower = center - width
#             upper = center + width
            
#             payoffs = calculate_butterfly_payoff(paths, lower, center, upper)
#             returns = payoffs / (2 * (center - lower))  # Approximate cost of butterfly
            
#             sortino = calculate_sortino_ratio(returns)
#             if sortino > best_sortino:
#                 best_sortino = sortino
#                 optimal_strikes = (lower, center, upper)
                
#     return optimal_strikes, best_sortino

# # Add this code after your existing plots and before the option chain verification
# print("\nRunning Monte Carlo Simulation for Optimal Strikes...")

# # Simulate price paths for 30 days
# paths = simulate_paths(current_price, mean_return, std_dev, 30)

# # Find optimal strikes using Monte Carlo
# optimal_strikes, best_sortino = optimize_butterfly_strikes(paths, current_price, mean_return, std_dev)

# print("\nMonte Carlo Optimized Butterfly Spread:")
# print(f"Lower Strike: ${optimal_strikes[0]:.2f}")
# print(f"Center Strike: ${optimal_strikes[1]:.2f}")
# print(f"Upper Strike: ${optimal_strikes[2]:.2f}")
# print(f"Sortino Ratio: {best_sortino:.4f}")

# # Add an additional plot for Monte Carlo paths
# plt.figure(figsize=(10, 6))
# plt.plot(paths[:100].T, alpha=0.1, color='blue')
# plt.axhline(y=current_price, color='r', linestyle='--')
# plt.title('Monte Carlo Price Paths (100 simulations)')
# plt.grid(True)
# plt.show()