bitcoin-model/model.py

import pandas as pd
import numpy as np
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
from scipy.signal import savgol_filter
from multiprocessing import Process, Pool


# Utility functions


def get_halving_dates():
    """Return known and projected Bitcoin halving dates"""
    return pd.to_datetime(
        [
            "2008-01-03",  # Bitcoin genesis block (treat as cycle start)
            "2012-11-28",  # First halving
            "2016-07-09",  # Second halving
            "2020-05-11",  # Third halving
            "2024-04-19",  # Fourth halving
            "2028-04-20",  # Fifth halving (projected)
        ]
    )


def get_cycle_position(date, halving_dates):
    """
    Calculate position in halving cycle (0 to 1) for a given date.
    0 represents a halving event, 1 represents just before the next halving.
    """
    # Convert date to datetime if it's not already
    date = pd.to_datetime(date)

    # Find the most recent halving before this date
    prev_halving = halving_dates[halving_dates <= date].max()
    if pd.isna(prev_halving):
        return 0.0  # For dates before first halving

    # Find next halving
    future_halvings = halving_dates[halving_dates > date]
    if len(future_halvings) == 0:
        # For dates after last known halving, use same cycle length as last known cycle
        last_cycle_length = (halving_dates[-1] - halving_dates[-2]).days
        days_since_halving = (date - halving_dates[-1]).days
        return min(days_since_halving / last_cycle_length, 1.0)

    next_halving = future_halvings.min()

    # Calculate position as fraction between halvings
    days_since_halving = (date - prev_halving).days
    cycle_length = (next_halving - prev_halving).days
    return min(days_since_halving / cycle_length, 1.0)


def format_price(x, p):
    """Format large numbers in K, M, B format with appropriate precision"""
    if abs(x) >= 1e9:
        return f"${x/1e9:.1f}B"
    if abs(x) >= 1e6:
        return f"${x/1e6:.1f}M"
    if abs(x) >= 1e3:
        return f"${x/1e3:.1f}K"
    if abs(x) >= 1:
        return f"${x:.0f}"
    return f"${x:.2f}"  # For values less than $1, show cents


def get_nice_price_points(min_price, max_price):
    """
    Generate a reasonable set of price points for the y-axis that look clean
    and cover the range without cluttering the chart.
    """
    # Handle zero or negative prices
    min_price = max(min_price, 0.0001)  # Set minimum price to $0.0001

    log_min = np.floor(np.log10(min_price))
    log_max = np.ceil(np.log10(max_price))
    price_points = []

    # For very large ranges (spanning more than 4 orders of magnitude),
    # only use powers of 10 and mid-points
    if log_max - log_min > 4:
        for exp in range(int(log_min), int(log_max + 1)):
            base = 10**exp
            # Add main power of 10
            if min_price <= base <= max_price:
                price_points.append(base)
            # Add mid-point if range is large enough
            if min_price <= base * 5 <= max_price and exp > log_min:
                price_points.append(base * 5)
    else:
        # For smaller ranges, use 1, 2, 5 sequence
        for exp in range(int(log_min), int(log_max + 1)):
            for mult in [1, 2, 5]:
                point = mult * 10**exp
                if min_price <= point <= max_price:
                    price_points.append(point)

    return np.array(price_points)


# Analysis functions


def analyze_trends(df):
    """
    Analyze Bitcoin price trends using log returns with asymmetric dampening.
    """
    df = df.copy()

    # Get halving dates and calculate cycle position
    halving_dates = get_halving_dates()
    df["Cycle_Position"] = df["Date"].apply(
        lambda x: get_cycle_position(x, halving_dates)
    )
    df["Cycle_Days"] = (df["Cycle_Position"] * 4 * 365).round().astype(int)

    # Calculate log returns
    df["Log_Price"] = np.log(df["Close"])
    df["Log_Return"] = df["Log_Price"].diff()

    # Group by position in cycle
    position_returns = df.groupby("Cycle_Days")["Log_Return"].mean()

    # Simple moving average smoothing
    window = 60
    smoothed_returns = position_returns.rolling(
        window=window,
        center=True,
        min_periods=int(window / 2),
    ).mean()

    # Fill any NaN values at the edges
    smoothed_returns = smoothed_returns.fillna(method="bfill").fillna(method="ffill")

    # Apply asymmetric dampening to extreme values
    returns_std = smoothed_returns.std()

    def asymmetric_dampen(x):
        if x > 2 * returns_std:
            return x * 0.6  # Stronger dampening for positive extremes
        elif x < -2 * returns_std:
            return x * 0.7  # Slightly less dampening for negative extremes
        return x

    smoothed_returns = smoothed_returns.map(asymmetric_dampen)

    # Additional dampening based on absolute magnitude
    magnitude_factor = 0.9  # Global dampening factor
    smoothed_returns = smoothed_returns * magnitude_factor

    return smoothed_returns


def calculate_volatility(df, short_window=30, medium_window=90, long_window=180):
    """
    Final volatility calculation with precise period adjustments.
    """
    df = df.copy()
    if "Log_Return" not in df.columns:
        df["Log_Price"] = np.log(df["Close"])
        df["Log_Return"] = df["Log_Price"].diff()

    # Base volatility calculation
    short_vol = df["Log_Return"].ewm(span=short_window, adjust=False).std().iloc[-1]
    medium_vol = df["Log_Return"].ewm(span=medium_window, adjust=False).std().iloc[-1]
    long_vol = df["Log_Return"].ewm(span=long_window, adjust=False).std().iloc[-1]

    # Standard weights
    base_vol = 0.2 * short_vol + 0.5 * medium_vol + 0.3 * long_vol

    # Precise period-specific scaling
    start_date = df["Date"].min()
    if start_date >= pd.Timestamp("2020-01-01"):
        base_adjustment = 0.64  # Slightly increased
    elif start_date >= pd.Timestamp("2016-07-09"):
        base_adjustment = 0.67  # Slightly increased
    elif start_date >= pd.Timestamp("2015-01-01"):
        base_adjustment = 0.69  # Increased for mid period
    else:
        # Early period with less aggressive scaling
        base_adjustment = 0.70  # Fixed value for stability

    return base_vol * base_adjustment


