# ============= HYPOTHESIS TESTING IN FINANCE =============

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import ttest_ind, mannwhitneyu, f_oneway, chi2_contingency
import matplotlib.pyplot as plt
import seaborn as sns

# Create sample financial data
np.random.seed(42)

# Customer data with different segments
premium_customers = pd.DataFrame({
    'segment': 'Premium',
    'annual_income': np.random.normal(120000, 30000, 200).clip(50000, 250000),
    'credit_score': np.random.normal(750, 40, 200).clip(650, 850),
    'monthly_spend': np.random.normal(3000, 800, 200).clip(500, 6000),
    'default_rate': np.random.beta(1, 10, 200)  # Low default rate
})

standard_customers = pd.DataFrame({
    'segment': 'Standard',
    'annual_income': np.random.normal(70000, 20000, 300).clip(30000, 150000),
    'credit_score': np.random.normal(680, 50, 300).clip(550, 800),
    'monthly_spend': np.random.normal(1800, 600, 300).clip(200, 4000),
    'default_rate': np.random.beta(3, 8, 300)  # Medium default rate
})

basic_customers = pd.DataFrame({
    'segment': 'Basic',
    'annual_income': np.random.normal(40000, 15000, 200).clip(15000, 80000),
    'credit_score': np.random.normal(620, 50, 200).clip(500, 750),
    'monthly_spend': np.random.normal(900, 400, 200).clip(100, 2500),
    'default_rate': np.random.beta(5, 5, 200)  # Higher default rate
})

# Combine data
all_customers = pd.concat([premium_customers, standard_customers, basic_customers], ignore_index=True)

class FinancialHypothesisTester:
    """Perform hypothesis tests for financial data."""
    
    def __init__(self, data):
        self.data = data
        self.results = {}
    
    def test_segment_income_difference(self):
        """Test if income differs significantly between segments."""
        premium_income = self.data[self.data['segment'] == 'Premium']['annual_income']
        standard_income = self.data[self.data['segment'] == 'Standard']['annual_income']
        basic_income = self.data[self.data['segment'] == 'Basic']['annual_income']
        
        # T-test: Premium vs Standard
        t_stat, p_value = ttest_ind(premium_income, standard_income)
        
        self.results['income_premium_vs_standard'] = {
            'test': 'Independent T-Test',
            't_statistic': t_stat,
            'p_value': p_value,
            'significant': p_value < 0.05,
            'interpretation': f"Premium customers have {'significantly different' if p_value < 0.05 else 'similar'} income compared to Standard"
        }
        
        # T-test: Standard vs Basic
        t_stat, p_value = ttest_ind(standard_income, basic_income)
        
        self.results['income_standard_vs_basic'] = {
            'test': 'Independent T-Test',
            't_statistic': t_stat,
            'p_value': p_value,
            'significant': p_value < 0.05,
            'interpretation': f"Standard customers have {'significantly different' if p_value < 0.05 else 'similar'} income compared to Basic"
        }
        
        # ANOVA: All three segments
        f_stat, p_value = f_oneway(premium_income, standard_income, basic_income)
        
        self.results['income_all_segments'] = {
            'test': 'ANOVA',
            'f_statistic': f_stat,
            'p_value': p_value,
            'significant': p_value < 0.05,
            'interpretation': f"There is {'a significant' if p_value < 0.05 else 'no'} difference in income across all segments"
        }
        
        return self.results
    
    def test_credit_score_difference(self):
        """Test if credit scores differ significantly between segments."""
        premium_scores = self.data[self.data['segment'] == 'Premium']['credit_score']
        standard_scores = self.data[self.data['segment'] == 'Standard']['credit_score']
        basic_scores = self.data[self.data['segment'] == 'Basic']['credit_score']
        
        # Mann-Whitney U (non-parametric, for non-normal distributions)
        u_stat, p_value = mannwhitneyu(premium_scores, standard_scores)
        
        self.results['credit_premium_vs_standard'] = {
            'test': 'Mann-Whitney U',
            'u_statistic': u_stat,
            'p_value': p_value,
            'significant': p_value < 0.05,
            'interpretation': f"Premium customers have {'significantly different' if p_value < 0.05 else 'similar'} credit scores compared to Standard"
        }
        
        return self.results
    
    def test_default_rate_difference(self):
        """Test if default rates differ significantly between segments."""
        premium_defaults = self.data[self.data['segment'] == 'Premium']['default_rate']
        standard_defaults = self.data[self.data['segment'] == 'Standard']['default_rate']
        basic_defaults = self.data[self.data['segment'] == 'Basic']['default_rate']
        
        # Kruskal-Wallis test (non-parametric ANOVA)
        h_stat, p_value = stats.kruskal(premium_defaults, standard_defaults, basic_defaults)
        
        self.results['default_rates'] = {
            'test': 'Kruskal-Wallis',
            'h_statistic': h_stat,
            'p_value': p_value,
            'significant': p_value < 0.05,
            'interpretation': f"There is {'a significant' if p_value < 0.05 else 'no'} difference in default rates across segments",
            'segment_means': {
                'Premium': premium_defaults.mean(),
                'Standard': standard_defaults.mean(),
                'Basic': basic_defaults.mean()
            }
        }
        
        return self.results
    
    def test_correlation_significance(self, var1, var2):
        """Test if correlation between two variables is significant."""
        correlation, p_value = stats.pearsonr(self.data[var1], self.data[var2])
        
        self.results[f'correlation_{var1}_{var2}'] = {
            'test': 'Pearson Correlation',
            'correlation': correlation,
            'p_value': p_value,
            'significant': p_value < 0.05,
            'interpretation': f"There is {'a significant' if p_value < 0.05 else 'no'} correlation between {var1} and {var2}"
        }
        
        return self.results
    
    def run_all_tests(self):
        """Run all hypothesis tests."""
        self.test_segment_income_difference()
        self.test_credit_score_difference()
        self.test_default_rate_difference()
        self.test_correlation_significance('annual_income', 'credit_score')
        self.test_correlation_significance('monthly_spend', 'default_rate')
        
        return self.results
    
    def report_results(self):
        """Generate a comprehensive report of test results."""
        print("\n" + "="*80)
        print("HYPOTHESIS TESTING REPORT")
        print("="*80)
        
        # Run tests if not already run
        if not self.results:
            self.run_all_tests()
        
        for test_name, result in self.results.items():
            print(f"\n📊 {test_name.replace('_', ' ').title()}")
            print(f"  Test: {result.get('test', 'N/A')}")
            print(f"  Statistic: {result.get('t_statistic', result.get('u_statistic', result.get('f_statistic', result.get('correlation', 'N/A')))):.4f}")
            print(f"  P-value: {result['p_value']:.6f}")
            print(f"  Significant: {'✅ YES' if result['significant'] else '❌ NO'} (p < 0.05)")
            print(f"  Interpretation: {result['interpretation']}")
            
            if 'segment_means' in result:
                print("  Segment Means:")
                for segment, mean in result['segment_means'].items():
                    print(f"    {segment}: {mean:.4f}")
        
        # Summary
        significant_tests = sum(1 for r in self.results.values() if r.get('significant', False))
        total_tests = len(self.results)
        
        print("\n" + "-"*80)
        print(f"SUMMARY: {significant_tests}/{total_tests} tests showed statistically significant results")

