Monte Carlo Analysis

import sys
sys.path.insert(0, '../')

from src.model import StadiumEconomicModel
from src.simulation import BeerPriceControlSimulator
import numpy as np
import pandas as pd
import plotly.graph_objects as go
import plotly.express as px
from scipy import stats

# Parameter distributions (from literature and reasonable bounds)
# CRITICAL: cross_price_elasticity is the most uncertain and influential parameter
param_distributions = {
    'cross_price_elasticity': {'dist': 'uniform', 'low': 0.0, 'high': 0.3},  # Key uncertainty!
    'drinker_share': {'dist': 'uniform', 'low': 0.30, 'high': 0.50},  # Lenk et al. 2010: ~40%
    'ticket_cost': {'dist': 'uniform', 'low': 2.0, 'high': 5.0},
    'beer_cost': {'dist': 'uniform', 'low': 4.0, 'high': 6.0},
    'experience_degradation_cost': {'dist': 'uniform', 'low': 40, 'high': 100},
    'crime_cost_per_beer': {'dist': 'uniform', 'low': 1.5, 'high': 3.5},
    'health_cost_per_beer': {'dist': 'uniform', 'low': 1.0, 'high': 2.0},
}

print("Parameter Ranges:")
print()
print("CRITICAL ASSUMPTIONS:")
print(f"  cross_price_elasticity: [{param_distributions['cross_price_elasticity']['low']}, {param_distributions['cross_price_elasticity']['high']}]")
print("    (0.0 = Leisten one-way, 0.1 = baseline, 0.3 = strong complementarity)")
print(f"  drinker_share: [{param_distributions['drinker_share']['low']}, {param_distributions['drinker_share']['high']}]")
print("    (Lenk et al. 2010 reports ~40%; we test 30-50% range)")
print()
print("Other parameters:")
for param, dist_info in param_distributions.items():
    if param not in ['cross_price_elasticity', 'drinker_share']:
        print(f"  {param}: [{dist_info['low']}, {dist_info['high']}]")

from src.model import ConsumerType

np.random.seed(42)  # Reproducibility
n_simulations = 1000

results = []

for i in range(n_simulations):
    # Sample parameters
    params = {}
    for param, dist_info in param_distributions.items():
        params[param] = np.random.uniform(dist_info['low'], dist_info['high'])
    
    # Create consumer types with sampled drinker_share
    drinker_share = params['drinker_share']
    custom_types = [
        ConsumerType(
            name="Non-Drinker",
            share=1.0 - drinker_share,
            alpha_beer=1.0,
            alpha_experience=3.0,
            income=200.0,
        ),
        ConsumerType(
            name="Drinker",
            share=drinker_share,
            alpha_beer=43.75,
            alpha_experience=2.5,
            income=200.0,
        ),
    ]
    
    # Create model with sampled parameters
    try:
        model = StadiumEconomicModel(
            consumer_types=custom_types,
            cross_price_elasticity=params['cross_price_elasticity'],
            ticket_cost=params['ticket_cost'],
            beer_cost=params['beer_cost'],
            experience_degradation_cost=params['experience_degradation_cost']
        )
        
        # Temporarily set external costs for welfare calculation
        model.external_costs['crime'] = params['crime_cost_per_beer']
        model.external_costs['health'] = params['health_cost_per_beer']
        
        # Baseline: beer at $12.50
        t_current, b_current, r_current = model.optimal_pricing(beer_price_control=12.5)
        
        # $7 ceiling
        t_ceiling, b_ceiling, r_ceiling = model.optimal_pricing(beer_price_control=7.0)
        
        # Calculate welfare (using model's internal external costs)
        welfare_current = model.social_welfare(t_current, b_current)
        welfare_ceiling = model.social_welfare(t_ceiling, b_ceiling)
        
        # Store results
        results.append({
            **params,
            'ticket_current': t_current,
            'ticket_ceiling': t_ceiling,
            'ticket_change': t_ceiling - t_current,
            'ticket_pct_change': (t_ceiling / t_current - 1) * 100,
            'attendance_current': r_current['attendance'],
            'attendance_ceiling': r_ceiling['attendance'],
            'attendance_change': r_ceiling['attendance'] - r_current['attendance'],
            'beers_current': r_current['total_beers'],
            'beers_ceiling': r_ceiling['total_beers'],
            'beers_change': r_ceiling['total_beers'] - r_current['total_beers'],
            'beers_per_fan_current': r_current['beers_per_fan'],
            'beers_per_fan_ceiling': r_ceiling['beers_per_fan'],
            'profit_current': r_current['profit'],
            'profit_ceiling': r_ceiling['profit'],
            'profit_change': r_ceiling['profit'] - r_current['profit'],
            'cs_current': welfare_current['consumer_surplus'],
            'cs_ceiling': welfare_ceiling['consumer_surplus'],
            'cs_change': welfare_ceiling['consumer_surplus'] - welfare_current['consumer_surplus'],
            'ps_current': welfare_current['producer_surplus'],
            'ps_ceiling': welfare_ceiling['producer_surplus'],
            'ps_change': welfare_ceiling['producer_surplus'] - welfare_current['producer_surplus'],
            'sw_current': welfare_current['social_welfare'],
            'sw_ceiling': welfare_ceiling['social_welfare'],
            'sw_change': welfare_ceiling['social_welfare'] - welfare_current['social_welfare'],
            'externality_current': welfare_current['externality_cost'],
            'externality_ceiling': welfare_ceiling['externality_cost'],
            'externality_change': welfare_ceiling['externality_cost'] - welfare_current['externality_cost'],
        })
    except Exception as e:
        # Skip failed optimizations
        continue