# Era definitions
era_adjustments = {
    "early": {
        "start_date": pd.Timestamp("2013-01-01"),
        "end_date": pd.Timestamp("2017-12-10"),
        "volatility_scale": 0.71,  # Slight increase
        "trend_scale": 0.75,
        "skew_scale": 1.0,
    },
    "transition": {
        "start_date": pd.Timestamp("2017-12-10"),
        "end_date": pd.Timestamp("2020-01-01"),
        "volatility_scale": 0.69,  # Slight increase
        "trend_scale": 0.80,
        "skew_scale": 1.0,
    },
    "mature": {
        "start_date": pd.Timestamp("2020-01-01"),
        "end_date": pd.Timestamp("2100-01-01"),
        "volatility_scale": 0.67,  # Slight increase
        "trend_scale": 0.85,
        "skew_scale": 1.0,
    },
}


def adjust_trend_expectations(expected_returns, cycle_position):
    """
    Simple trend adjustment.
    """
    if cycle_position > 0.75:
        damping_factor = 0.70
    else:
        damping_factor = 0.85

    return expected_returns * damping_factor


def get_projection_adjustments(days_forward, current_cycle_position):
    """
    Final projection adjustments with precise uncertainty scaling.
    """
    adjustments = np.ones(days_forward)

    # Fixed base uncertainty with slight cycle variation
    base_uncertainty = 0.016  # Standard rate
    if current_cycle_position > 0.75:
        base_uncertainty *= 1.1  # 10% increase late cycle

    for i in range(days_forward):
        # Conservative growth with fixed cap
        time_factor = min(1 + (i / 365) * base_uncertainty, 1.055)  # Lower cap

        # Simpler cycle factors
        cycle_position = (current_cycle_position + i / 1460) % 1
        if cycle_position > 0.75:
            cycle_factor = 0.94
        else:
            cycle_factor = 0.96

        adjustments[i] = time_factor * cycle_factor

    return adjustments


def project_prices(
    df, days_forward=365, simulations=1000, confidence_levels=[0.95, 0.68]
):
    """
    Project future Bitcoin prices with simplified calibration.
    """
    df = df.copy()
    df["Log_Price"] = np.log(df["Close"])
    df["Log_Return"] = df["Log_Price"].diff()

    # Get current cycle position
    halving_dates = get_halving_dates()
    current_date = df["Date"].max()
    cycle_position = get_cycle_position(current_date, halving_dates)
    current_cycle_days = int(cycle_position * 4 * 365)

    # Current price and date
    last_price = df["Close"].iloc[-1]
    last_date = df["Date"].iloc[-1]

    # Generate dates for projection
    future_dates = pd.date_range(
        start=last_date + timedelta(days=1), periods=days_forward, freq="D"
    )

    # Calculate expected returns with cycle boundary handling
    future_cycle_days = [
        (current_cycle_days + i) % (4 * 365) for i in range(days_forward)
    ]
    cycle_trends = analyze_trends(df)

    # Get base expected returns
    expected_returns = np.array(
        [cycle_trends.get(day, cycle_trends.mean()) for day in future_cycle_days]
    )

    # Apply trend adjustments
    expected_returns = adjust_trend_expectations(expected_returns, cycle_position)

    # Calculate base volatility
    base_volatility = calculate_volatility(df)

    # Get era adjustments
    current_era = None
    for era, params in era_adjustments.items():
        if (
            df["Date"].min() >= params["start_date"]
            and df["Date"].min() < params["end_date"]
        ):
            current_era = era
            break
    if current_era is None:
        current_era = "mature"

    era_params = era_adjustments[current_era]

    # Get projection adjustments for scaling uncertainty over time
    projection_adjustments = get_projection_adjustments(days_forward, cycle_position)

    # Run Monte Carlo simulation
    np.random.seed(42)
    simulated_paths = np.zeros((days_forward, simulations))

    for sim in range(simulations):
        # Apply era-specific adjustments
        drift = expected_returns * era_params["trend_scale"]
        vol = base_volatility * era_params["volatility_scale"]

        # Scale volatility by projection adjustments
        time_scaled_vol = vol * projection_adjustments

        # Generate returns with time-varying volatility
        returns = np.random.normal(loc=drift, scale=time_scaled_vol, size=days_forward)

        # Calculate price path
        cumulative_returns = np.cumsum(returns)
        price_path = last_price * np.exp(cumulative_returns)
        simulated_paths[:, sim] = price_path

    # Calculate results
    results = pd.DataFrame(index=future_dates)
    results["Median"] = np.percentile(simulated_paths, 50, axis=1)

    # Calculate Expected_Trend using adjusted drift
    cumulative_drift = np.cumsum(drift)
    results["Expected_Trend"] = last_price * np.exp(cumulative_drift)

    # Calculate confidence intervals
    for level in confidence_levels:
        lower_percentile = (1 - level) * 100 / 2
        upper_percentile = 100 - lower_percentile

        results[f"Lower_{int(level*100)}"] = np.percentile(
            simulated_paths, lower_percentile, axis=1
        )
        results[f"Upper_{int(level*100)}"] = np.percentile(
            simulated_paths, upper_percentile, axis=1
        )

    return results