# Run hypothesis tests
tester = FinancialHypothesisTester(all_customers)
tester.run_all_tests()
tester.report_results()

# Visualize the differences
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Income by segment
sns.boxplot(data=all_customers, x='segment', y='annual_income', ax=axes[0, 0])
axes[0, 0].set_title('Income Distribution by Customer Segment', fontsize=14)
axes[0, 0].set_ylabel('Annual Income ($)')

# Credit score by segment
sns.boxplot(data=all_customers, x='segment', y='credit_score', ax=axes[0, 1])
axes[0, 1].set_title('Credit Score Distribution by Customer Segment', fontsize=14)
axes[0, 1].set_ylabel('Credit Score')

# Default rate by segment
sns.boxplot(data=all_customers, x='segment', y='default_rate', ax=axes[1, 0])
axes[1, 0].set_title('Default Rate by Customer Segment', fontsize=14)
axes[1, 0].set_ylabel('Default Rate')

# Correlation: Income vs Credit Score
sns.scatterplot(data=all_customers, x='annual_income', y='credit_score', hue='segment', alpha=0.6, ax=axes[1, 1])
axes[1, 1].set_title('Income vs Credit Score by Segment', fontsize=14)
axes[1, 1].set_xlabel('Annual Income ($)')
axes[1, 1].set_ylabel('Credit Score')

plt.tight_layout()
plt.savefig('hypothesis_testing_visualization.png', dpi=300)
plt.show()

# ============= CONFIDENCE INTERVALS =============

