Performative scenario optimization solutions converge to classical stochastic optimization solutions in the limit where the decision-influence on the data-generating process vanishes, analogous to how entropic optimal transport converges to classical optimal transport as regularization approaches zero.

For a performative scenario optimization problem parameterized by a decision variable θ ∈ Θ ⊆ ℝ^d, let D(θ) denote the decision-dependent data distribution and D₀ denote the baseline (decision-independent) distribution. Define the influence parameter ε ≥ 0 such that D_ε(θ) = (1-ε)·D(θ) + ε·D₀ interpolates between full performativity (ε=0) and no performativity (ε=1). The hypothesis states: as ε → 1 (equivalently, as the Wasserstein distance W₂(D_ε(θ), D₀) → 0 uniformly in θ), the performative stable solution θ_PS(ε) converges in norm to the classical stochastic optimization solution θSO = argmin{θ} E_{D₀}[ℓ(θ, Z)], with convergence rate ‖θ_PS(ε) - θ_SO‖ = O((1-ε)^α) for some α > 0. This is falsifiable: if the limit does not exist, if convergence is non-monotone, or if the rate is slower than any polynomial, the hypothesis is refuted.

PHASE 1 — Analytical Verification (Days 1–14): Derive closed-form solutions for 3 tractable problem classes: (A) quadratic loss with Gaussian performative distribution, (B) linear regression with mean-shift performativity, (C) portfolio optimization with covariance-shift performativity. Verify convergence analytically and compute exact rates.

PHASE 2 — Numerical Convergence Study (Days 15–35): Implement gradient-based repeated risk minimization (RRM) and performative gradient descent across ε ∈ {0.0, 0.1, 0.2, ..., 0.9, 0.95, 0.99, 0.999, 0.9999}. Run 5 problem families × 10 random seeds × 20 ε values = 1000 experiments. Measure ‖θ_PS(ε) - θ_SO‖₂ at convergence.

PHASE 3 — Rate Estimation (Days 36–50): Fit power-law model ‖error‖ = C·(1-ε)^α via log-log regression. Bootstrap confidence intervals (n=1000 resamples). Compare α estimates across problem families.

PHASE 4 — Stress Testing (Days 51–65): Test near-boundary conditions: β approaching μ/L, non-smooth losses (Huber), heavy-tailed distributions (Student-t with ν=3), and high-dimensional θ (d ∈ {10, 100, 1000}).

PHASE 5 — Analogy Formalization (Days 66–80): Construct explicit mathematical mapping between performativity parameter ε and entropic OT regularization λ. Verify whether convergence rates match under this mapping.

# Core Implementation Architecture

import jax
import jax.numpy as jnp
from jax import grad, jit, vmap
import optax
from dataclasses import dataclass
from typing import Callable, Tuple
import numpy as np

@dataclass
class PerformativeProblem:
    """Parameterized performative optimization problem."""
    loss_fn: Callable          # ℓ(θ, z) -> scalar
    dist_map: Callable         # θ -> distribution parameters
    base_dist_params: dict     # D₀ parameters
    theta_dim: int
    epsilon: float = 0.0       # performativity parameter (0=full, 1=classical)
    
    def interpolated_dist_params(self, theta):
        """D_ε(θ) = (1-ε)·D(θ) + ε·D₀"""
        perf_params = self.dist_map(theta)
        # Interpolate distribution parameters
        return jax.tree_util.tree_map(
            lambda p, b: (1 - self.epsilon) * p + self.epsilon * b,
            perf_params, self.base_dist_params
        )
    
    def expected_loss(self, theta, n_samples=1000, key=jax.random.PRNGKey(0)):
        """Monte Carlo estimate of E_{D_ε(θ)}[ℓ(θ, Z)]"""
        dist_params = self.interpolated_dist_params(theta)
        samples = sample_distribution(dist_params, n_samples, key)
        return jnp.mean(vmap(lambda z: self.loss_fn(theta, z))(samples))

def repeated_risk_minimization(problem: PerformativeProblem, 
                                theta_init: jnp.ndarray,
                                max_iter: int = 10000,
                                tol: float = 1e-8,
                                lr: float = 0.01) -> Tuple[jnp.ndarray, int]:
    """
    RRM: θ_{t+1} = argmin_θ E_{D(θ_t)}[ℓ(θ, Z)]
    Returns: (theta_star, n_iterations)
    """
    optimizer = optax.adam(lr)
    theta = theta_init.copy()
    opt_state = optimizer.init(theta)
    
    for t in range(max_iter):
        theta_prev = theta.copy()
        # Fix distribution at current theta, minimize over theta
        fixed_dist_params = problem.interpolated_dist_params(theta)
        
        def fixed_loss(th):
            samples = sample_distribution(fixed_dist_params, n_samples=500)
            return jnp.mean(vmap(lambda z: problem.loss_fn(th, z))(samples))
        
        # Inner optimization loop
        theta = inner_minimize(fixed_loss, theta, n_steps=100, lr=lr)
        
        if jnp.norm(theta - theta_prev) < tol:
            return theta, t
    
    return theta, max_iter  # Did not converge

