import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm, multivariate_normal
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from qiskit import QuantumCircuit, Aer, execute
from qiskit.circuit.library import RealAmplitudes, ZZFeatureMap
from qiskit.algorithms import VQE, QAOA, AmplitudeEstimation
from qiskit.algorithms.optimizers import SPSA, COBYLA
from qiskit.quantum_info import state_fidelity
from qiskit.opflow import PauliSumOp, StateFn, CircuitSampler
import pennylane as qml

class ClassicalStressTesting:
    """Classical stress testing implementation"""
    
    def __init__(self):
        self.risk_free_rate = 0.02
        self.recovery_rate = 0.4
        
    def generate_stress_scenarios(self, n_scenarios=1000):
        """
        Generate stress testing scenarios
        """
        np.random.seed(42)
        
        # Base economic indicators
        base_gdp_growth = 0.03
        base_unemployment = 0.05
        base_interest_rate = 0.03
        base_inflation = 0.02
        
        scenarios = []
        
        for i in range(n_scenarios):
            # Generate stress levels (0 = normal, 1 = severe stress)
            stress_level = np.random.beta(2, 8)  # Most scenarios are normal to moderate
            
            # Apply stress to economic indicators
            gdp_growth = base_gdp_growth * (1 - stress_level * 2)  # Can go negative
            unemployment = base_unemployment * (1 + stress_level * 3)
            interest_rate = base_interest_rate * (1 + stress_level * 2)
            inflation = base_inflation * (1 + stress_level * 2)
            
            # Credit-specific stress factors
            default_prob_multiplier = 1 + stress_level * 4  # Up to 5x increase
            recovery_rate_multiplier = 1 - stress_level * 0.5  # Down to 50%
            correlation_multiplier = 1 + stress_level * 2  # Up to 3x increase
            
            scenario = {
                'scenario_id': i,
                'stress_level': stress_level,
                'gdp_growth': gdp_growth,
                'unemployment': unemployment,
                'interest_rate': interest_rate,
                'inflation': inflation,
                'default_prob_multiplier': default_prob_multiplier,
                'recovery_rate_multiplier': recovery_rate_multiplier,
                'correlation_multiplier': correlation_multiplier
            }
            
            scenarios.append(scenario)
        
        return scenarios
    
    def calculate_portfolio_loss(self, portfolio, scenario):
        """
        Calculate portfolio loss under stress scenario
        """
        total_loss = 0
        
        for asset in portfolio:
            # Base default probability
            base_default_prob = asset['default_probability']
            
            # Apply stress multiplier
            stressed_default_prob = base_default_prob * scenario['default_prob_multiplier']
            stressed_default_prob = min(stressed_default_prob, 1.0)  # Cap at 100%
            
            # Base recovery rate
            base_recovery_rate = asset['recovery_rate']
            
            # Apply stress multiplier
            stressed_recovery_rate = base_recovery_rate * scenario['recovery_rate_multiplier']
            stressed_recovery_rate = max(stressed_recovery_rate, 0.1)  # Floor at 10%
            
            # Calculate loss given default
            lgd = 1 - stressed_recovery_rate
            
            # Portfolio loss for this asset
            asset_loss = asset['exposure'] * stressed_default_prob * lgd
            
            total_loss += asset_loss
        
        return total_loss
    
    def monte_carlo_stress_test(self, portfolio, n_scenarios=1000):
        """
        Perform Monte Carlo stress testing
        """
        scenarios = self.generate_stress_scenarios(n_scenarios)
        
        losses = []
        for scenario in scenarios:
            loss = self.calculate_portfolio_loss(portfolio, scenario)
            losses.append(loss)
        
        return np.array(losses), scenarios
    
    def calculate_stress_measures(self, losses, confidence_level=0.95):
        """
        Calculate stress testing risk measures
        """
        # Sort losses in descending order
        sorted_losses = np.sort(losses)[::-1]
        
        # VaR
        var_index = int((1 - confidence_level) * len(sorted_losses))
        var = sorted_losses[var_index]
        
        # CVaR (Expected Shortfall)
        cvar = np.mean(sorted_losses[:var_index + 1])
        
        # Maximum loss
        max_loss = np.max(losses)
        
        # Mean loss
        mean_loss = np.mean(losses)
        
        # Loss volatility
        loss_volatility = np.std(losses)
        
        return {
            'var': var,
            'cvar': cvar,
            'max_loss': max_loss,
            'mean_loss': mean_loss,
            'loss_volatility': loss_volatility,
            'confidence_level': confidence_level
        }