def analyze_bitcoin_prices(csv_path):
    """
    Analyze Bitcoin price data to calculate volatility and growth rates.
    """
    # Read CSV with proper data types
    df = pd.read_csv(csv_path, parse_dates=[0])

    # Print first few rows of raw data to inspect
    print("\nFirst few rows of raw data:")
    print(df.head())

    # Print data info to see types and non-null counts
    print("\nDataset Info:")
    print(df.info())

    # Convert price columns to float and handle any potential formatting issues
    numeric_columns = ["Price", "Open", "High", "Low", "Vol."]  # Added Volume
    for col in numeric_columns:
        # Remove any commas and 'K'/'M' suffixes
        df[col] = df[col].astype(str).str.replace(",", "")
        # Convert K to thousands
        df[col] = df[col].str.replace("K", "e3")
        # Convert M to millions
        df[col] = df[col].str.replace("M", "e6")
        # Convert B to billions
        df[col] = df[col].str.replace("B", "e9")
        # Convert to numeric
        df[col] = pd.to_numeric(df[col], errors="coerce")

    # Rename columns for clarity
    df.columns = ["Date", "Close", "Open", "High", "Low", "Volume", "Change"]

    # Sort by date in ascending order
    df = df.sort_values("Date")

    # Print summary statistics after conversion
    print("\nPrice Summary After Conversion:")
    print(df[["Close", "Open", "High", "Low", "Volume"]].describe())

    # Calculate daily returns
    df["Daily_Return"] = df["Close"].pct_change()

    # Print first few daily returns to verify calculation
    print("\nFirst few daily returns:")
    print(df[["Date", "Close", "Daily_Return"]].head())

    # Check for any infinite or NaN values
    print("\nInfinite or NaN value counts:")
    print(df.isna().sum())

    # Calculate metrics using 365 days for annualization
    analysis = {
        "period_start": df["Date"].min().strftime("%Y-%m-%d"),
        "period_end": df["Date"].max().strftime("%Y-%m-%d"),
        "total_days": len(df),
        "daily_volatility": df["Daily_Return"].std(),
        "annualized_volatility": df["Daily_Return"].std() * np.sqrt(365),
        "total_return": (df["Close"].iloc[-1] / df["Close"].iloc[0] - 1) * 100,
        "average_daily_return": df["Daily_Return"].mean() * 100,
        "average_annual_return": ((1 + df["Daily_Return"].mean()) ** 365 - 1) * 100,
        "min_price": df["Low"].min(),
        "max_price": df["High"].max(),
        "avg_price": df["Close"].mean(),
        "start_price": df["Close"].iloc[0],
        "end_price": df["Close"].iloc[-1],
    }

    # Calculate rolling metrics
    df["Rolling_Volatility_30d"] = df["Daily_Return"].rolling(
        window=30
    ).std() * np.sqrt(365)
    df["Rolling_Return_30d"] = df["Close"].pct_change(periods=30) * 100

    return analysis, df


# Main plotting functions


def create_plots(df, start=None, end=None, project_days=365):
    """
    Create enhanced plots including market maturity visualization.
    """
    # Filter data based on date range
    mask = pd.Series(True, index=df.index)
    if start:
        mask &= df["Date"] >= pd.to_datetime(start)
    if end:
        mask &= df["Date"] <= pd.to_datetime(end)

    plot_df = df[mask].copy()

    if len(plot_df) == 0:
        raise ValueError("No data found for the specified date range")

    # Generate projections
    projections = project_prices(plot_df, days_forward=project_days)

    # Set up the style
    plt.style.use("seaborn-v0_8")

    # Create figure with adjusted size for additional subplot
    fig = plt.figure(figsize=(15, 15))  # Increased height to accommodate new subplot

    # Date range for titles
    hist_date_range = f" ({plot_df['Date'].min().strftime('%Y-%m-%d')} to {plot_df['Date'].max().strftime('%Y-%m-%d')})"

    # 1. Price history and projections (log scale)
    ax1 = plt.subplot(4, 1, 1)

    # Plot historical prices
    ax1.semilogy(plot_df["Date"], plot_df["Close"], "b-", label="Historical Price")

    # Plot projections
    ax1.semilogy(
        projections.index,
        projections["Expected_Trend"],
        "--",
        color="purple",
        label="Expected Trend",
    )
    ax1.semilogy(
        projections.index,
        projections["Median"],
        ":",
        color="green",
        label="Simulated Median",
    )
    ax1.fill_between(
        projections.index,
        projections["Lower_95"],
        projections["Upper_95"],
        alpha=0.2,
        color="orange",
        label="95% Confidence Interval",
    )
    ax1.fill_between(
        projections.index,
        projections["Lower_68"],
        projections["Upper_68"],
        alpha=0.3,
        color="green",
        label="68% Confidence Interval",
    )

    # Customize y-axis
    ax1.yaxis.set_major_formatter(plt.FuncFormatter(format_price))
    min_price = min(plot_df["Low"].min(), projections["Lower_95"].min())
    max_price = max(plot_df["High"].max(), projections["Upper_95"].max())
    price_points = get_nice_price_points(min_price, max_price)
    ax1.set_yticks(price_points)
    ax1.tick_params(axis="y", labelsize=8)
    ax1.margins(y=0.02)
    ax1.grid(True, which="major", linestyle="-", alpha=0.5)
    ax1.grid(True, which="minor", linestyle=":", alpha=0.2)
    ax1.set_title("Bitcoin Price History and Projections (Log Scale)" + hist_date_range)
    ax1.legend(fontsize=8)

    # 3. Rolling volatility
    ax3 = plt.subplot(5, 1, 2)
    ax3.plot(
        plot_df["Date"],
        plot_df["Rolling_Volatility_30d"],
        "r-",
        label="30-Day Rolling Volatility",
    )
    ax3.set_title("30-Day Rolling Volatility (Annualized)" + hist_date_range)
    ax3.set_ylabel("Volatility")
    ax3.grid(True)
    ax3.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: "{:.0%}".format(y)))
    ax3.legend()

    # 4. Returns distribution
    ax4 = plt.subplot(5, 1, 3)
    returns_mean = plot_df["Daily_Return"].mean()
    returns_std = plot_df["Daily_Return"].std()
    filtered_returns = plot_df["Daily_Return"][
        (plot_df["Daily_Return"] > returns_mean - 5 * returns_std)
        & (plot_df["Daily_Return"] < returns_mean + 5 * returns_std)
    ]

    sns.histplot(filtered_returns, bins=100, ax=ax4)
    ax4.set_title(
        "Distribution of Daily Returns (Excluding Extreme Outliers)" + hist_date_range
    )
    ax4.set_xlabel("Daily Return")
    ax4.set_ylabel("Count")
    ax4.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: "{:.0%}".format(x)))

    # Add mean line
    ax4.axvline(filtered_returns.mean(), color="r", linestyle="dashed", linewidth=1)
    ax4.text(
        filtered_returns.mean(),
        ax4.get_ylim()[1],
        "Mean",
        rotation=90,
        va="top",
        ha="right",
    )

    # 5. Projection ranges
    ax5 = plt.subplot(5, 1, 4)
    timepoints = np.array(range(30, project_days, 30))
    timepoints = timepoints[timepoints <= project_days]

    ranges = []
    labels = []
    positions = []

    for t in timepoints:
        idx = t - 1
        ranges.extend(
            [
                projections["Lower_95"].iloc[idx],
                projections["Lower_68"].iloc[idx],
                projections["Median"].iloc[idx],
                projections["Upper_68"].iloc[idx],
                projections["Upper_95"].iloc[idx],
            ]
        )
        labels.extend(["95% Lower", "68% Lower", "Median", "68% Upper", "95% Upper"])
        positions.extend([t] * 5)

    ax5.scatter(positions, ranges, alpha=0.6)

    for t in timepoints:
        idx = positions.index(t)
        ax5.plot([t] * 5, ranges[idx : idx + 5], "k-", alpha=0.3)

    ax5.set_yscale("log")
    min_price = min(ranges)
    max_price = max(ranges)
    price_points = get_nice_price_points(min_price, max_price)
    ax5.set_yticks(price_points)
    ax5.yaxis.set_major_formatter(plt.FuncFormatter(format_price))
    ax5.set_title("Projected Price Ranges at Future Timepoints")
    ax5.set_xlabel("Days Forward")
    ax5.set_ylabel("Price (USD)")
    ax5.grid(True, alpha=0.3)
    ax5.set_xticks(timepoints)

    # Adjust layout
    plt.tight_layout(h_pad=1.0)  # Increased spacing between subplots

    # Save the plot
    start_str = start if start else plot_df["Date"].min().strftime("%Y-%m-%d")
    end_str = end if end else plot_df["Date"].max().strftime("%Y-%m-%d")
    filename = f"bitcoin_analysis_{start_str}_to_{end_str}_with_projections.png"
    plt.savefig(filename, dpi=300, bbox_inches="tight")
    plt.close()

    return projections


