chrishwiggins · September 27, 2024 14:55 · Sep 27, 2024
diff --git a/4186.py b/4186.py
@@ -0,0 +1,85 @@
+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}%"
+)
No results found