df = pd.DataFrame(results)
print(f"Successful simulations: {len(df)}/{n_simulations}")
print(f"Failed: {n_simulations - len(df)}")

# Summary statistics
print("TICKET PRICE CHANGE ($7 Beer Ceiling):")
print(f"  Mean:   ${df['ticket_change'].mean():.2f}")
print(f"  Median: ${df['ticket_change'].median():.2f}")
print(f"  Std:    ${df['ticket_change'].std():.2f}")
print(f"  5th %:  ${df['ticket_change'].quantile(0.05):.2f}")
print(f"  95th %: ${df['ticket_change'].quantile(0.95):.2f}")
print()

print("WELFARE CHANGES:")
print(f"  Consumer surplus: ${df['cs_change'].mean()/1e6:.2f}M ± ${df['cs_change'].std()/1e6:.2f}M")
print(f"  Producer surplus: ${df['ps_change'].mean()/1e6:.2f}M ± ${df['ps_change'].std()/1e6:.2f}M")
print(f"  Social welfare:   ${df['sw_change'].mean()/1e6:.2f}M ± ${df['sw_change'].std()/1e6:.2f}M")
print()

# Probability statements
prob_ticket_rise = (df['ticket_change'] > 0).mean()
prob_profit_fall = (df['profit_change'] < 0).mean()
prob_welfare_rise = (df['sw_change'] > 0).mean()

print("PROBABILITY:")
print(f"  Tickets rise: {prob_ticket_rise:.1%}")
print(f"  Stadium profit falls: {prob_profit_fall:.1%}")
print(f"  Social welfare rises: {prob_welfare_rise:.1%}")

# Histogram of ticket changes
fig1 = go.Figure()
fig1.add_trace(go.Histogram(
    x=df['ticket_change'],
    nbinsx=50,
    name='Ticket Change',
    marker_color='#003087'
))

fig1.add_vline(
    x=df['ticket_change'].median(),
    line_dash="dash",
    annotation_text=f"Median: ${df['ticket_change'].median():.2f}"
)

fig1.update_layout(
    title='Distribution of Ticket Price Changes (Monte Carlo)',
    xaxis_title='Ticket Price Change ($)',
    yaxis_title='Frequency',
    height=400
)

fig1.show()

# Welfare components
fig2 = go.Figure()

fig2.add_trace(go.Box(
    y=df['cs_change'] / 1e6,
    name='Consumer Surplus',
    marker_color='lightblue'
))

fig2.add_trace(go.Box(
    y=df['ps_change'] / 1e6,
    name='Producer Surplus',
    marker_color='lightgreen'
))

fig2.add_trace(go.Box(
    y=df['sw_change'] / 1e6,
    name='Social Welfare',
    marker_color='orange'
))

fig2.update_layout(
    title='Welfare Changes from $7 Beer Ceiling (Distribution)',
    yaxis_title='Change ($ Millions per game)',
    height=500
)

fig2.show()

# Correlation analysis - cross_price_elasticity and drinker_share are KEY drivers
correlations = df[['cross_price_elasticity', 'drinker_share', 'ticket_cost',
                    'beer_cost', 'experience_degradation_cost',
                    'crime_cost_per_beer', 'health_cost_per_beer']].corrwith(df['ticket_change'])

print("Correlation with Ticket Price Change:")
print()
for param in correlations.index:
    corr = correlations[param]
    marker = "**" if abs(corr) > 0.3 else ""
    print(f"  {marker}{param:<30} {corr:+.3f}{marker}")
print()

print("Interpretation:")
print("  Positive correlation: Higher parameter → larger ticket increase")
print("  Negative correlation: Higher parameter → smaller ticket increase")
print()
print("** indicates strong influence (|r| > 0.3)")
print()

# Show how cross_price_elasticity affects results
print("\nCROSS-PRICE ELASTICITY IMPACT:")
low_cross = df[df['cross_price_elasticity'] < 0.05]
mid_cross = df[(df['cross_price_elasticity'] >= 0.05) & (df['cross_price_elasticity'] < 0.15)]
high_cross = df[df['cross_price_elasticity'] >= 0.15]

print(f"  Low (0-0.05, Leisten-like):  Ticket change = ${low_cross['ticket_change'].mean():+.2f}")
print(f"  Medium (0.05-0.15, baseline): Ticket change = ${mid_cross['ticket_change'].mean():+.2f}")
print(f"  High (0.15-0.30):             Ticket change = ${high_cross['ticket_change'].mean():+.2f}")