class FinancialConfidenceIntervals:
    """Calculate confidence intervals for financial metrics."""
    
    def __init__(self, data):
        self.data = data
        self.results = {}
    
    def confidence_interval_mean(self, column, confidence=0.95):
        """Calculate confidence interval for the mean of a column."""
        n = len(self.data[column])
        mean = self.data[column].mean()
        std = self.data[column].std()
        
        # Standard error
        se = std / np.sqrt(n)
        
        # Z-score for confidence level
        z_score = stats.norm.ppf((1 + confidence) / 2)
        
        # Confidence interval
        ci_lower = mean - z_score * se
        ci_upper = mean + z_score * se
        
        self.results[f'ci_mean_{column}'] = {
            'mean': mean,
            'std': std,
            'n': n,
            'confidence': confidence,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper,
            'interpretation': f"We are {confidence*100:.0f}% confident that the true mean of {column} is between {ci_lower:.2f} and {ci_upper:.2f}"
        }
        
        return self.results[f'ci_mean_{column}']
    
    def confidence_interval_proportion(self, column, value, confidence=0.95):
        """Calculate confidence interval for a proportion."""
        n = len(self.data[column])
        p = (self.data[column] == value).mean()
        
        # Standard error for proportion
        se = np.sqrt(p * (1 - p) / n)
        
        # Z-score for confidence level
        z_score = stats.norm.ppf((1 + confidence) / 2)
        
        # Confidence interval
        ci_lower = p - z_score * se
        ci_upper = p + z_score * se
        
        # Clamp to [0, 1]
        ci_lower = max(0, ci_lower)
        ci_upper = min(1, ci_upper)
        
        self.results[f'ci_proportion_{column}_{value}'] = {
            'proportion': p,
            'n': n,
            'confidence': confidence,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper,
            'interpretation': f"We are {confidence*100:.0f}% confident that the true proportion of {column} = {value} is between {ci_lower:.2%} and {ci_upper:.2%}"
        }
        
        return self.results[f'ci_proportion_{column}_{value}']
    
    def bootstrap_ci(self, column, statistic='mean', n_bootstrap=1000, confidence=0.95):
        """Calculate bootstrap confidence interval."""
        values = self.data[column].values
        bootstrap_stats = []
        
        for _ in range(n_bootstrap):
            # Sample with replacement
            sample = np.random.choice(values, size=len(values), replace=True)
            
            if statistic == 'mean':
                bootstrap_stats.append(sample.mean())
            elif statistic == 'median':
                bootstrap_stats.append(np.median(sample))
            elif statistic == 'std':
                bootstrap_stats.append(sample.std())
        
        # Percentile confidence interval
        alpha = 1 - confidence
        ci_lower = np.percentile(bootstrap_stats, alpha/2 * 100)
        ci_upper = np.percentile(bootstrap_stats, (1 - alpha/2) * 100)
        
        self.results[f'bootstrap_ci_{column}_{statistic}'] = {
            'statistic': statistic,
            'bootstrap_mean': np.mean(bootstrap_stats),
            'n_bootstrap': n_bootstrap,
            'confidence': confidence,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper,
            'interpretation': f"Bootstrap CI ({confidence*100:.0f}%) for {statistic} of {column}: [{ci_lower:.2f}, {ci_upper:.2f}]"
        }
        
        return self.results[f'bootstrap_ci_{column}_{statistic}']
    
    def report_results(self):
        """Generate report of confidence intervals."""
        print("\n" + "="*80)
        print("CONFIDENCE INTERVAL REPORT")
        print("="*80)
        
        for result_name, result in self.results.items():
            print(f"\n📊 {result_name.replace('_', ' ').title()}")
            print(f"  {result['interpretation']}")
            if 'mean' in result:
                print(f"  Sample mean: {result['mean']:.2f}")
            if 'proportion' in result:
                print(f"  Sample proportion: {result['proportion']:.2%}")

# Calculate confidence intervals
ci_calculator = FinancialConfidenceIntervals(all_customers)

# Mean confidence intervals
ci_calculator.confidence_interval_mean('annual_income', confidence=0.95)
ci_calculator.confidence_interval_mean('credit_score', confidence=0.95)
ci_calculator.confidence_interval_mean('monthly_spend', confidence=0.95)

# Proportion confidence intervals
ci_calculator.confidence_interval_proportion('segment', 'Premium', confidence=0.95)

# Bootstrap confidence intervals
ci_calculator.bootstrap_ci('annual_income', statistic='mean', n_bootstrap=1000)
ci_calculator.bootstrap_ci('annual_income', statistic='median', n_bootstrap=1000)

ci_calculator.report_results()

# ============= MULTIVARIATE CORRELATION ANALYSIS =============

# Select numerical columns for multivariate analysis
numeric_cols = ['annual_income', 'credit_score', 'monthly_spend', 'default_rate', 'age']
correlation_matrix = all_customers[numeric_cols].corr()

print("\n" + "="*80)
print("CORRELATION MATRIX")
print("="*80)
print(correlation_matrix.round(3))

# Create detailed correlation visualization
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# 1. Correlation heatmap with annotations
sns.heatmap(correlation_matrix, annot=True, cmap='RdBu_r', center=0, 
            fmt='.3f', square=True, ax=axes[0, 0],
            cbar_kws={'shrink': 0.8})
axes[0, 0].set_title('Correlation Heatmap - Customer Features', fontsize=14)

# 2. Pairplot (scatter matrix)
# Create a simplified pairplot manually for selected columns
selected_cols = ['annual_income', 'credit_score', 'monthly_spend']
subset = all_customers[selected_cols]

# Create scatter matrix
from pandas.plotting import scatter_matrix
scatter_matrix(subset, alpha=0.3, figsize=(8, 8), diagonal='hist', ax=axes[0, 1])
axes[0, 1].set_title('Scatter Matrix - Customer Features', fontsize=14)

# 3. Correlation dendrogram (hierarchical clustering of variables)
from scipy.cluster import hierarchy
correlation_linkage = hierarchy.linkage(correlation_matrix, method='average')
dendro = hierarchy.dendrogram(correlation_linkage, ax=axes[1, 0], 
                             labels=correlation_matrix.columns, 
                             orientation='top')
axes[1, 0].set_title('Variable Clustering Dendrogram', fontsize=14)

# 4. Correlation with segment (using point-biserial correlation)
from scipy.stats import pointbiserialr

# Convert segment to numeric
segment_mapping = {'Premium': 2, 'Standard': 1, 'Basic': 0}
all_customers['segment_numeric'] = all_customers['segment'].map(segment_mapping)

correlations_with_segment = {}
for col in numeric_cols:
    corr, p_val = pointbiserialr(all_customers['segment_numeric'], all_customers[col])
    correlations_with_segment[col] = {'correlation': corr, 'p_value': p_val}

# Plot correlations with segment
seg_corr_values = [correlations_with_segment[col]['correlation'] for col in numeric_cols]
bars = axes[1, 1].bar(numeric_cols, seg_corr_values, color='steelblue')
axes[1, 1].set_title('Correlation with Customer Segment', fontsize=14)
axes[1, 1].set_ylabel('Point-Biserial Correlation')
axes[1, 1].axhline(0, color='black', linestyle='-', alpha=0.3)
axes[1, 1].axhline(0.2, color='red', linestyle='--', alpha=0.5, label='Moderate correlation')
axes[1, 1].axhline(-0.2, color='red', linestyle='--', alpha=0.5)
axes[1, 1].legend()
axes[1, 1].tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.savefig('multivariate_analysis.png', dpi=300)
plt.show()