def visualize_cycle_patterns(df, cycle_returns, cycle_volatility):
    """
    Create enhanced visualization of Bitcoin's behavior across halving cycles.
    """
    plt.style.use("seaborn-v0_8")
    fig = plt.figure(figsize=(15, 15))

    # Create a 3x1 subplot grid with different heights
    gs = plt.GridSpec(3, 1, height_ratios=[2, 1, 2], hspace=0.3)

    # Plot 1: Returns across cycle with confidence bands
    ax1 = plt.subplot(gs[0])

    # Convert days to percentage through cycle
    x_points = np.array(cycle_returns.index) / (4 * 365) * 100

    # Calculate rolling mean and standard deviation for confidence bands
    window = 30  # 30-day window
    rolling_mean = pd.Series(cycle_returns.values).rolling(window=window).mean()
    rolling_std = pd.Series(cycle_returns.values).rolling(window=window).std()

    # Plot confidence bands
    ax1.fill_between(
        x_points,
        (rolling_mean - 2 * rolling_std) * 100,
        (rolling_mean + 2 * rolling_std) * 100,
        alpha=0.2,
        color="blue",
        label="95% Confidence",
    )
    ax1.fill_between(
        x_points,
        (rolling_mean - rolling_std) * 100,
        (rolling_mean + rolling_std) * 100,
        alpha=0.3,
        color="blue",
        label="68% Confidence",
    )

    # Plot average returns
    ax1.plot(
        x_points,
        cycle_returns.values * 100,
        "b-",
        label="Average Daily Return",
        linewidth=2,
    )
    ax1.axhline(y=0, color="gray", linestyle="--", alpha=0.5)

    # Add vertical lines for each year in cycle
    for year in range(1, 4):
        ax1.axvline(x=year * 25, color="gray", linestyle=":", alpha=0.3)
        ax1.text(
            year * 25,
            ax1.get_ylim()[1],
            f"Year {year}",
            rotation=90,
            va="top",
            ha="right",
            alpha=0.7,
        )

    # Highlight halving points
    ax1.axvline(x=0, color="red", linestyle="--", alpha=0.5, label="Halving Event")
    ax1.axvline(x=100, color="red", linestyle="--", alpha=0.5)

    ax1.set_title("Bitcoin Return Patterns Across Halving Cycle", pad=20)
    ax1.set_xlabel("Position in Cycle (%)")
    ax1.set_ylabel("Average Daily Return (%)")
    ax1.grid(True, alpha=0.3)
    ax1.legend(loc="upper right")

    # Plot 2: Volatility across cycle
    ax2 = plt.subplot(gs[1])

    # Calculate rolling volatility confidence bands
    vol_mean = pd.Series(cycle_volatility.values).rolling(window=window).mean()
    vol_std = pd.Series(cycle_volatility.values).rolling(window=window).std()

    # Plot volatility with confidence bands
    annualized_factor = np.sqrt(365) * 100
    ax2.fill_between(
        x_points,
        (vol_mean - 2 * vol_std) * annualized_factor,
        (vol_mean + 2 * vol_std) * annualized_factor,
        alpha=0.2,
        color="red",
        label="95% Confidence",
    )
    ax2.plot(
        x_points,
        cycle_volatility.values * annualized_factor,
        "r-",
        label="Annualized Volatility",
        linewidth=2,
    )

    # Add year markers
    for year in range(1, 4):
        ax2.axvline(x=year * 25, color="gray", linestyle=":", alpha=0.3)

    ax2.axvline(x=0, color="red", linestyle="--", alpha=0.5)
    ax2.axvline(x=100, color="red", linestyle="--", alpha=0.5)

    ax2.set_xlabel("Position in Cycle (%)")
    ax2.set_ylabel("Volatility (%)")
    ax2.grid(True, alpha=0.3)
    ax2.legend(loc="upper right")

    # Plot 3: Average price trajectory within cycles
    ax3 = plt.subplot(gs[2])

    # Define a color scheme for cycles
    cycle_colors = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd"]

    # Calculate average price path for each cycle
    halving_dates = get_halving_dates()
    cycles = []

    for i in range(len(halving_dates) - 1):
        cycle_start = halving_dates[i]
        cycle_end = halving_dates[i + 1]
        cycle_data = df[(df["Date"] >= cycle_start) & (df["Date"] < cycle_end)].copy()

        if len(cycle_data) > 0:
            cycle_data["Cycle_Pct"] = (
                (cycle_data["Date"] - cycle_start).dt.total_seconds()
                / (cycle_end - cycle_start).total_seconds()
                * 100
            )
            cycle_data["Normalized_Price"] = (
                cycle_data["Close"] / cycle_data["Close"].iloc[0]
            )
            cycles.append(cycle_data)

    # Plot each historical cycle with distinct colors
    for i, cycle in enumerate(cycles):
        ax3.semilogy(
            cycle["Cycle_Pct"],
            cycle["Normalized_Price"],
            color=cycle_colors[i],
            alpha=0.7,
            label=f'Cycle {i+1} ({cycle["Date"].iloc[0].strftime("%Y")}-{cycle["Date"].iloc[-1].strftime("%Y")})',
        )

    # Calculate and plot average cycle
    if cycles:
        avg_cycle = pd.concat(
            [c.set_index("Cycle_Pct")["Normalized_Price"] for c in cycles], axis=1
        )
        avg_cycle_mean = avg_cycle.mean(axis=1)
        avg_cycle_std = avg_cycle.std(axis=1)

        ax3.semilogy(
            avg_cycle_mean.index,
            avg_cycle_mean.values,
            "k-",
            linewidth=2,
            label="Average Cycle",
        )
        ax3.fill_between(
            avg_cycle_mean.index,
            avg_cycle_mean * np.exp(-2 * avg_cycle_std),
            avg_cycle_mean * np.exp(2 * avg_cycle_std),
            alpha=0.2,
            color="gray",
        )

    # Add year markers
    for year in range(1, 4):
        ax3.axvline(x=year * 25, color="gray", linestyle=":", alpha=0.3)

    ax3.axvline(x=0, color="red", linestyle="--", alpha=0.5)
    ax3.axvline(x=100, color="red", linestyle="--", alpha=0.5)

    ax3.set_title("Price Performance Across Cycles (Normalized)", pad=20)
    ax3.set_xlabel("Position in Cycle (%)")
    ax3.set_ylabel("Price (Relative to Cycle Start)")
    ax3.grid(True, alpha=0.3)
    ax3.legend(loc="center left", bbox_to_anchor=(1.02, 0.5))

    # Add current cycle position marker on all plots
    current_position = get_cycle_position(df["Date"].max(), halving_dates) * 100
    for ax in [ax1, ax2, ax3]:
        ax.axvline(
            x=current_position,
            color="green",
            linestyle="-",
            alpha=0.5,
            label="Current Position",
        )

    # Main title for the figure
    fig.suptitle("Bitcoin Halving Cycle Analysis", fontsize=16, y=0.95)

    # Adjust layout to prevent legend cutoff
    plt.tight_layout()

    # Save the plot
    plt.savefig("bitcoin_cycle_patterns.png", dpi=300, bbox_inches="tight")
    plt.close()