# Show how drinker_share affects results
print("\nDRINKER SHARE IMPACT:")
low_drinker = df[df['drinker_share'] < 0.35]
mid_drinker = df[(df['drinker_share'] >= 0.35) & (df['drinker_share'] < 0.45)]
high_drinker = df[df['drinker_share'] >= 0.45]

print(f"  Low (30-35%):    Beers/fan change = {low_drinker['beers_per_fan_ceiling'].mean() - low_drinker['beers_per_fan_current'].mean():+.2f}")
print(f"  Medium (35-45%): Beers/fan change = {mid_drinker['beers_per_fan_ceiling'].mean() - mid_drinker['beers_per_fan_current'].mean():+.2f}")
print(f"  High (45-50%):   Beers/fan change = {high_drinker['beers_per_fan_ceiling'].mean() - high_drinker['beers_per_fan_current'].mean():+.2f}")

print("WELFARE DECOMPOSITION (Mean across simulations):")
print()
print("Current ($12.50 beer):")
print(f"  Consumer Surplus:  ${df['cs_current'].mean()/1e6:.2f}M")
print(f"  Producer Surplus:  ${df['ps_current'].mean()/1e6:.2f}M")
print(f"  Externality Cost:  ${df['externality_current'].mean()/1e3:.0f}k")
print(f"  Social Welfare:    ${df['sw_current'].mean()/1e6:.2f}M")
print()

# Calculate shares
total_surplus_current = df['cs_current'].mean() + df['ps_current'].mean()
cs_share = df['cs_current'].mean() / total_surplus_current * 100
ps_share = df['ps_current'].mean() / total_surplus_current * 100

print(f"  Consumer share of surplus: {cs_share:.1f}%")
print(f"  Producer share of surplus: {ps_share:.1f}%")
print()

print("With $7 Ceiling:")
print(f"  Consumer Surplus:  ${df['cs_ceiling'].mean()/1e6:.2f}M")
print(f"  Producer Surplus:  ${df['ps_ceiling'].mean()/1e6:.2f}M")
print(f"  Externality Cost:  ${df['externality_ceiling'].mean()/1e3:.0f}k")
print(f"  Social Welfare:    ${df['sw_ceiling'].mean()/1e6:.2f}M")
print()

print("Changes:")
print(f"  Consumer Surplus:  ${df['cs_change'].mean()/1e6:+.2f}M ({df['cs_change'].mean()/df['cs_current'].mean()*100:+.1f}%)")
print(f"  Producer Surplus:  ${df['ps_change'].mean()/1e6:+.2f}M ({df['ps_change'].mean()/df['ps_current'].mean()*100:+.1f}%)")
print(f"  Social Welfare:    ${df['sw_change'].mean()/1e6:+.2f}M ({df['sw_change'].mean()/df['sw_current'].mean()*100:+.1f}%)")

# Calculate probability of each group being better off
prob_consumers_win = (df['cs_change'] > 0).mean()
prob_stadium_loses = (df['ps_change'] < 0).mean()
prob_society_wins = (df['sw_change'] > 0).mean()

print("PROBABILITY OF WELFARE GAINS FROM $7 CEILING:")
print()
print(f"  Consumers gain:       {prob_consumers_win:.1%}")
print(f"  Stadium loses:        {prob_stadium_loses:.1%}")
print(f"  Net social gain:      {prob_society_wins:.1%}")
print()

# When does society gain?
winners = df[df['sw_change'] > 0]
losers = df[df['sw_change'] < 0]

if len(winners) > 0:
    print(f"When society GAINS ({len(winners)} scenarios):")
    print(f"  Average CS gain:  ${winners['cs_change'].mean()/1e6:.2f}M")
    print(f"  Average PS loss:  ${winners['ps_change'].mean()/1e6:.2f}M")
    print(f"  Net social gain:  ${winners['sw_change'].mean()/1e6:.2f}M")
    print()

if len(losers) > 0:
    print(f"When society LOSES ({len(losers)} scenarios):")
    print(f"  Average CS change: ${losers['cs_change'].mean()/1e6:.2f}M")
    print(f"  Average PS loss:   ${losers['ps_change'].mean()/1e6:.2f}M")
    print(f"  Net social loss:   ${losers['sw_change'].mean()/1e6:.2f}M")

fig3 = px.scatter(
    df,
    x='ticket_change',
    y='sw_change',
    color='experience_degradation_cost',
    title='Ticket Price Response vs Social Welfare Change',
    labels={
        'ticket_change': 'Ticket Price Increase ($)',
        'sw_change': 'Social Welfare Change ($)',
        'experience_degradation_cost': 'Internalized Cost (α)'
    },
    height=500
)

fig3.add_hline(y=0, line_dash="dash", line_color="gray")
fig3.add_vline(x=0, line_dash="dash", line_color="gray")

fig3.show()