# Print correlation insights
print("\n" + "="*80)
print("MULTIVARIATE ANALYSIS INSIGHTS")
print("="*80)

# Find strongest correlations
corr_pairs = []
for i in range(len(correlation_matrix.columns)):
    for j in range(i+1, len(correlation_matrix.columns)):
        corr_pairs.append({
            'pair': f"{correlation_matrix.columns[i]} - {correlation_matrix.columns[j]}",
            'correlation': correlation_matrix.iloc[i, j]
        })

corr_pairs = sorted(corr_pairs, key=lambda x: abs(x['correlation']), reverse=True)

print("\nStrongest Correlations:")
for pair in corr_pairs[:5]:
    print(f"  • {pair['pair']}: {pair['correlation']:.3f}")

# Correlation with segment
print("\nCorrelation with Customer Segment:")
for col, data in correlations_with_segment.items():
    print(f"  • {col}: {data['correlation']:.3f} (p={data['p_value']:.4f})")
    if data['p_value'] < 0.05:
        print(f"    → Significant relationship with segment")

# ============= PRINCIPAL COMPONENT ANALYSIS =============

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# Prepare data for PCA
pca_data = all_customers[numeric_cols].copy()

# Scale the data
scaler = StandardScaler()
pca_data_scaled = scaler.fit_transform(pca_data)

# Perform PCA
pca = PCA()
pca_result = pca.fit_transform(pca_data_scaled)

# Create PCA results DataFrame
pca_df = pd.DataFrame(pca_result, columns=[f'PC{i+1}' for i in range(len(numeric_cols))])
pca_df['segment'] = all_customers['segment'].values

# Analyze PCA results
explained_variance = pca.explained_variance_ratio_
cumulative_variance = np.cumsum(explained_variance)

print("\n" + "="*80)
print("PRINCIPAL COMPONENT ANALYSIS (PCA)")
print("="*80)

print("\nExplained Variance Ratio:")
for i, var in enumerate(explained_variance):
    print(f"  PC{i+1}: {var:.2%} (Cumulative: {cumulative_variance[i]:.2%})")

print(f"\nTotal variance explained by first 2 components: {cumulative_variance[1]:.2%}")
print(f"Total variance explained by first 3 components: {cumulative_variance[2]:.2%}")

# PCA loadings
loadings = pd.DataFrame(
    pca.components_.T,
    columns=[f'PC{i+1}' for i in range(len(numeric_cols))],
    index=numeric_cols
)

print("\nPCA Loadings (Feature Contributions):")
print(loadings.round(3))

# Visualize PCA
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# 1. Scree plot
axes[0, 0].bar(range(1, len(explained_variance)+1), explained_variance, alpha=0.7, label='Individual')
axes[0, 0].plot(range(1, len(explained_variance)+1), cumulative_variance, 'ro-', label='Cumulative')
axes[0, 0].set_xlabel('Principal Component')
axes[0, 0].set_ylabel('Explained Variance Ratio')
axes[0, 0].set_title('Scree Plot', fontsize=14)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. PCA scatter plot (PC1 vs PC2)
for segment in ['Premium', 'Standard', 'Basic']:
    segment_data = pca_df[pca_df['segment'] == segment]
    axes[0, 1].scatter(segment_data['PC1'], segment_data['PC2'], 
                       label=segment, alpha=0.6, s=50)
axes[0, 1].set_xlabel(f'PC1 ({explained_variance[0]:.1%} variance)')
axes[0, 1].set_ylabel(f'PC2 ({explained_variance[1]:.1%} variance)')
axes[0, 1].set_title('PCA - Customer Segmentation', fontsize=14)
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# 3. PCA loadings heatmap
sns.heatmap(loadings, annot=True, cmap='RdBu_r', center=0, fmt='.2f', 
            ax=axes[1, 0], cbar_kws={'shrink': 0.8})
axes[1, 0].set_title('PCA Loadings', fontsize=14)

# 4. Biplot (PCA with feature vectors)
# Simplified biplot: show first two PCs with feature vectors
for col in numeric_cols:
    axes[1, 1].arrow(0, 0, 
                     loadings.loc[col, 'PC1'] * 3,
                     loadings.loc[col, 'PC2'] * 3,
                     head_width=0.1, head_length=0.1, fc='red', ec='red')
    axes[1, 1].text(loadings.loc[col, 'PC1'] * 3.2,
                    loadings.loc[col, 'PC2'] * 3.2,
                    col, fontsize=10)

# Add data points
for segment in ['Premium', 'Standard', 'Basic']:
    segment_data = pca_df[pca_df['segment'] == segment]
    axes[1, 1].scatter(segment_data['PC1'], segment_data['PC2'], 
                       label=segment, alpha=0.3, s=30)

axes[1, 1].set_xlabel('PC1')
axes[1, 1].set_ylabel('PC2')
axes[1, 1].set_title('Biplot - PC1 vs PC2', fontsize=14)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('pca_analysis.png', dpi=300)
plt.show()

# PCA insights
print("\n" + "="*80)
print("PCA INSIGHTS")
print("="*80)