def create_backtest_plot(
    df, backtest_date="2020-05-11", start_date="2012-11-28", project_days=1650
):
    """
    Create a plot comparing actual price history against model projections from a historical date.
    Returns both the projections and performance metrics.

    Args:
        df: DataFrame with historical price data
        backtest_date: Date to start the backtest from
        start_date: Date to start considering historical data
        project_days: Number of days to project forward from backtest date

    Returns:
        tuple: (projections DataFrame, metrics dictionary)
    """
    # Convert dates to datetime
    backtest_date = pd.to_datetime(backtest_date)
    start_date = pd.to_datetime(start_date)

    # Validate dates
    if start_date >= backtest_date:
        raise ValueError("start_date must be earlier than backtest_date")

    # Clean the data: remove rows with zero or invalid prices and filter by date
    df = df[(df["Close"] > 0) & (df["Date"] >= start_date)].copy()

    # Split data into training (before backtest date) and validation (after backtest date)
    training_df = df[df["Date"] <= backtest_date].copy()
    validation_df = df[df["Date"] > backtest_date].copy()

    # Check if we have enough data
    if len(training_df) < 30:  # Require at least 30 days of training data
        raise ValueError("Insufficient training data before backtest date")

    # Generate historical projections using only training data
    historical_projections = project_prices(training_df, days_forward=project_days)

    # Set up the plot
    plt.style.use("seaborn-v0_8")
    fig, ax = plt.figure(figsize=(15, 10)), plt.gca()

    # Plot training data
    heading_label = f'Historical Price (Training: {start_date.strftime("%Y-%m-%d")} to {backtest_date.strftime("%Y-%m-%d")})'

    ax.semilogy(
        training_df["Date"],
        training_df["Close"],
        "b-",
        label=heading_label,
        alpha=0.7,
    )

    # Plot validation data
    ax.semilogy(
        validation_df["Date"],
        validation_df["Close"],
        "g-",
        label=f'Actual Price (Validation: {backtest_date.strftime("%Y-%m-%d")} onwards)',
        linewidth=2,
    )

    # Plot projections
    ax.semilogy(
        historical_projections.index,
        historical_projections["Expected_Trend"],
        "--",
        color="purple",
        label="Model Projection (Expected)",
    )
    ax.semilogy(
        historical_projections.index,
        historical_projections["Median"],
        ":",
        color="orange",
        label="Model Projection (Median)",
    )

    # Add confidence intervals
    ax.fill_between(
        historical_projections.index,
        historical_projections["Lower_95"],
        historical_projections["Upper_95"],
        alpha=0.2,
        color="orange",
        label="95% Confidence Interval",
    )
    ax.fill_between(
        historical_projections.index,
        historical_projections["Lower_68"],
        historical_projections["Upper_68"],
        alpha=0.3,
        color="green",
        label="68% Confidence Interval",
    )

    # Customize y-axis
    ax.yaxis.set_major_formatter(plt.FuncFormatter(format_price))

    # Set custom y-axis ticks
    min_price = min(
        df["Low"].min(),
        historical_projections["Lower_95"].min(),
        0.0001,  # Set minimum price floor
    )
    max_price = max(df["High"].max(), historical_projections["Upper_95"].max())
    price_points = get_nice_price_points(min_price, max_price)
    ax.set_yticks(price_points)

    # Add halving lines
    halving_dates = get_halving_dates()
    relevant_halvings = halving_dates[
        (halving_dates >= start_date) & (halving_dates <= validation_df["Date"].max())
    ]
    for date in relevant_halvings:
        ax.axvline(date, color="red", linestyle="--", alpha=0.3)
        ax.text(
            date,
            ax.get_ylim()[1],
            "Halving",
            rotation=90,
            va="top",
            ha="right",
            alpha=0.7,
        )

    # Calculate model performance metrics
    metrics = {}
    if len(validation_df) > 0:
        # Create a common date range for comparison
        actual_prices = validation_df.set_index("Date")["Close"]
        common_dates = actual_prices.index.intersection(historical_projections.index)

        if len(common_dates) > 0:
            actual_aligned = actual_prices[common_dates]
            projections_aligned = historical_projections.loc[common_dates]

            # Calculate metrics using aligned data
            metrics = {
                "mape": np.mean(
                    np.abs(
                        (actual_aligned - projections_aligned["Expected_Trend"])
                        / actual_aligned
                    )
                )
                * 100,
                "rmse": np.sqrt(
                    np.mean(
                        (actual_aligned - projections_aligned["Expected_Trend"]) ** 2
                    )
                ),
                "max_error": np.max(
                    np.abs(actual_aligned - projections_aligned["Expected_Trend"])
                ),
                "coverage_95": np.mean(
                    (actual_aligned >= projections_aligned["Lower_95"])
                    & (actual_aligned <= projections_aligned["Upper_95"])
                )
                * 100,
                "coverage_68": np.mean(
                    (actual_aligned >= projections_aligned["Lower_68"])
                    & (actual_aligned <= projections_aligned["Upper_68"])
                )
                * 100,
            }

            # Add metrics to plot
            metrics_text = (
                f"Model Performance Metrics:\n"
                f"MAPE: {metrics['mape']:.1f}%\n"
                f"RMSE: ${metrics['rmse']:,.0f}\n"
                f"Max Error: ${metrics['max_error']:,.0f}\n"
                f"95% CI Coverage: {metrics['coverage_95']:.1f}%\n"
                f"68% CI Coverage: {metrics['coverage_68']:.1f}%"
            )
            ax.text(
                0.02,
                0.98,
                metrics_text,
                transform=ax.transAxes,
                verticalalignment="top",
                bbox=dict(facecolor="white", alpha=0.8),
            )

    # Customize plot
    ax.set_title(
        f'Bitcoin Price: Model Backtest\nTraining: {start_date.strftime("%Y-%m-%d")} to {backtest_date.strftime("%Y-%m-%d")}'
    )
    ax.set_xlabel("Date")
    ax.set_ylabel("Price (USD)")
    ax.grid(True, which="major", linestyle="-", alpha=0.5)
    ax.grid(True, which="minor", linestyle=":", alpha=0.2)
    ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5))

    # Adjust layout and save
    plt.tight_layout()
    filename = f'bitcoin_backtest_{start_date.strftime("%Y%m%d")}_to_{backtest_date.strftime("%Y%m%d")}.png'
    plt.savefig(filename, dpi=300, bbox_inches="tight")
    plt.close()

    return historical_projections, metrics