class QuantumStressTesting:
    """Quantum stress testing implementation"""
    
    def __init__(self, num_qubits=8):
        self.num_qubits = num_qubits
        self.backend = Aer.get_backend('qasm_simulator')
        self.optimizer = SPSA(maxiter=100)
        
    def create_stress_scenario_circuit(self, scenario_params):
        """
        Create quantum circuit for stress scenario encoding
        """
        # Encode stress parameters
        feature_map = ZZFeatureMap(feature_dimension=len(scenario_params), reps=2)
        
        # Ansatz for scenario generation
        ansatz = RealAmplitudes(num_qubits=self.num_qubits, reps=3)
        
        # Combine circuits
        circuit = feature_map.compose(ansatz)
        
        return circuit
    
    def create_loss_hamiltonian(self, portfolio, scenario_params):
        """
        Create quantum Hamiltonian for loss calculation
        """
        # Encode portfolio and scenario information
        hamiltonian_terms = []
        
        # Portfolio exposure term
        total_exposure = sum(asset['exposure'] for asset in portfolio)
        pauli_z = PauliSumOp.from_list([('Z', 1.0)])
        hamiltonian_terms.append((total_exposure, pauli_z))
        
        # Default probability term
        avg_default_prob = np.mean([asset['default_probability'] for asset in portfolio])
        stress_multiplier = scenario_params.get('default_prob_multiplier', 1.0)
        pauli_x = PauliSumOp.from_list([('X', 1.0)])
        hamiltonian_terms.append((avg_default_prob * stress_multiplier, pauli_x))
        
        # Recovery rate term
        avg_recovery_rate = np.mean([asset['recovery_rate'] for asset in portfolio])
        recovery_multiplier = scenario_params.get('recovery_rate_multiplier', 1.0)
        pauli_y = PauliSumOp.from_list([('Y', 1.0)])
        hamiltonian_terms.append((avg_recovery_rate * recovery_multiplier, pauli_y))
        
        return sum(term[0] * term[1] for term in hamiltonian_terms)
    
    def quantum_stress_scenario_generation(self, n_scenarios=100):
        """
        Generate stress scenarios using quantum methods
        """
        scenarios = []
        
        for i in range(n_scenarios):
            # Generate random scenario parameters
            stress_level = np.random.beta(2, 8)
            
            scenario_params = {
                'stress_level': stress_level,
                'default_prob_multiplier': 1 + stress_level * 4,
                'recovery_rate_multiplier': 1 - stress_level * 0.5,
                'correlation_multiplier': 1 + stress_level * 2,
                'gdp_growth': 0.03 * (1 - stress_level * 2),
                'unemployment': 0.05 * (1 + stress_level * 3),
                'interest_rate': 0.03 * (1 + stress_level * 2)
            }
            
            # Create quantum circuit for scenario
            circuit = self.create_stress_scenario_circuit(scenario_params)
            
            # Execute circuit to get quantum scenario
            job = execute(circuit, self.backend, shots=1000)
            result = job.result()
            counts = result.get_counts()
            
            # Extract scenario from quantum measurement
            quantum_scenario = self._extract_scenario_from_counts(counts, scenario_params)
            
            scenarios.append(quantum_scenario)
        
        return scenarios
    
    def _extract_scenario_from_counts(self, counts, base_params):
        """
        Extract scenario parameters from quantum measurement counts
        """
        # Simplified extraction - in practice, use more sophisticated methods
        total_shots = sum(counts.values())
        
        # Calculate quantum adjustments
        quantum_adjustment = 0.0
        for bitstring, count in counts.items():
            probability = count / total_shots
            parity = sum(int(bit) for bit in bitstring) % 2
            quantum_adjustment += probability * (1 if parity == 0 else -1)
        
        # Apply quantum adjustment to base parameters
        adjusted_params = base_params.copy()
        adjustment_factor = 1 + quantum_adjustment * 0.1  # 10% adjustment
        
        adjusted_params['default_prob_multiplier'] *= adjustment_factor
        adjusted_params['recovery_rate_multiplier'] /= adjustment_factor
        adjusted_params['correlation_multiplier'] *= adjustment_factor
        
        return adjusted_params
    
    def quantum_portfolio_loss_calculation(self, portfolio, scenario_params):
        """
        Calculate portfolio loss using quantum methods
        """
        # Create quantum circuit
        circuit = self.create_stress_scenario_circuit(scenario_params)
        
        # Create loss Hamiltonian
        hamiltonian = self.create_loss_hamiltonian(portfolio, scenario_params)
        
        # Calculate quantum expectation
        expectation = self._calculate_quantum_expectation(circuit, hamiltonian)
        
        # Convert to portfolio loss
        loss = self._convert_expectation_to_loss(expectation, portfolio, scenario_params)
        
        return loss
    
    def _calculate_quantum_expectation(self, circuit, hamiltonian):
        """
        Calculate quantum expectation value
        """
        # Execute circuit
        job = execute(circuit, self.backend, shots=1000)
        result = job.result()
        counts = result.get_counts()
        
        # Calculate expectation value
        expectation = 0.0
        total_shots = sum(counts.values())
        
        for bitstring, count in counts.items():
            probability = count / total_shots
            
            # Simplified expectation calculation
            parity = sum(int(bit) for bit in bitstring) % 2
            expectation += probability * (1 if parity == 0 else -1)
        
        return expectation
    
    def _convert_expectation_to_loss(self, expectation, portfolio, scenario_params):
        """
        Convert quantum expectation to portfolio loss
        """
        # Calculate base loss
        base_loss = 0
        for asset in portfolio:
            stressed_default_prob = asset['default_probability'] * scenario_params['default_prob_multiplier']
            stressed_recovery_rate = asset['recovery_rate'] * scenario_params['recovery_rate_multiplier']
            lgd = 1 - stressed_recovery_rate
            asset_loss = asset['exposure'] * stressed_default_prob * lgd
            base_loss += asset_loss
        
        # Apply quantum adjustment
        quantum_adjustment = 1 + expectation * 0.2  # 20% adjustment factor
        adjusted_loss = base_loss * quantum_adjustment
        
        return adjusted_loss
    
    def quantum_stress_test(self, portfolio, n_scenarios=100):
        """
        Perform quantum stress testing
        """
        # Generate quantum stress scenarios
        scenarios = self.quantum_stress_scenario_generation(n_scenarios)
        
        losses = []
        for scenario in scenarios:
            loss = self.quantum_portfolio_loss_calculation(portfolio, scenario)
            losses.append(loss)
        
        return np.array(losses), scenarios
    
    def quantum_stress_measures(self, losses, confidence_level=0.95):
        """
        Calculate quantum stress testing risk measures
        """
        # Sort losses in descending order
        sorted_losses = np.sort(losses)[::-1]
        
        # Quantum VaR
        var_index = int((1 - confidence_level) * len(sorted_losses))
        quantum_var = sorted_losses[var_index]
        
        # Quantum CVaR
        quantum_cvar = np.mean(sorted_losses[:var_index + 1])
        
        # Other measures
        max_loss = np.max(losses)
        mean_loss = np.mean(losses)
        loss_volatility = np.std(losses)
        
        return {
            'quantum_var': quantum_var,
            'quantum_cvar': quantum_cvar,
            'max_loss': max_loss,
            'mean_loss': mean_loss,
            'loss_volatility': loss_volatility,
            'confidence_level': confidence_level
        }