def convergence_experiment(problem_family: str,
                           epsilon_values: list,
                           n_seeds: int = 10) -> dict:
    """
    Main experiment: measure ||θ*_PS(ε) - θ*_SO||₂ across ε values.
    """
    results = {'epsilon': [], 'error': [], 'seed': [], 'family': []}
    
    # Solve classical problem once
    classical_problem = create_problem(problem_family, epsilon=1.0)
    theta_SO = solve_classical(classical_problem)
    
    for eps in epsilon_values:
        for seed in range(n_seeds):
            key = jax.random.PRNGKey(seed)
            theta_init = jax.random.normal(key, (classical_problem.theta_dim,))
            
            perf_problem = create_problem(problem_family, epsilon=eps)
            theta_PS, n_iter = repeated_risk_minimization(
                perf_problem, theta_init
            )
            
            error = float(jnp.norm(theta_PS - theta_SO))
            results['epsilon'].append(eps)
            results['error'].append(error)
            results['seed'].append(seed)
            results['family'].append(problem_family)
    
    return results

def estimate_convergence_rate(results: dict) -> dict:
    """
    Fit power law: log(error) = log(C) + α·log(1-ε)
    Returns: {'alpha': float, 'alpha_ci': (float, float), 'r_squared': float}
    """
    eps = np.array(results['epsilon'])
    errors = np.array(results['error'])
    
    # Filter to high-ε regime
    mask = eps >= 0.9
    log_x = np.log(1 - eps[mask])
    log_y = np.log(errors[mask] + 1e-12)
    
    # OLS fit
    A = np.column_stack([np.ones_like(log_x), log_x])
    coeffs, residuals, _, _ = np.linalg.lstsq(A, log_y, rcond=None)
    log_C, alpha = coeffs
    
    # Bootstrap CI
    alpha_bootstrap = []
    for _ in range(1000):
        idx = np.random.choice(len(log_x), len(log_x), replace=True)
        c, _, _, _ = np.linalg.lstsq(A[idx], log_y[idx], rcond=None)
        alpha_bootstrap.append(c[1])
    
    r_squared = 1 - np.var(log_y - A @ coeffs) / np.var(log_y)
    
    return {
        'alpha': alpha,
        'alpha_ci': (np.percentile(alpha_bootstrap, 2.5),
                     np.percentile(alpha_bootstrap, 97.5)),
        'r_squared': r_squared,
        'C': np.exp(log_C)
    }

# ANALYTICAL VERIFICATION: Quadratic case
# ℓ(θ,z) = ½||θ-z||², D(θ) = N(Aθ+b, Σ)
# θ*_PS = (I-A)⁻¹b  [requires ||A|| < 1]
# D_ε(θ) = N(ε·(Aθ+b) + (1-ε)·b, Σ) [simplified interpolation]
# Wait: D_ε(θ) = N((1-ε)·(Aθ+b) + ε·b, Σ) = N((1-ε)Aθ+b, Σ)
# θ*_PS(ε) = (I-(1-ε)A)⁻¹b
# As ε→1: θ*_PS(ε) → (I-0)⁻¹b = b = θ*_SO ✓
# Rate: ||θ*_PS(ε) - θ*_SO|| = ||(I-(1-ε)A)⁻¹b - b||
#      ≈ ||(1-ε)A·b|| = O(1-ε) for small (1-ε)

def analytical_quadratic_solution(A, b, epsilon):
    """Closed-form θ*_PS(ε) for quadratic performative problem."""
    I = jnp.eye(A.shape[0])
    return jnp.linalg.solve(I - (1 - epsilon) * A, b)

# MAIN EXECUTION
if __name__ == "__main__":
    epsilon_grid = [0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 
                    0.95, 0.99, 0.999, 0.9999]
    problem_families = ['quadratic', 'linear_regression', 
                        'portfolio', 'strategic_classification',
                        'credit_scoring']
    
    all_results = {}
    for family in problem_families:
        print(f"Running experiments for {family}...")
        results = convergence_experiment(family, epsilon_grid, n_seeds=10)
        rate_info = estimate_convergence_rate(results)
        all_results[family] = {'data': results, 'rate': rate_info}
        print(f"  α = {rate_info['alpha']:.3f} "
              f"[{rate_info['alpha_ci'][0]:.3f}, {rate_info['alpha_ci'][1]:.3f}], "
              f"R² = {rate_info['r_squared']:.4f}")
    
    # Check success criteria
    convergence_check = all(
        min(r['data']['error'][-10:]) < 1e-3  # last 10 = highest ε
        for r in all_results.values()
    )
    print(f"\nConvergence criterion met: {convergence_check}")