def run_projection(args):
    df, start = args
    projections = create_plots(df, start=start, project_days=365 * 4)


def run_projections(df):
    # Create main projection
    projection_starts = [
        "2013-01-01",
        "2014-01-01",
        "2015-01-01",
        "2016-07-09",
    ]
    args = [(df, start) for start in projection_starts]
    with Pool() as pool:
        pool.map(run_projection, args)


def run_single_backtest(args):
    """
    Run a single backtest with the given parameters.
    Must be defined at module level for multiprocessing.

    Args:
        args: tuple of (params dict, DataFrame)
    """
    params, df = args
    try:
        # Create a copy of params without the description
        backtest_params = params.copy()
        backtest_params.pop("description", None)

        projections, metrics = create_backtest_plot(df, **backtest_params)

        # Ensure metrics has all required keys with default values
        if metrics is None:
            metrics = {}

        default_metrics = {
            "mape": 0.0,
            "rmse": 0.0,
            "max_error": 0.0,
            "coverage_95": 0.0,
            "coverage_68": 0.0,
        }

        # Update metrics with defaults for any missing keys
        metrics = {**default_metrics, **metrics}

        return {
            "params": params,
            "projections": projections,
            "metrics": metrics,
            "success": True,
        }
    except Exception as e:
        print(
            f"Error in backtest for period {params['description']}: {str(e)}"
        )  # Debug print
        return {"params": params, "error": str(e), "success": False}