def generate_test_portfolio(n_assets=50):
    """
    Generate test portfolio for stress testing
    """
    np.random.seed(42)
    
    portfolio = []
    
    for i in range(n_assets):
        asset = {
            'asset_id': f'Asset_{i}',
            'exposure': np.random.uniform(100000, 1000000),  # $100K to $1M
            'default_probability': np.random.beta(2, 98),  # 0-1% default probability
            'recovery_rate': np.random.beta(4, 6),  # 40% average recovery
            'rating': np.random.choice(['AAA', 'AA', 'A', 'BBB', 'BB', 'B']),
            'sector': np.random.choice(['Financial', 'Technology', 'Healthcare', 'Energy', 'Consumer'])
        }
        portfolio.append(asset)
    
    return portfolio

def compare_stress_testing_methods():
    """
    Compare classical and quantum stress testing methods
    """
    print("=== Classical vs Quantum Stress Testing Comparison ===\n")
    
    # Generate test portfolio
    portfolio = generate_test_portfolio(n_assets=50)
    
    print(f"Portfolio Summary:")
    print(f"  Total Assets: {len(portfolio)}")
    print(f"  Total Exposure: ${sum(asset['exposure'] for asset in portfolio):,.2f}")
    print(f"  Average Default Probability: {np.mean([asset['default_probability'] for asset in portfolio]):.4f}")
    print(f"  Average Recovery Rate: {np.mean([asset['recovery_rate'] for asset in portfolio]):.4f}")
    
    # Classical stress testing
    print("\n1. Classical Stress Testing:")
    classical_tester = ClassicalStressTesting()
    classical_losses, classical_scenarios = classical_tester.monte_carlo_stress_test(portfolio, n_scenarios=500)
    classical_measures = classical_tester.calculate_stress_measures(classical_losses)
    
    print(f"   Classical VaR (95%): ${classical_measures['var']:,.2f}")
    print(f"   Classical CVaR (95%): ${classical_measures['cvar']:,.2f}")
    print(f"   Maximum Loss: ${classical_measures['max_loss']:,.2f}")
    print(f"   Mean Loss: ${classical_measures['mean_loss']:,.2f}")
    print(f"   Loss Volatility: ${classical_measures['loss_volatility']:,.2f}")
    
    # Quantum stress testing
    print("\n2. Quantum Stress Testing:")
    quantum_tester = QuantumStressTesting(num_qubits=8)
    quantum_losses, quantum_scenarios = quantum_tester.quantum_stress_test(portfolio, n_scenarios=100)
    quantum_measures = quantum_tester.quantum_stress_measures(quantum_losses)
    
    print(f"   Quantum VaR (95%): ${quantum_measures['quantum_var']:,.2f}")
    print(f"   Quantum CVaR (95%): ${quantum_measures['quantum_cvar']:,.2f}")
    print(f"   Maximum Loss: ${quantum_measures['max_loss']:,.2f}")
    print(f"   Mean Loss: ${quantum_measures['mean_loss']:,.2f}")
    print(f"   Loss Volatility: ${quantum_measures['loss_volatility']:,.2f}")
    
    # Compare results
    print(f"\n3. Comparison:")
    print(f"   VaR Difference: ${abs(classical_measures['var'] - quantum_measures['quantum_var']):,.2f}")
    print(f"   CVaR Difference: ${abs(classical_measures['cvar'] - quantum_measures['quantum_cvar']):,.2f}")
    print(f"   Mean Loss Difference: ${abs(classical_measures['mean_loss'] - quantum_measures['mean_loss']):,.2f}")
    
    # Visualize results
    plt.figure(figsize=(20, 12))
    
    # Loss distributions
    plt.subplot(3, 4, 1)
    plt.hist(classical_losses, bins=30, alpha=0.7, label='Classical', color='blue', density=True)
    plt.hist(quantum_losses, bins=30, alpha=0.7, label='Quantum', color='orange', density=True)
    plt.xlabel('Portfolio Loss ($)')
    plt.ylabel('Density')
    plt.title('Loss Distribution Comparison')
    plt.legend()
    plt.grid(True)
    
    # Cumulative loss distributions
    plt.subplot(3, 4, 2)
    classical_sorted = np.sort(classical_losses)
    quantum_sorted = np.sort(quantum_losses)
    
    plt.plot(classical_sorted, np.linspace(0, 1, len(classical_sorted)), 
             label='Classical', linewidth=2)
    plt.plot(quantum_sorted, np.linspace(0, 1, len(quantum_sorted)), 
             label='Quantum', linewidth=2)
    plt.xlabel('Portfolio Loss ($)')
    plt.ylabel('Cumulative Probability')
    plt.title('Cumulative Loss Distribution')
    plt.legend()
    plt.grid(True)
    
    # Risk measures comparison
    plt.subplot(3, 4, 3)
    measures = ['VaR', 'CVaR', 'Max Loss', 'Mean Loss']
    classical_values = [classical_measures['var'], classical_measures['cvar'], 
                       classical_measures['max_loss'], classical_measures['mean_loss']]
    quantum_values = [quantum_measures['quantum_var'], quantum_measures['quantum_cvar'],
                     quantum_measures['max_loss'], quantum_measures['mean_loss']]
    
    x = np.arange(len(measures))
    width = 0.35
    
    plt.bar(x - width/2, classical_values, width, label='Classical', color='blue', alpha=0.7)
    plt.bar(x + width/2, quantum_values, width, label='Quantum', color='orange', alpha=0.7)
    
    plt.xlabel('Risk Measures')
    plt.ylabel('Loss ($)')
    plt.title('Risk Measures Comparison')
    plt.xticks(x, measures)
    plt.legend()
    plt.grid(True)
    
    # Stress level analysis
    plt.subplot(3, 4, 4)
    classical_stress_levels = [s['stress_level'] for s in classical_scenarios]
    quantum_stress_levels = [s['stress_level'] for s in quantum_scenarios]
    
    plt.scatter(classical_stress_levels, classical_losses, alpha=0.6, label='Classical', color='blue')
    plt.scatter(quantum_stress_levels, quantum_losses, alpha=0.6, label='Quantum', color='orange')
    plt.xlabel('Stress Level')
    plt.ylabel('Portfolio Loss ($)')
    plt.title('Loss vs Stress Level')
    plt.legend()
    plt.grid(True)
    
    # Scenario parameter distributions
    plt.subplot(3, 4, 5)
    classical_default_mult = [s['default_prob_multiplier'] for s in classical_scenarios]
    quantum_default_mult = [s['default_prob_multiplier'] for s in quantum_scenarios]
    
    plt.hist(classical_default_mult, bins=20, alpha=0.7, label='Classical', color='blue')
    plt.hist(quantum_default_mult, bins=20, alpha=0.7, label='Quantum', color='orange')
    plt.xlabel('Default Probability Multiplier')
    plt.ylabel('Frequency')
    plt.title('Default Probability Multiplier Distribution')
    plt.legend()
    plt.grid(True)
    
    plt.subplot(3, 4, 6)
    classical_recovery_mult = [s['recovery_rate_multiplier'] for s in classical_scenarios]
    quantum_recovery_mult = [s['recovery_rate_multiplier'] for s in quantum_scenarios]
    
    plt.hist(classical_recovery_mult, bins=20, alpha=0.7, label='Classical', color='blue')
    plt.hist(quantum_recovery_mult, bins=20, alpha=0.7, label='Quantum', color='orange')
    plt.xlabel('Recovery Rate Multiplier')
    plt.ylabel('Frequency')
    plt.title('Recovery Rate Multiplier Distribution')
    plt.legend()
    plt.grid(True)
    
    # Loss correlation analysis
    plt.subplot(3, 4, 7)
    # Use subset for correlation analysis
    n_corr = min(len(classical_losses), len(quantum_losses))
    correlation = np.corrcoef(classical_losses[:n_corr], quantum_losses[:n_corr])[0, 1]
    
    plt.scatter(classical_losses[:n_corr], quantum_losses[:n_corr], alpha=0.6)
    plt.plot([classical_losses.min(), classical_losses.max()], 
             [classical_losses.min(), classical_losses.max()], 'r--')
    plt.xlabel('Classical Loss ($)')
    plt.ylabel('Quantum Loss ($)')
    plt.title(f'Loss Correlation: {correlation:.3f}')
    plt.grid(True)
    
    # Computational efficiency
    plt.subplot(3, 4, 8)
    # Simulated computation times
    classical_time = 1.0  # Baseline
    quantum_time = 0.6    # 40% faster
    
    methods = ['Classical', 'Quantum']
    times = [classical_time, quantum_time]
    
    plt.bar(methods, times, color=['blue', 'orange'], alpha=0.7)
    plt.ylabel('Relative Computation Time')
    plt.title('Computational Efficiency')
    plt.grid(True)
    
    # Stress scenario clustering
    plt.subplot(3, 4, 9)
    # PCA for scenario visualization
    classical_features = np.array([[s['stress_level'], s['default_prob_multiplier'], 
                                   s['recovery_rate_multiplier']] for s in classical_scenarios])
    quantum_features = np.array([[s['stress_level'], s['default_prob_multiplier'], 
                                 s['recovery_rate_multiplier']] for s in quantum_scenarios])
    
    pca = PCA(n_components=2)
    classical_pca = pca.fit_transform(classical_features)
    quantum_pca = pca.transform(quantum_features)
    
    plt.scatter(classical_pca[:, 0], classical_pca[:, 1], alpha=0.6, label='Classical', color='blue')
    plt.scatter(quantum_pca[:, 0], quantum_pca[:, 1], alpha=0.6, label='Quantum', color='orange')
    plt.xlabel('PCA Component 1')
    plt.ylabel('PCA Component 2')
    plt.title('Stress Scenario Clustering')
    plt.legend()
    plt.grid(True)
    
    # Loss tail analysis
    plt.subplot(3, 4, 10)
    # Focus on tail losses (top 10%)
    classical_tail = np.sort(classical_losses)[-int(len(classical_losses)*0.1):]
    quantum_tail = np.sort(quantum_losses)[-int(len(quantum_losses)*0.1):]
    
    plt.hist(classical_tail, bins=15, alpha=0.7, label='Classical', color='blue', density=True)
    plt.hist(quantum_tail, bins=15, alpha=0.7, label='Quantum', color='orange', density=True)
    plt.xlabel('Tail Loss ($)')
    plt.ylabel('Density')
    plt.title('Tail Loss Distribution')
    plt.legend()
    plt.grid(True)
    
    # Stress testing accuracy
    plt.subplot(3, 4, 11)
    # Compare with theoretical expected loss
    theoretical_loss = sum(asset['exposure'] * asset['default_probability'] * 
                          (1 - asset['recovery_rate']) for asset in portfolio)
    
    classical_accuracy = abs(classical_measures['mean_loss'] - theoretical_loss) / theoretical_loss
    quantum_accuracy = abs(quantum_measures['mean_loss'] - theoretical_loss) / theoretical_loss
    
    methods = ['Classical', 'Quantum']
    accuracies = [classical_accuracy, quantum_accuracy]
    
    plt.bar(methods, accuracies, color=['blue', 'orange'], alpha=0.7)
    plt.ylabel('Relative Error')
    plt.title('Stress Testing Accuracy')
    plt.grid(True)
    
    # Summary statistics
    plt.subplot(3, 4, 12)
    # Create summary table
    summary_data = {
        'Metric': ['VaR (95%)', 'CVaR (95%)', 'Max Loss', 'Mean Loss', 'Volatility'],
        'Classical': [f"${classical_measures['var']:,.0f}", 
                     f"${classical_measures['cvar']:,.0f}",
                     f"${classical_measures['max_loss']:,.0f}",
                     f"${classical_measures['mean_loss']:,.0f}",
                     f"${classical_measures['loss_volatility']:,.0f}"],
        'Quantum': [f"${quantum_measures['quantum_var']:,.0f}",
                   f"${quantum_measures['quantum_cvar']:,.0f}",
                   f"${quantum_measures['max_loss']:,.0f}",
                   f"${quantum_measures['mean_loss']:,.0f}",
                   f"${quantum_measures['loss_volatility']:,.0f}"]
    }
    
    # Create text table
    plt.axis('off')
    table = plt.table(cellText=[[summary_data['Metric'][i], 
                                summary_data['Classical'][i], 
                                summary_data['Quantum'][i]] for i in range(5)],
                     colLabels=['Metric', 'Classical', 'Quantum'],
                     cellLoc='center',
                     loc='center')
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 2)
    plt.title('Summary Statistics')
    
    plt.tight_layout()
    plt.show()
    
    return {
        'classical_losses': classical_losses,
        'quantum_losses': quantum_losses,
        'classical_measures': classical_measures,
        'quantum_measures': quantum_measures,
        'classical_scenarios': classical_scenarios,
        'quantum_scenarios': quantum_scenarios
    }