# Most important features for PC1
pc1_loadings = loadings['PC1'].abs().sort_values(ascending=False)
print(f"\nPC1 (Explains {explained_variance[0]:.1%} of variance):")
print("  Top features:")
for feature, loading in pc1_loadings.head(3).items():
    print(f"    • {feature}: {loading:.3f}")

# Most important features for PC2
pc2_loadings = loadings['PC2'].abs().sort_values(ascending=False)
print(f"\nPC2 (Explains {explained_variance[1]:.1%} of variance):")
print("  Top features:")
for feature, loading in pc2_loadings.head(3).items():
    print(f"    • {feature}: {loading:.3f}")

# ============= TIME SERIES DECOMPOSITION =============

from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
import warnings
warnings.filterwarnings('ignore')

# Create daily transaction data
np.random.seed(42)
dates = pd.date_range(start='2023-01-01', end='2024-12-31', freq='D')
n_days = len(dates)

# Create a series with trend, seasonality, and noise
trend = np.linspace(1000, 5000, n_days)
seasonality = 500 * np.sin(2 * np.pi * dates.dayofyear / 365)
seasonality += 200 * np.sin(4 * np.pi * dates.dayofyear / 365)  # Weekly pattern
noise = np.random.normal(0, 200, n_days)

# Add day-of-week effect
day_of_week = dates.dayofweek
weekday_effect = np.array([100, 150, 200, 250, 300, 100, 50])[day_of_week]

transaction_volume = trend + seasonality + noise + weekday_effect
transaction_volume = np.maximum(transaction_volume, 0)  # No negative values

# Create DataFrame
ts_data = pd.DataFrame({
    'date': dates,
    'transaction_volume': transaction_volume
})
ts_data.set_index('date', inplace=True)

print("\n" + "="*80)
print("TIME SERIES ANALYSIS")
print("="*80)

# 1. Stationarity test (Augmented Dickey-Fuller)
adf_result = adfuller(ts_data['transaction_volume'], autolag='AIC')
print(f"\nStationarity Test (ADF):")
print(f"  Test Statistic: {adf_result[0]:.4f}")
print(f"  P-value: {adf_result[1]:.4f}")
print(f"  Critical Values:")
for key, value in adf_result[4].items():
    print(f"    {key}: {value:.4f}")
print(f"  {'Stationary' if adf_result[1] < 0.05 else 'Non-Stationary'}")

# 2. Decompose the time series
decomposition = seasonal_decompose(ts_data['transaction_volume'], model='additive', period=365)

# Extract components
trend_component = decomposition.trend
seasonal_component = decomposition.seasonal
residual_component = decomposition.resid

# Visualize decomposition
fig, axes = plt.subplots(4, 1, figsize=(15, 12))

# Original series
axes[0].plot(ts_data.index, ts_data['transaction_volume'], color='black', alpha=0.7)
axes[0].set_title('Original Transaction Volume', fontsize=14)
axes[0].set_ylabel('Volume')
axes[0].grid(True, alpha=0.3)

# Trend
axes[1].plot(ts_data.index, trend_component, color='blue', linewidth=2)
axes[1].set_title('Trend Component', fontsize=14)
axes[1].set_ylabel('Trend')
axes[1].grid(True, alpha=0.3)

# Seasonal
axes[2].plot(ts_data.index, seasonal_component, color='green', linewidth=2)
axes[2].set_title('Seasonal Component (365-day cycle)', fontsize=14)
axes[2].set_ylabel('Seasonal')
axes[2].grid(True, alpha=0.3)

# Residual
axes[3].plot(ts_data.index, residual_component, color='red', alpha=0.7)
axes[3].set_title('Residual (Noise)', fontsize=14)
axes[3].set_ylabel('Residual')
axes[3].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('time_series_decomposition.png', dpi=300)
plt.show()

# 3. Weekly pattern analysis
ts_data['day_of_week'] = ts_data.index.dayofweek
day_names = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']

weekly_pattern = ts_data.groupby('day_of_week')['transaction_volume'].agg(['mean', 'std'])

print("\nWeekly Pattern Analysis:")
print(weekly_pattern)

# 4. Monthly pattern analysis
ts_data['month'] = ts_data.index.month
monthly_pattern = ts_data.groupby('month')['transaction_volume'].agg(['mean', 'std'])

print("\nMonthly Pattern Analysis:")
print(monthly_pattern)

# 5. Detect seasonality strength
seasonal_strength = 1 - (np.var(residual_component.dropna()) / np.var(ts_data['transaction_volume'].dropna()))
print(f"\nSeasonality Strength: {seasonal_strength:.2%}")
print(f"  • If close to 1: Strong seasonality")
print(f"  • If close to 0: Weak or no seasonality")

# 6. Identify peaks and troughs
seasonal_peak = seasonal_component.max()
seasonal_trough = seasonal_component.min()
seasonal_amplitude = seasonal_peak - seasonal_trough

print(f"\nSeasonal Characteristics:")
print(f"  Peak seasonality: {seasonal_peak:.0f}")
print(f"  Trough seasonality: {seasonal_trough:.0f}")
print(f"  Amplitude: {seasonal_amplitude:.0f}")

# ============= AUTOCORRELATION ANALYSIS =============

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Calculate autocorrelation
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# 1. Autocorrelation Function (ACF)
plot_acf(ts_data['transaction_volume'].dropna(), lags=60, ax=axes[0, 0])
axes[0, 0].set_title('Autocorrelation Function (ACF)', fontsize=14)

