import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, confusion_matrix
import statsmodels.api as sm
from scipy import stats

# Create the dataset from the table in Image 1
data = {
    'x1': [-62.8, 3.3, -120.8, -18.1, -3.8, -61.2, -20.3, -194.5, 20.8, -106.1, -39.4, -164.1, -308.9, 7.2, -118.3, -185.9, -34.6, -27.9, -48.2, -49.2, -19.2, -18.1, -98.0, -129.0, -4.0, -8.7, -59.2, -13.1, -38.0, -57.9, -8.8, -64.7, -11.4, 43.0, 47.0, -3.3, 35.0, 46.7, 20.8, 33.0, 26.1, 68.6, 37.3, 59.0, 49.6, 12.5, 37.3, 35.3, 48.5, 18.1, 31.4, 21.5, 8.5, 40.6, 34.6, 19.9, 17.4, 54.7, 53.5, 35.9, 39.4, 53.1, 39.8, 59.5, 16.3, 21.7],
    'x2': [-89.5, -3.5, -103.2, -28.8, -50.6, -56.2, -17.4, -25.8, -4.3, -22.9, -35.7, -17.7, -65.8, -22.6, -34.2, -280.0, -19.4, 6.3, 6.8, -17.2, -36.7, -6.5, -20.8, -14.2, -15.8, -36.3, -12.8, -17.6, 1.6, 0.7, -9.1, -0.4, 4.8, 16.4, 16.0, 4.0, 20.8, 12.6, 12.5, 23.6, 10.4, 13.8, 33.4, 23.1, 23.8, 7.0, 34.1, 4.2, 25.1, 13.5, 15.7, -14.4, 5.8, 5.8, 26.4, 26.7, 12.6, 14.6, 20.6, 26.4, 30.5, 7.1, 13.8, 7.0, 20.4, -7.8],
    'x3': [1.7, 1.1, 2.5, 1.1, 0.9, 1.7, 1.0, 0.5, 1.0, 1.5, 1.2, 1.3, 0.8, 2.0, 1.5, 6.7, 3.4, 1.3, 1.6, 0.3, 0.8, 0.9, 1.7, 1.3, 2.1, 2.8, 2.1, 0.9, 1.2, 0.8, 0.9, 0.1, 0.9, 1.3, 1.9, 2.7, 1.9, 0.9, 2.4, 1.5, 2.1, 1.6, 3.5, 5.5, 1.9, 1.8, 1.5, 0.9, 2.6, 4.0, 1.9, 1.0, 1.5, 1.8, 1.8, 2.3, 1.3, 1.7, 1.1, 2.0, 1.9, 1.9, 1.2, 2.0, 1.0, 1.6],
    'y': [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
}

df = pd.DataFrame(data)
print("Dataset shape:", df.shape)
print("\nData summary:")
print(df.describe())

# Check class distribution
print("\nBankruptcy status distribution:")
print(df['y'].value_counts())

# 1. Build logistic regression model with all three variables X1, X2, X3
# Add constant term for statsmodels
X = sm.add_constant(df[['x1', 'x2', 'x3']])
y = df['y']

# Fit model using statsmodels for detailed statistics
model_full = sm.Logit(y, X).fit()
print("\nFull model summary:")
print(model_full.summary())

# 2. Test significance of X3 using likelihood ratio test
X_reduced = sm.add_constant(df[['x1', 'x2']])
model_reduced = sm.Logit(y, X_reduced).fit()

# Likelihood ratio test
lr_stat = -2 * (model_reduced.llf - model_full.llf)
p_value = 1 - stats.chi2.cdf(lr_stat, df=1)  # df is the difference in parameters (1 here)

print("\nLikelihood ratio test for X3:")
print(f"LR statistic: {lr_stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"X3 is {'significant' if p_value < 0.05 else 'not significant'} at α=0.05 level")

# If X3 is not significant, build a model with just X1 and X2
print("\nReduced model summary (X1 and X2 only):")
print(model_reduced.summary())

# Calculate predicted probabilities for all companies
df['pred_prob'] = model_reduced.predict(X_reduced)

# Fit using sklearn for alternative metrics
X_sk = df[['x1', 'x2']]
model_sk = LogisticRegression(random_state=42)
model_sk.fit(X_sk, y)

print("\nClassification metrics:")
y_pred = model_sk.predict(X_sk)
print(classification_report(y, y_pred))

# 3. Predict probability for the specific company
specific_company = {'x1': 48.8, 'x2': -10.5, 'x3': 1.8}

# For reduced model (X1, X2) - FIXED THIS PART
X_new_reduced = pd.DataFrame({'const': [1], 'x1': [specific_company['x1']], 'x2': [specific_company['x2']]})
prob_reduced = model_reduced.predict(X_new_reduced)[0]

# For full model (X1, X2, X3) - FIXED THIS PART
X_new_full = pd.DataFrame({
    'const': [1], 
    'x1': [specific_company['x1']], 
    'x2': [specific_company['x2']],
    'x3': [specific_company['x3']]
})
prob_full = model_full.predict(X_new_full)[0]

print("\nPrediction for specific company (x1=48.8, x2=-10.5, x3=1.8):")
print(f"Probability of survival (full model): {prob_full:.4f}")
print(f"Probability of survival (reduced model): {prob_reduced:.4f}")
print(f"Difference in probabilities: {abs(prob_full - prob_reduced):.4f}")

# Visualize the data and decision boundary
plt.figure(figsize=(12, 10))

# Plot 1: X1 vs X2 colored by Y
plt.subplot(2, 2, 1)
plt.scatter(df[df['y']==0]['x1'], df[df['y']==0]['x2'], c='red', label='Bankrupt (Y=0)')
plt.scatter(df[df['y']==1]['x1'], df[df['y']==1]['x2'], c='blue', label='Solvent (Y=1)')
plt.scatter(specific_company['x1'], specific_company['x2'], c='green', marker='*', s=200, label='Test Company')
plt.xlabel('X1 (Retained Earnings/Total Assets)')
plt.ylabel('X2 (Profit before Interest & Tax/Total Assets)')
plt.title('Financial Indicators X1 vs X2')
plt.legend()
plt.grid(True)

# Plot 2: Predicted probabilities
plt.subplot(2, 2, 2)
y_probs = df['pred_prob']
plt.scatter(df[df['y']==0].index, df[df['y']==0]['pred_prob'], c='red', label='Bankrupt (Y=0)')
plt.scatter(df[df['y']==1].index, df[df['y']==1]['pred_prob'], c='blue', label='Solvent (Y=1)')
plt.axhline(y=0.5, color='k', linestyle='--', label='Decision Boundary (0.5)')
plt.xlabel('Company Index')
plt.ylabel('Predicted Probability of Solvency P(Y=1)')
plt.title('Predicted Probabilities')
plt.legend()
plt.grid(True)

# Plot 3: X1 distribution by class
plt.subplot(2, 2, 3)
plt.hist(df[df['y']==0]['x1'], alpha=0.5, color='red', bins=10, label='Bankrupt (Y=0)')
plt.hist(df[df['y']==1]['x1'], alpha=0.5, color='blue', bins=10, label='Solvent (Y=1)')
plt.axvline(x=specific_company['x1'], color='green', linestyle='--', label='Test Company')
plt.xlabel('X1 (Retained Earnings/Total Assets)')
plt.ylabel('Frequency')
plt.title('Distribution of X1 by Bankruptcy Status')
plt.legend()
plt.grid(True)

# Plot 4: X2 distribution by class
plt.subplot(2, 2, 4)
plt.hist(df[df['y']==0]['x2'], alpha=0.5, color='red', bins=10, label='Bankrupt (Y=0)')
plt.hist(df[df['y']==1]['x2'], alpha=0.5, color='blue', bins=10, label='Solvent (Y=1)')
plt.axvline(x=specific_company['x2'], color='green', linestyle='--', label='Test Company')
plt.xlabel('X2 (Profit before Interest & Tax/Total Assets)')
plt.ylabel('Frequency')
plt.title('Distribution of X2 by Bankruptcy Status')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.savefig('bankruptcy_analysis.png')
plt.show()