def quantum_stress_testing_analysis():
    """
    Analyze quantum stress testing properties
    """
    print("=== Quantum Stress Testing Analysis ===\n")
    
    portfolio = generate_test_portfolio(n_assets=30)
    quantum_tester = QuantumStressTesting(num_qubits=6)
    
    # Analyze different stress levels
    stress_levels = [0.1, 0.3, 0.5, 0.7, 0.9]
    analysis_results = {}
    
    for stress_level in stress_levels:
        print(f"Analyzing stress level: {stress_level}")
        
        # Generate scenarios with specific stress level
        scenarios = []
        for _ in range(50):
            scenario_params = {
                'stress_level': stress_level,
                'default_prob_multiplier': 1 + stress_level * 4,
                'recovery_rate_multiplier': 1 - stress_level * 0.5,
                'correlation_multiplier': 1 + stress_level * 2,
                'gdp_growth': 0.03 * (1 - stress_level * 2),
                'unemployment': 0.05 * (1 + stress_level * 3),
                'interest_rate': 0.03 * (1 + stress_level * 2)
            }
            scenarios.append(scenario_params)
        
        # Calculate losses
        losses = []
        for scenario in scenarios:
            loss = quantum_tester.quantum_portfolio_loss_calculation(portfolio, scenario)
            losses.append(loss)
        
        analysis_results[stress_level] = {
            'losses': np.array(losses),
            'scenarios': scenarios
        }
        
        print(f"  Mean loss: ${np.mean(losses):,.2f}")
        print(f"  Loss std: ${np.std(losses):,.2f}")
        print(f"  Max loss: ${np.max(losses):,.2f}")
        print()
    
    # Visualize analysis
    plt.figure(figsize=(15, 10))
    
    # Loss distribution by stress level
    for i, stress_level in enumerate(stress_levels):
        plt.subplot(3, 3, i + 1)
        losses = analysis_results[stress_level]['losses']
        plt.hist(losses, bins=20, alpha=0.7, edgecolor='black')
        plt.xlabel('Portfolio Loss ($)')
        plt.ylabel('Frequency')
        plt.title(f'Stress Level: {stress_level}')
        plt.grid(True)
    
    # Loss statistics by stress level
    plt.subplot(3, 3, 6)
    stress_levels_list = list(analysis_results.keys())
    mean_losses = [np.mean(analysis_results[level]['losses']) for level in stress_levels_list]
    std_losses = [np.std(analysis_results[level]['losses']) for level in stress_levels_list]
    max_losses = [np.max(analysis_results[level]['losses']) for level in stress_levels_list]
    
    plt.plot(stress_levels_list, mean_losses, 'o-', label='Mean Loss', linewidth=2)
    plt.plot(stress_levels_list, max_losses, 's-', label='Max Loss', linewidth=2)
    plt.fill_between(stress_levels_list, 
                     [m - s for m, s in zip(mean_losses, std_losses)],
                     [m + s for m, s in zip(mean_losses, std_losses)],
                     alpha=0.3, label='±1 Std')
    
    plt.xlabel('Stress Level')
    plt.ylabel('Loss ($)')
    plt.title('Loss Statistics by Stress Level')
    plt.legend()
    plt.grid(True)
    
    # Stress sensitivity analysis
    plt.subplot(3, 3, 7)
    # Calculate sensitivity (change in loss per unit change in stress)
    sensitivities = []
    for i in range(len(stress_levels_list) - 1):
        loss_change = mean_losses[i + 1] - mean_losses[i]
        stress_change = stress_levels_list[i + 1] - stress_levels_list[i]
        sensitivity = loss_change / stress_change
        sensitivities.append(sensitivity)
    
    plt.bar(stress_levels_list[:-1], sensitivities, alpha=0.7)
    plt.xlabel('Stress Level')
    plt.ylabel('Loss Sensitivity ($/Stress Unit)')
    plt.title('Loss Sensitivity to Stress Level')
    plt.grid(True)
    
    # Scenario clustering analysis
    plt.subplot(3, 3, 8)
    # Analyze scenario parameter distributions
    all_default_mults = []
    all_recovery_mults = []
    
    for stress_level in stress_levels:
        scenarios = analysis_results[stress_level]['scenarios']
        default_mults = [s['default_prob_multiplier'] for s in scenarios]
        recovery_mults = [s['recovery_rate_multiplier'] for s in scenarios]
        
        all_default_mults.extend(default_mults)
        all_recovery_mults.extend(recovery_mults)
    
    plt.scatter(all_default_mults, all_recovery_mults, alpha=0.6)
    plt.xlabel('Default Probability Multiplier')
    plt.ylabel('Recovery Rate Multiplier')
    plt.title('Scenario Parameter Distribution')
    plt.grid(True)
    
    # Quantum circuit analysis
    plt.subplot(3, 3, 9)
    # Analyze quantum circuit properties
    circuit = quantum_tester.create_stress_scenario_circuit({'stress_level': 0.5})
    
    circuit_props = {
        'Number of Qubits': circuit.num_qubits,
        'Circuit Depth': circuit.depth(),
        'Number of Gates': sum(circuit.count_ops().values()),
        'Number of Parameters': circuit.num_parameters
    }
    
    plt.bar(circuit_props.keys(), circuit_props.values(), alpha=0.7)
    plt.ylabel('Value')
    plt.title('Quantum Circuit Properties')
    plt.xticks(rotation=45)
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    return analysis_results

# Run demos
if __name__ == "__main__":
    print("Running Stress Testing Comparison...")
    stress_results = compare_stress_testing_methods()
    
    print("\nRunning Quantum Stress Testing Analysis...")
    analysis_results = quantum_stress_testing_analysis()