def run_systematic_backtests(df, validation_years=1, min_training_years=8):
    """
    Run a comprehensive suite of backtests with consistent validation periods.
    Uses sliding windows for both start and end dates.
    """
    # Convert years to days
    validation_days = validation_years * 365
    min_training_days = min_training_years * 365

    # Define start date for reliable data
    mature_start = pd.Timestamp("2013-01-01")
    last_possible_start = df["Date"].max() - pd.Timedelta(
        days=min_training_days + validation_days
    )
    end_date = df["Date"].max() - pd.Timedelta(days=validation_days)

    if mature_start >= last_possible_start:
        raise ValueError(
            f"Insufficient data for backtesting with current parameters:\n"
            f"- Data range: {mature_start} to {df['Date'].max()}\n"
            f"- Minimum training period: {min_training_years} years\n"
            f"- Validation period: {validation_years} years"
        )

    old_backtests = [
        {
            "start_date": "2016-07-09",
            "backtest_date": "2024-04-19",
            "project_days": validation_days,
            "description": "Second until fourth halving",
        },
        {
            "start_date": "2013-01-01",  # Includes pre-futures for cycle learning
            "backtest_date": "2020-05-11",
            "project_days": validation_days,
            "description": "Post-Futures Window with two cycles of training",
        },
        {
            "start_date": "2014-01-01",
            "backtest_date": "2021-12-31",
            "project_days": validation_days,
            "description": "Cross-Regime Test with two cycles of training",
        },
        {
            "start_date": "2015-01-01",
            "backtest_date": "2022-01-01",
            "project_days": validation_days,
            "description": "Recent Window focusing on post-2022 behavior",
        },
    ]
    backtest_periods = []
    backtest_periods.extend(old_backtests)

    # Generate backtest periods with sliding windows
    window_start = mature_start
    step = pd.Timedelta(days=180)  # 6 month steps

    while window_start <= last_possible_start:
        backtest_date = window_start + pd.Timedelta(days=min_training_days)

        backtest_periods.append(
            {
                "start_date": window_start.strftime("%Y-%m-%d"),
                "backtest_date": backtest_date.strftime("%Y-%m-%d"),
                "project_days": validation_days,
                "description": f"Training {window_start.strftime('%Y-%m-%d')} to {backtest_date.strftime('%Y-%m-%d')}",
            }
        )
        window_start += step

    # Add specific periods of interest
    special_periods = []

    # Halving-based periods
    halving_dates = get_halving_dates()
    relevant_halvings = [
        h
        for h in halving_dates
        if h < end_date and h > (mature_start + pd.Timedelta(days=min_training_days))
    ]

    for halving in relevant_halvings:
        earliest_start = halving - pd.Timedelta(days=min_training_days)
        if earliest_start >= mature_start:
            special_periods.append(
                {
                    "start_date": earliest_start.strftime("%Y-%m-%d"),
                    "backtest_date": halving.strftime("%Y-%m-%d"),
                    "project_days": validation_days,
                    "description": f"Pre-halving {halving.strftime('%Y')}",
                }
            )

    # Market structure change periods
    important_dates = [
        ("2017-12-01", "Post-futures introduction"),
        ("2020-03-01", "Post-COVID crash"),
        ("2021-11-01", "Post-2021 peak"),
    ]

    for date, description in important_dates:
        test_date = pd.Timestamp(date)
        if test_date < end_date:
            earliest_start = test_date - pd.Timedelta(days=min_training_days)
            if earliest_start >= mature_start:
                special_periods.append(
                    {
                        "start_date": earliest_start.strftime("%Y-%m-%d"),
                        "backtest_date": date,
                        "project_days": validation_days,
                        "description": description,
                    }
                )

    # Combine and remove any duplicates
    all_periods = backtest_periods + special_periods
    unique_periods = []
    seen_dates = set()
    for period in all_periods:
        key = f"{period['start_date']}_{period['backtest_date']}"
        if key not in seen_dates:
            unique_periods.append(period)
            seen_dates.add(key)

    if not unique_periods:
        raise ValueError("No valid backtest periods found with current parameters")

    # Sort periods by backtest date for clearer analysis
    unique_periods.sort(key=lambda x: pd.Timestamp(x["backtest_date"]))

    print(f"\nRunning backtests with:")
    print(
        f"- Start dates range: {unique_periods[0]['start_date']} to {unique_periods[-1]['start_date']}"
    )
    print(
        f"- Backtest dates range: {unique_periods[0]['backtest_date']} to {unique_periods[-1]['backtest_date']}"
    )
    print(f"- Minimum training period: {min_training_years} years")
    print(f"- Validation period: {validation_years} years")
    print(f"- Number of test periods: {len(unique_periods)}")

    print("\nTest periods:")
    for period in unique_periods:
        print(f"- {period['description']}")

    # Create args tuples with params and DataFrame
    args = [(params, df) for params in unique_periods]

    # Use multiprocessing
    with Pool() as pool:
        results = pool.map(run_single_backtest, args)

    # Analyze results
    successful_tests = [r for r in results if r["success"]]
    failed_tests = [r for r in results if not r["success"]]

    # Define stress periods
    stress_periods = {
        # COVID crash and recovery
        ("2020-03-01", "2020-09-01"): "COVID crash period",
        # 2021 peak and subsequent crash
        ("2021-11-01", "2022-06-01"): "2021 peak aftermath",
        # Add more stress periods as needed
    }

    def is_stress_period(test_date):
        """Check if a test date falls in any stress period"""
        test_date = pd.Timestamp(test_date)
        for (start, end), _ in stress_periods.items():
            if pd.Timestamp(start) <= test_date <= pd.Timestamp(end):
                return True
        return False

    # Categorize results
    normal_periods = []
    stress_periods_results = []

    for result in successful_tests:
        if is_stress_period(result["params"]["backtest_date"]):
            stress_periods_results.append(result)
        else:
            normal_periods.append(result)

    # Calculate metrics for each category
    def calculate_category_metrics(results):
        if not results:
            return None
        return {
            "count": len(results),
            "mape": np.mean([r["metrics"]["mape"] for r in results]),
            "rmse": np.mean([r["metrics"]["rmse"] for r in results]),
            "max_error": np.mean([r["metrics"]["max_error"] for r in results]),
            "coverage_95": np.mean([r["metrics"]["coverage_95"] for r in results]),
            "coverage_68": np.mean([r["metrics"]["coverage_68"] for r in results]),
        }

    normal_metrics = calculate_category_metrics(normal_periods)
    stress_metrics = calculate_category_metrics(stress_periods_results)

    # Write detailed results
    with open("bitcoin_backtest_results_summary.txt", "w") as f:
        f.write("Systematic Backtest Results\n")
        f.write("==========================\n\n")

        f.write("Configuration:\n")
        f.write(f"- Minimum training period: {min_training_years} years\n")
        f.write(f"- Validation period: {validation_years} years\n")
        f.write(
            f"- Start dates range: {unique_periods[0]['start_date']} to {unique_periods[-1]['start_date']}\n"
        )
        f.write(
            f"- Backtest dates range: {unique_periods[0]['backtest_date']} to {unique_periods[-1]['backtest_date']}\n"
        )
        f.write(f"- Number of test periods: {len(unique_periods)}\n\n")

        # Normal Periods
        f.write("Normal Market Periods\n")
        f.write("====================\n")
        f.write(f"Number of periods: {len(normal_periods)}\n\n")

        for result in normal_periods:
            f.write("\n" + "=" * 50 + "\n")
            f.write(f"Period: {result['params']['description']}\n")
            f.write(
                f"Training: {result['params']['start_date']} to {result['params']['backtest_date']}\n"
            )
            f.write(
                f"Validation: {result['params']['backtest_date']} to {pd.Timestamp(result['params']['backtest_date']) + pd.Timedelta(days=validation_years*365):%Y-%m-%d}\n"
            )
            f.write("\nMetrics:\n")
            for metric, value in result["metrics"].items():
                if metric in ["mape", "coverage_95", "coverage_68"]:
                    f.write(f"- {metric}: {value:.1f}%\n")
                else:
                    f.write(f"- {metric}: ${value:,.0f}\n")
            f.write("\n")

        if normal_metrics:
            f.write("\nNormal Periods Aggregate Metrics:\n")
            f.write(f"MAPE: {normal_metrics['mape']:.1f}%\n")
            f.write(f"RMSE: ${normal_metrics['rmse']:,.0f}\n")
            f.write(f"Average Max Error: ${normal_metrics['max_error']:,.0f}\n")
            f.write(f"95% CI Coverage: {normal_metrics['coverage_95']:.1f}%\n")
            f.write(f"68% CI Coverage: {normal_metrics['coverage_68']:.1f}%\n")

        # Stress Periods
        f.write("\n\nStress Periods\n")
        f.write("=============\n")
        f.write(f"Number of periods: {len(stress_periods_results)}\n\n")

        for result in stress_periods_results:
            f.write("\n" + "=" * 50 + "\n")
            f.write(f"Period: {result['params']['description']}\n")
            f.write(
                f"Training: {result['params']['start_date']} to {result['params']['backtest_date']}\n"
            )
            f.write(
                f"Validation: {result['params']['backtest_date']} to {pd.Timestamp(result['params']['backtest_date']) + pd.Timedelta(days=validation_years*365):%Y-%m-%d}\n"
            )
            f.write("\nMetrics:\n")
            for metric, value in result["metrics"].items():
                if metric in ["mape", "coverage_95", "coverage_68"]:
                    f.write(f"- {metric}: {value:.1f}%\n")
                else:
                    f.write(f"- {metric}: ${value:,.0f}\n")
            f.write("\n")

        if stress_metrics:
            f.write("\nStress Periods Aggregate Metrics:\n")
            f.write(f"MAPE: {stress_metrics['mape']:.1f}%\n")
            f.write(f"RMSE: ${stress_metrics['rmse']:,.0f}\n")
            f.write(f"Average Max Error: ${stress_metrics['max_error']:,.0f}\n")
            f.write(f"95% CI Coverage: {stress_metrics['coverage_95']:.1f}%\n")
            f.write(f"68% CI Coverage: {stress_metrics['coverage_68']:.1f}%\n")

    return (
        normal_metrics,
        stress_metrics,
        normal_periods,
        stress_periods_results,
        failed_tests,
    )


# if __name__ == "__main__":
#    analysis, df = analyze_bitcoin_prices("prices.csv")
#    procs = []
#
#    for proc in procs:
#        proc.join()
#


if __name__ == "__main__":
    analysis, df = analyze_bitcoin_prices("prices.csv")
    run_projections(df)
    normal_metrics, stress_metrics, normal_results, stress_results, failed_tests = (
        run_systematic_backtests(df)
    )

    print("\nAggregate Metrics:")
    print(f"Total backtests run: {normal_metrics['count']}")
    print(f"Successful tests: {len(normal_results)}")
    print(f"Failed tests: {len(failed_tests)}")
    print("\nAverage Performance:")
    print(f"MAPE: {normal_metrics['mape']:.1f}%")
    print(f"RMSE: ${normal_metrics['rmse']:,.0f}")
    print(f"95% CI Coverage: {normal_metrics['coverage_95']:.1f}%")
    print(f"68% CI Coverage: {normal_metrics['coverage_68']:.1f}%")

    print("\nFailed Tests:")
    for test in failed_tests:
        print(f"Period: {test['params']['description']}")
        print(f"Error: {test['error']}\n")