# 2. Partial Autocorrelation Function (PACF)
plot_pacf(ts_data['transaction_volume'].dropna(), lags=60, ax=axes[0, 1])
axes[0, 1].set_title('Partial Autocorrelation Function (PACF)', fontsize=14)

# 3. Autocorrelation at different lags
from statsmodels.tsa.stattools import acf

lags_to_check = [1, 7, 14, 30, 60, 90, 365]
acf_values = acf(ts_data['transaction_volume'].dropna(), nlags=365, fft=False)

for lag in lags_to_check:
    axes[1, 0].bar(lag, acf_values[lag], color='steelblue', width=5)

axes[1, 0].set_xlabel('Lag')
axes[1, 0].set_ylabel('Autocorrelation')
axes[1, 0].set_title('Selected Lag Autocorrelations', fontsize=14)
axes[1, 0].axhline(0, color='black', linestyle='-')
axes[1, 0].axhline(1.96/np.sqrt(len(ts_data)), color='red', linestyle='--', alpha=0.5, label='95% CI')
axes[1, 0].axhline(-1.96/np.sqrt(len(ts_data)), color='red', linestyle='--', alpha=0.5)
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 4. Heatmap of autocorrelation matrix
from scipy.linalg import toeplitz
acf_vals = acf(ts_data['transaction_volume'].dropna(), nlags=30)
acf_matrix = toeplitz(acf_vals)

sns.heatmap(acf_matrix, cmap='RdBu_r', center=0, ax=axes[1, 1],
            cbar_kws={'shrink': 0.8}, square=True)
axes[1, 1].set_title('Autocorrelation Matrix (30 lags)', fontsize=14)

plt.tight_layout()
plt.savefig('autocorrelation_analysis.png', dpi=300)
plt.show()

# Interpret autocorrelation
print("\n" + "="*80)
print("AUTOCORRELATION INTERPRETATION")
print("="*80)

print("\n🔍 Key Observations:")
for lag in lags_to_check:
    if lag < len(acf_values):
        value = acf_values[lag]
        if abs(value) > 1.96/np.sqrt(len(ts_data)):
            significance = "significant"
            if value > 0.3:
                strength = "strong"
            elif value > 0.1:
                strength = "moderate"
            else:
                strength = "weak"
            print(f"  • Lag {lag}: {value:.3f} - {strength} {significance} positive autocorrelation")
        else:
            print(f"  • Lag {lag}: {value:.3f} - not significant")

# Seasonal period detection
# Find peaks in autocorrelation that correspond to seasonal patterns
seasonal_lags = []
for lag in range(30, 400):
    if lag < len(acf_values) and acf_values[lag] > 0.2:
        if len(seasonal_lags) == 0 or lag - seasonal_lags[-1] > 10:
            seasonal_lags.append(lag)

print(f"\n📊 Detected seasonal lags: {seasonal_lags[:5]}")
if seasonal_lags:
    print(f"  Likely seasonal period: {seasonal_lags[0]} days")

# ============= ADVANCED VISUALIZATION =============

# Create additional data for comprehensive visualization
all_customers['age'] = np.random.normal(45, 15, len(all_customers)).clip(18, 80).astype(int)
all_customers['monthly_spend'] = all_customers['monthly_spend'] + np.random.normal(0, 200, len(all_customers)).clip(-500, 500)

# 1. Violin Plots
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# Income distribution by segment
sns.violinplot(data=all_customers, x='segment', y='annual_income', ax=axes[0, 0])
axes[0, 0].set_title('Income Distribution by Segment (Violin Plot)', fontsize=14)
axes[0, 0].set_ylabel('Annual Income ($)')

# Credit score distribution by segment
sns.violinplot(data=all_customers, x='segment', y='credit_score', ax=axes[0, 1])
axes[0, 1].set_title('Credit Score Distribution by Segment (Violin Plot)', fontsize=14)
axes[0, 1].set_ylabel('Credit Score')

# Monthly spend distribution by segment
sns.violinplot(data=all_customers, x='segment', y='monthly_spend', ax=axes[1, 0])
axes[1, 0].set_title('Monthly Spend Distribution by Segment (Violin Plot)', fontsize=14)
axes[1, 0].set_ylabel('Monthly Spend ($)')

# Age distribution by segment
sns.violinplot(data=all_customers, x='segment', y='age', ax=axes[1, 1])
axes[1, 1].set_title('Age Distribution by Segment (Violin Plot)', fontsize=14)
axes[1, 1].set_ylabel('Age')

plt.tight_layout()
plt.savefig('violin_plots.png', dpi=300)
plt.show()

# ============= PARALLEL COORDINATES PLOT =============

from pandas.plotting import parallel_coordinates

# Prepare data for parallel coordinates
parallel_data = all_customers[['segment', 'annual_income', 'credit_score', 'monthly_spend', 'default_rate']].copy()

# Scale the data for better visualization
for col in ['annual_income', 'credit_score', 'monthly_spend', 'default_rate']:
    parallel_data[col] = (parallel_data[col] - parallel_data[col].min()) / (parallel_data[col].max() - parallel_data[col].min())

# Create parallel coordinates plot
fig, ax = plt.subplots(figsize=(14, 8))

