1625 lines
55 KiB
Python
1625 lines
55 KiB
Python
import pandas as pd
|
|
import numpy as np
|
|
from datetime import timedelta
|
|
import matplotlib.pyplot as plt
|
|
import seaborn as sns
|
|
from multiprocessing import 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_adaptive_volatility(
|
|
df,
|
|
short_window=30,
|
|
medium_window=90,
|
|
long_window=180,
|
|
vol_clip_min=0.5,
|
|
vol_clip_max=2.0,
|
|
):
|
|
"""
|
|
Calculate volatility with adaptive window sizes based on market conditions.
|
|
Returns a single volatility value for the most recent period.
|
|
|
|
Incorporates long-term volatility as a stability baseline and additional
|
|
reference point for regime detection.
|
|
"""
|
|
df = df.copy()
|
|
df["Log_Return"] = np.log(df["Close"]).diff()
|
|
|
|
# Remove any NaN values that could cause issues
|
|
df = df.dropna()
|
|
|
|
if len(df) < long_window:
|
|
# Not enough data, fall back to simple volatility
|
|
return df["Log_Return"].std()
|
|
|
|
# Get recent data for efficiency
|
|
lookback = max(long_window * 2, 360) # Use enough data for stable estimates
|
|
recent_df = df.iloc[-lookback:].copy() if len(df) > lookback else df.copy()
|
|
|
|
try:
|
|
# Initial volatility estimate using base windows
|
|
short_vol = recent_df["Log_Return"].ewm(span=short_window, adjust=False).std()
|
|
medium_vol = recent_df["Log_Return"].ewm(span=medium_window, adjust=False).std()
|
|
long_vol = recent_df["Log_Return"].ewm(span=long_window, adjust=False).std()
|
|
|
|
# Ensure we have valid volatility values
|
|
if short_vol.iloc[-1] == 0 or np.isnan(short_vol.iloc[-1]):
|
|
return df["Log_Return"].std() # Fallback to simple volatility
|
|
|
|
# Calculate regime indicators for recent period
|
|
medium_vol_mean = medium_vol.rolling(min(90, len(recent_df))).mean()
|
|
long_vol_mean = long_vol.rolling(min(180, len(recent_df))).mean()
|
|
|
|
if medium_vol_mean.iloc[-1] == 0:
|
|
vol_regime = pd.Series([1.0] * len(recent_df))
|
|
else:
|
|
# Compare short-term to both medium and long-term volatility
|
|
medium_regime = short_vol / medium_vol_mean
|
|
long_regime = short_vol / long_vol_mean
|
|
|
|
# Use the more conservative (higher) regime indicator
|
|
vol_regime = pd.concat([medium_regime, long_regime], axis=1).max(axis=1)
|
|
|
|
vol_regime = vol_regime.clip(vol_clip_min, vol_clip_max)
|
|
|
|
# Get most recent regime reading
|
|
latest_regime = vol_regime.iloc[-1]
|
|
|
|
# Adjust window sizes based on current regime
|
|
adj_factor = 1 / latest_regime
|
|
adj_short = max(10, int(short_window * adj_factor)) # Minimum window of 10
|
|
adj_medium = max(30, int(medium_window * adj_factor))
|
|
adj_long = max(60, int(long_window * adj_factor))
|
|
|
|
# Calculate final volatilities using adjusted windows
|
|
final_short = recent_df["Log_Return"].iloc[-adj_short:].std()
|
|
final_medium = recent_df["Log_Return"].iloc[-adj_medium:].std()
|
|
final_long = recent_df["Log_Return"].iloc[-adj_long:].std()
|
|
|
|
# If any volatility measure is NaN or 0, fall back to simple volatility
|
|
if np.isnan([final_short, final_medium, final_long]).any() or 0 in [
|
|
final_short,
|
|
final_medium,
|
|
final_long,
|
|
]:
|
|
return df["Log_Return"].std()
|
|
|
|
# Calculate regime-based weights, now incorporating long-term volatility
|
|
high_vol_weight = (latest_regime - vol_clip_min) / (vol_clip_max - vol_clip_min)
|
|
base_weights = np.array([0.2, 0.5, 0.3]) # Short, medium, long weights
|
|
stress_weights = np.array(
|
|
[0.4, 0.4, 0.2]
|
|
) # More weight on short-term during stress
|
|
|
|
# Interpolate between base and stress weights
|
|
weights = (
|
|
base_weights * (1 - high_vol_weight) + stress_weights * high_vol_weight
|
|
)
|
|
|
|
# Calculate final volatility using all three timeframes
|
|
final_vol = (
|
|
final_short * weights[0]
|
|
+ final_medium * weights[1]
|
|
+ final_long * weights[2]
|
|
)
|
|
|
|
# Add uncertainty adjustment based on regime changes
|
|
regime_change = abs(vol_regime.diff()).fillna(0)
|
|
regime_change_mean = regime_change.rolling(5, min_periods=1).mean().iloc[-1]
|
|
if regime_change_mean == 0:
|
|
uncertainty_adjustment = 1.0
|
|
else:
|
|
regime_change_zscore = regime_change.iloc[-1] / regime_change_mean
|
|
uncertainty_adjustment = 1 + np.clip(regime_change_zscore / 2, 0, 0.5)
|
|
|
|
return max(final_vol * uncertainty_adjustment, df["Log_Return"].std() * 0.5)
|
|
|
|
except Exception as e:
|
|
print(f"Error in adaptive volatility calculation: {e}")
|
|
# Fall back to simple volatility calculation
|
|
return df["Log_Return"].std()
|
|
|
|
|
|
def calculate_volatility(df, short_window=30, medium_window=90, long_window=180):
|
|
"""
|
|
Calculate volatility using adaptive windows and era-specific adjustments.
|
|
Returns a single volatility value.
|
|
"""
|
|
if len(df) < 30:
|
|
return 0.02 # Return a reasonable default for very short periods
|
|
|
|
try:
|
|
# Calculate adaptive volatility
|
|
base_vol = calculate_adaptive_volatility(
|
|
df,
|
|
short_window=short_window,
|
|
medium_window=medium_window,
|
|
long_window=long_window,
|
|
)
|
|
|
|
if np.isnan(base_vol) or base_vol == 0:
|
|
base_vol = df["Close"].pct_change().std()
|
|
|
|
# Era-specific adjustments
|
|
start_date = df["Date"].min()
|
|
if start_date >= pd.Timestamp("2020-01-01"):
|
|
base_adjustment = 0.64
|
|
elif start_date >= pd.Timestamp("2016-07-09"):
|
|
base_adjustment = 0.67
|
|
elif start_date >= pd.Timestamp("2015-01-01"):
|
|
base_adjustment = 0.69
|
|
else:
|
|
base_adjustment = 0.70
|
|
|
|
return max(
|
|
base_vol * base_adjustment, 0.01
|
|
) # Ensure we never return 0 volatility
|
|
|
|
except Exception as e:
|
|
print(f"Error in volatility calculation: {e}")
|
|
# Fall back to simple volatility with minimum floor
|
|
return max(df["Close"].pct_change().std(), 0.01)
|
|
|
|
|
|
# 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 calculate_market_conditions(df, lookback_window=180):
|
|
"""
|
|
Calculate market condition metrics to inform uncertainty scaling.
|
|
"""
|
|
df = df.copy() # Avoid modifying original dataframe
|
|
metrics = {}
|
|
|
|
# Use log returns for stability
|
|
df["Log_Return"] = np.log(df["Close"]).diff()
|
|
|
|
# Handle initial NaN values
|
|
df["Log_Return"] = df["Log_Return"].fillna(method="bfill")
|
|
|
|
# Recent vs historical volatility ratio
|
|
recent_vol = max(df["Log_Return"].tail(30).std(), 1e-6) # Prevent division by zero
|
|
historical_vol = max(df["Log_Return"].tail(lookback_window).std(), 1e-6)
|
|
metrics["vol_ratio"] = recent_vol / historical_vol
|
|
|
|
# Trend strength using log prices
|
|
log_prices = np.log(df["Close"])
|
|
ma50 = log_prices.rolling(50, min_periods=1).mean()
|
|
ma200 = log_prices.rolling(200, min_periods=1).mean()
|
|
metrics["trend_strength"] = (ma50.iloc[-1] - ma200.iloc[-1]) / historical_vol
|
|
|
|
# Drawdown intensity
|
|
rolling_max = df["Close"].rolling(lookback_window, min_periods=1).max()
|
|
current_drawdown = df["Close"].iloc[-1] / rolling_max.iloc[-1] - 1
|
|
metrics["drawdown"] = abs(min(current_drawdown, 0))
|
|
|
|
return metrics
|
|
|
|
|
|
def get_projection_adjustments(days_forward, current_cycle_position, df):
|
|
"""
|
|
Enhanced projection adjustments with dynamic uncertainty scaling.
|
|
"""
|
|
adjustments = np.ones(days_forward)
|
|
|
|
# Get market condition metrics
|
|
conditions = calculate_market_conditions(df)
|
|
|
|
# Base uncertainty varies with market conditions
|
|
base_uncertainty = 0.016 # Standard rate
|
|
|
|
# Increase uncertainty if volatility is unusually high or low
|
|
vol_factor = 1 + 0.2 * abs(1 - conditions["vol_ratio"])
|
|
|
|
# Increase uncertainty during strong trends (both up and down)
|
|
trend_factor = 1 + 0.15 * abs(conditions["trend_strength"])
|
|
|
|
# Increase uncertainty during significant drawdowns
|
|
drawdown_factor = 1 + 0.25 * conditions["drawdown"]
|
|
|
|
# Combine factors with cycle position
|
|
if current_cycle_position > 0.75:
|
|
cycle_factor = 1.15 # Higher uncertainty late in cycle
|
|
else:
|
|
cycle_factor = 1.0
|
|
|
|
# Calculate time-varying uncertainty
|
|
for i in range(days_forward):
|
|
# Conservative growth with cycle and condition awareness
|
|
time_factor = min(
|
|
1
|
|
+ (i / 365)
|
|
* base_uncertainty
|
|
* vol_factor
|
|
* trend_factor
|
|
* drawdown_factor,
|
|
1.20,
|
|
)
|
|
|
|
# Update cycle position for this future point
|
|
cycle_position = (current_cycle_position + i / 1460) % 1
|
|
if cycle_position > 0.75:
|
|
local_cycle_factor = 1.15
|
|
else:
|
|
local_cycle_factor = 1.0
|
|
|
|
# Apply all factors including the initial cycle factor
|
|
adjustments[i] = time_factor * local_cycle_factor * cycle_factor
|
|
|
|
# Add minimum floor to prevent overconfidence
|
|
adjustments[i] = max(adjustments[i], 1.02 + (i / 365) * 0.01)
|
|
|
|
return adjustments
|
|
|
|
|
|
def calculate_confidence_intervals(simulated_paths, confidence_levels=[0.95, 0.68]):
|
|
"""
|
|
Calculate confidence intervals with dynamic quantile selection based on market conditions.
|
|
"""
|
|
results = {}
|
|
|
|
for level in confidence_levels:
|
|
# Calculate standard error of the median
|
|
median_std = np.std(
|
|
[np.median(simulated_paths[:, i]) for i in range(simulated_paths.shape[1])]
|
|
)
|
|
|
|
# Adjust quantiles based on estimation uncertainty
|
|
adjustment = min(0.1, median_std / np.median(simulated_paths)) # Cap adjustment
|
|
|
|
# Widen intervals slightly when uncertainty is high
|
|
effective_level = level + (1 - level) * adjustment
|
|
|
|
lower_percentile = (1 - effective_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 project_prices(
|
|
df, days_forward=365, simulations=1000, confidence_levels=[0.95, 0.68]
|
|
):
|
|
"""
|
|
Modified projection function incorporating enhanced uncertainty estimation.
|
|
"""
|
|
df = df.copy()
|
|
df["Log_Price"] = np.log(df["Close"])
|
|
df["Log_Return"] = df["Log_Price"].diff()
|
|
|
|
# Get halving dates and 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 with market awareness
|
|
projection_adjustments = get_projection_adjustments(
|
|
days_forward, cycle_position, df
|
|
)
|
|
|
|
# Run Monte Carlo simulation
|
|
np.random.seed(42) # Restored for reproducibility
|
|
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 with dynamic confidence intervals
|
|
results = pd.DataFrame(index=future_dates)
|
|
results["Median"] = np.percentile(simulated_paths, 50, axis=1)
|
|
results["Expected_Trend"] = last_price * np.exp(np.cumsum(drift))
|
|
|
|
# Calculate confidence intervals with dynamic adjustment
|
|
ci_results = calculate_confidence_intervals(simulated_paths, confidence_levels)
|
|
for key, values in ci_results.items():
|
|
results[key] = values
|
|
|
|
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 and spacing
|
|
fig = plt.figure(figsize=(15, 15))
|
|
|
|
# Use GridSpec for better control over subplot spacing
|
|
gs = plt.GridSpec(5, 1, height_ratios=[3, 1.5, 1.5, 1.5, 2], hspace=0.4)
|
|
|
|
# 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')})"
|
|
|
|
# Calculate full date range including projections
|
|
full_date_range = pd.date_range(plot_df["Date"].min(), projections.index.max())
|
|
|
|
# 1. Price history and projections (log scale)
|
|
ax1 = fig.add_subplot(gs[0])
|
|
|
|
# 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)
|
|
|
|
# Set x-axis limits to full range
|
|
ax1.set_xlim(full_date_range[0], full_date_range[-1])
|
|
ax1.tick_params(axis="x", rotation=45)
|
|
|
|
# 3. Rolling volatility
|
|
ax3 = fig.add_subplot(gs[1])
|
|
ax3.plot(
|
|
plot_df["Date"],
|
|
plot_df["Rolling_Volatility_30d"],
|
|
"r-",
|
|
label="30-Day Rolling Volatility",
|
|
)
|
|
|
|
# Add empty space to match price plot x-axis
|
|
ax3.set_xlim(full_date_range[0], full_date_range[-1])
|
|
|
|
# Add vertical line to mark start of projections
|
|
ax3.axvline(plot_df["Date"].max(), color="gray", linestyle="--", alpha=0.5)
|
|
ax3.text(
|
|
plot_df["Date"].max(),
|
|
ax3.get_ylim()[1],
|
|
"Projection Start",
|
|
rotation=90,
|
|
va="top",
|
|
ha="right",
|
|
alpha=0.7,
|
|
)
|
|
|
|
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()
|
|
ax3.tick_params(axis="x", rotation=45)
|
|
|
|
# 4. Returns distribution
|
|
ax4 = fig.add_subplot(gs[2])
|
|
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 = fig.add_subplot(gs[3:]) # Use last two grid spaces for larger plot
|
|
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)
|
|
|
|
# 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"
|
|
|
|
# Use tight_layout with adjusted parameters
|
|
plt.tight_layout(pad=2.0)
|
|
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")
|
|
_, 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
|
|
_ = create_plots(df, start=start, project_days=365 * 4)
|
|
|
|
|
|
def run_projections(df):
|
|
# Create main projection
|
|
projection_starts = [
|
|
"2011-01-01",
|
|
"2012-01-01",
|
|
"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=2, 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("2011-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("\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")
|