4186.py

import numpy as np
import matplotlib.pyplot as plt


def simulate_market_shares(n_firms, n_periods, volatility):
    shares = np.random.dirichlet(np.ones(n_firms), n_periods)
    # Apply some smoothing to simulate persistence
    for t in range(1, n_periods):
        shares[t] = shares[t - 1] + volatility * (shares[t] - shares[t - 1])
    return shares


def calculate_expected_profits(shares, phi_star, Pi, Pi_star, delta):
    n_periods = shares.shape[0]
    expected_profits = np.zeros((n_periods, shares.shape[1]))

    for t in range(n_periods - 1, -1, -1):
        p_above_threshold = np.mean(np.all(shares[t:] >= phi_star, axis=1))
        for i in range(shares.shape[1]):
            if t == n_periods - 1:
                expected_profits[t, i] = shares[t, i] * (
                    p_above_threshold * Pi + (1 - p_above_threshold) * Pi_star
                )
            else:
                expected_profits[t, i] = (
                    shares[t, i]
                    * (p_above_threshold * Pi + (1 - p_above_threshold) * Pi_star)
                    + delta * expected_profits[t + 1, i]
                )

    return expected_profits


def check_equilibrium_conditions(shares, expected_profits, phi_star, Pi, Pi_star):
    conditions_met = np.ones(shares.shape, dtype=bool)

    for t in range(shares.shape[0]):
        for i in range(shares.shape[1]):
            if np.all(shares[t] >= phi_star):
                conditions_met[t, i] = expected_profits[t, i] >= (1 - shares[t, i]) * Pi
            else:
                conditions_met[t, i] = (
                    expected_profits[t, i] >= (1 - shares[t, i]) * Pi_star
                )

    return conditions_met


# Set parameters
n_firms = 3
n_periods = 100
volatility = 0.1
phi_star = 0.2
Pi = 100
Pi_star = 80
delta = 0.95

# Run simulation
shares = simulate_market_shares(n_firms, n_periods, volatility)
expected_profits = calculate_expected_profits(shares, phi_star, Pi, Pi_star, delta)
equilibrium_conditions = check_equilibrium_conditions(
    shares, expected_profits, phi_star, Pi, Pi_star
)

# Plot results
plt.figure(figsize=(12, 8))
plt.subplot(211)
for i in range(n_firms):
    plt.plot(shares[:, i], label=f"Firm {i+1}")
plt.axhline(y=phi_star, color="r", linestyle="--", label="ϕ*")
plt.title("Market Shares Over Time")
plt.legend()

plt.subplot(212)
for i in range(n_firms):
    plt.plot(equilibrium_conditions[:, i], label=f"Firm {i+1}")
plt.title("Equilibrium Conditions Met")
plt.legend()

plt.tight_layout()
plt.show()

print(
    f"Percentage of periods where all firms meet equilibrium conditions: {np.mean(np.all(equilibrium_conditions, axis=1)) * 100:.2f}%"
)