parallel_coordinates(parallel_data, 'segment', ax=ax, colormap='Set2')
ax.set_title('Parallel Coordinates - Customer Segments', fontsize=14)
ax.set_xlabel('Features')
ax.set_ylabel('Normalized Value')

plt.tight_layout()
plt.savefig('parallel_coordinates.png', dpi=300)
plt.show()

# ============= ADVANCED PAIRPLOT =============

import seaborn as sns

# Create pairplot with additional features
pairplot_cols = ['annual_income', 'credit_score', 'monthly_spend', 'default_rate', 'age']

# Use a subset for readability
subset_data = all_customers.sample(n=200, random_state=42)

g = sns.pairplot(subset_data[pairplot_cols], 
                 diag_kind='kde',
                 corner=True,
                 plot_kws={'alpha': 0.5, 's': 30})

g.fig.suptitle('Pairplot - Customer Features', y=1.02, fontsize=16)
plt.tight_layout()
plt.savefig('advanced_pairplot.png', dpi=300)
plt.show()

# ============= CUSTOMER SEGMENTATION =============

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# Prepare data for clustering
cluster_features = ['annual_income', 'credit_score', 'monthly_spend', 'age']

# Standardize the features
scaler = StandardScaler()
scaled_data = scaler.fit_transform(all_customers[cluster_features])

# Elbow method to find optimal K
inertias = []
K_range = range(2, 10)

for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(scaled_data)
    inertias.append(kmeans.inertia_)

# Perform final clustering with K=4
kmeans = KMeans(n_clusters=4, random_state=42, n_init=10)
cluster_labels = kmeans.fit_predict(scaled_data)
all_customers['cluster'] = cluster_labels

# Analyze clusters
cluster_analysis = all_customers.groupby('cluster')[cluster_features].mean()
cluster_counts = all_customers['cluster'].value_counts()

print("\n" + "="*80)
print("CUSTOMER SEGMENTATION - CLUSTER ANALYSIS")
print("="*80)

print("\nCluster Profiles:")
for cluster in sorted(cluster_analysis.index):
    print(f"\n🔵 Cluster {cluster}:")
    print(f"  Size: {cluster_counts[cluster]} customers ({cluster_counts[cluster]/len(all_customers)*100:.1f}%)")
    for feature in cluster_features:
        print(f"  {feature}: ${cluster_analysis.loc[cluster, feature]:,.2f}" if 'income' in feature or 'spend' in feature else f"  {feature}: {cluster_analysis.loc[cluster, feature]:.1f}")

# Visualize clusters
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# 1. Elbow curve
axes[0, 0].plot(K_range, inertias, 'bo-')
axes[0, 0].set_xlabel('Number of Clusters (K)')
axes[0, 0].set_ylabel('Inertia')
axes[0, 0].set_title('Elbow Method for Optimal K', fontsize=14)
axes[0, 0].grid(True, alpha=0.3)

# 2. Cluster distribution
cluster_counts.plot(kind='bar', ax=axes[0, 1], color=['#2ecc71', '#3498db', '#e74c3c', '#f39c12'])
axes[0, 1].set_title('Cluster Distribution', fontsize=14)
axes[0, 1].set_xlabel('Cluster')
axes[0, 1].set_ylabel('Count')

# 3. Cluster profiles (heatmap)
cluster_profile_heatmap = cluster_analysis.T
sns.heatmap(cluster_profile_heatmap, annot=True, fmt='.0f', ax=axes[1, 0],
            cmap='RdYlGn', cbar_kws={'shrink': 0.8})
axes[1, 0].set_title('Cluster Profiles Heatmap', fontsize=14)

# 4. 2D projection (using PCA)
from sklearn.decomposition import PCA
pca_2d = PCA(n_components=2)
pca_result_2d = pca_2d.fit_transform(scaled_data)

axes[1, 1].scatter(pca_result_2d[:, 0], pca_result_2d[:, 1], 
                  c=cluster_labels, cmap='Set2', alpha=0.6, s=50)
axes[1, 1].set_title('Customer Clusters (PCA Projection)', fontsize=14)
axes[1, 1].set_xlabel('PC1')
axes[1, 1].set_ylabel('PC2')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('customer_segmentation.png', dpi=300)
plt.show()

# Segment names
cluster_segments = {
    0: 'High Income, High Credit',
    1: 'Medium Income, Standard Credit',
    2: 'Low Income, Lower Credit',
    3: 'High Spenders'
}

print("\n" + "="*80)
print("CLUSTER SEGMENT NAMES")
print("="*80)
for cluster in sorted(cluster_analysis.index):
    print(f"\nCluster {cluster}: {cluster_segments.get(cluster, 'Unnamed')}")
    print(f"  Key characteristics:")
    for feature in cluster_features:
        if feature in ['annual_income', 'monthly_spend']:
            value = cluster_analysis.loc[cluster, feature]
            avg = all_customers[feature].mean()
            if value > avg * 1.2:
                direction = "high"
            elif value < avg * 0.8:
                direction = "low"
            else:
                direction = "average"
            print(f"    • {feature}: {direction} (${value:,.0f} vs avg ${avg:,.0f})")
        elif feature == 'credit_score':
            value = cluster_analysis.loc[cluster, feature]
            avg = all_customers[feature].mean()
            if value > avg * 1.1:
                direction = "high"
            elif value < avg * 0.9:
                direction = "low"
            else:
                direction = "average"
            print(f"    • {feature}: {direction} ({value:.0f} vs avg {avg:.0f})")

# ============= COMPREHENSIVE EDA REPORT =============

def generate_comprehensive_eda_report(data, dataset_name):
    """Generate a comprehensive EDA report."""
    
    print("\n" + "="*80)
    print(f"COMPREHENSIVE EDA REPORT: {dataset_name}")
    print("="*80)
    
    # 1. Data Overview
    print("\n📊 SECTION 1: DATA OVERVIEW")
    print("-"*40)
    print(f"  Dataset size: {data.shape[0]} rows, {data.shape[1]} columns")
    print(f"  Memory usage: {data.memory_usage(deep=True).sum() / 1024**2:.2f} MB")
    
    # 2. Data Quality
    print("\n✅ SECTION 2: DATA QUALITY")
    print("-"*40)
    missing = data.isnull().sum()
    missing_pct = missing / len(data) * 100
    if missing.sum() > 0:
        print("  Missing Data:")
        for col in data.columns:
            if missing[col] > 0:
                print(f"    • {col}: {missing[col]} missing ({missing_pct[col]:.1f}%)")
    else:
        print("  ✓ No missing data found")
    
    # 3. Summary Statistics
    print("\n📈 SECTION 3: SUMMARY STATISTICS")
    print("-"*40)
    numeric_cols = data.select_dtypes(include=[np.number]).columns
    for col in numeric_cols[:5]:  # Show first 5 numeric columns
        print(f"\n  {col}:")
        print(f"    Mean: {data[col].mean():.2f}")
        print(f"    Median: {data[col].median():.2f}")
        print(f"    Std Dev: {data[col].std():.2f}")
        print(f"    Min: {data[col].min():.2f}")
        print(f"    Max: {data[col].max():.2f}")
        print(f"    Skewness: {data[col].skew():.2f}")
        print(f"    Kurtosis: {data[col].kurtosis():.2f}")
    
    # 4. Key Insights
    print("\n💡 SECTION 4: KEY INSIGHTS")
    print("-"*40)
    
    # Distribution insights
    for col in numeric_cols:
        skewness = data[col].skew()
        if abs(skewness) > 1:
            print(f"  • {col} is highly skewed ({skewness:.2f}) - may need transformation")
        elif abs(skewness) > 0.5:
            print(f"  • {col} is moderately skewed ({skewness:.2f})")
    
    # Outlier insights
    print("\n  Outlier Detection (IQR method):")
    for col in numeric_cols:
        q1 = data[col].quantile(0.25)
        q3 = data[col].quantile(0.75)
        iqr = q3 - q1
        outliers = data[(data[col] < q1 - 1.5*iqr) | (data[col] > q3 + 1.5*iqr)]
        if len(outliers) > 0:
            print(f"    • {col}: {len(outliers)} outliers ({len(outliers)/len(data)*100:.1f}%)")
    
    # 5. Recommendations
    print("\n🎯 SECTION 5: RECOMMENDATIONS")
    print("-"*40)
    
    recommendations = []
    
    # Check for missing data
    if missing.sum() > 0:
        recommendations.append("• Handle missing data using appropriate imputation methods")
    
    # Check for skewness
    for col in numeric_cols:
        if abs(data[col].skew()) > 1:
            recommendations.append(f"• Apply transformation (log/sqrt) to {col} due to high skewness")
            break
    
    # Check for outliers
    for col in numeric_cols:
        q1 = data[col].quantile(0.25)
        q3 = data[col].quantile(0.75)
        iqr = q3 - q1
        outliers = data[(data[col] < q1 - 1.5*iqr) | (data[col] > q3 + 1.5*iqr)]
        if len(outliers) > 0:
            recommendations.append(f"• Consider winsorizing or transforming {col} to handle outliers")
            break
    
    # Check for high correlation
    if len(numeric_cols) > 1:
        corr_matrix = data[numeric_cols].corr()
        for i in range(len(corr_matrix.columns)):
            for j in range(i+1, len(corr_matrix.columns)):
                if abs(corr_matrix.iloc[i, j]) > 0.8:
                    recommendations.append(f"• Consider removing redundant features with high correlation ({corr_matrix.columns[i]} vs {corr_matrix.columns[j]}: {corr_matrix.iloc[i, j]:.2f})")
                    break
    
    for rec in recommendations:
        print(f"  {rec}")
    
    if not recommendations:
        print("  ✓ Data quality is good - ready for modeling")
    
    # 6. EDA Checklist
    print("\n📋 SECTION 6: EDA CHECKLIST")
    print("-"*40)
    checklist = [
        ("Data loaded and inspected", True),
        ("Summary statistics calculated", True),
        ("Missing data identified and analyzed", missing.sum() > 0),
        ("Outliers detected", any(abs(data[col].skew()) > 1 for col in numeric_cols)),
        ("Distributions visualized", True),
        ("Correlation analysis performed", len(numeric_cols) > 1),
        ("Time series patterns analyzed", 'date' in data.columns),
        ("Segmentation analysis performed", True),
        ("Key insights documented", True),
        ("Recommendations provided", True)
    ]
    
    for item, status in checklist:
        status_icon = "✅" if status else "⏳"
        print(f"  {status_icon} {item}")
    
    print("\n" + "="*80)
    print("END OF EDA REPORT")
    print("="*80)

# Generate comprehensive report
generate_comprehensive_eda_report(all_customers, "Customer Banking Data")