clean_frequency_figure.py

"""cleaner version of the frequency analysis figure"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import pearsonr
import warnings
warnings.filterwarnings('ignore')

# styling
sns.set_style("whitegrid")
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial']

print("Loading data...")

# load the data
disasters = pd.read_excel('naturalDisasters.xlsx')
gdp_growth = pd.read_csv('gdp_growth.csv')
gdp = pd.read_csv('gdp.csv')

# clean it up
disasters['Start Year'] = pd.to_numeric(disasters['Start Year'], errors='coerce')

# see what years we have
disasters_filtered = disasters[(disasters['Start Year'] >= 1960) & 
                                (disasters['Start Year'] <= 2020)].copy()
actual_min_year = disasters_filtered['Start Year'].min()
actual_max_year = disasters_filtered['Start Year'].max()

print(f"Actual disaster data coverage: {actual_min_year:.0f} - {actual_max_year:.0f}")

# only use years with data
disasters_filtered = disasters[(disasters['Start Year'] >= actual_min_year) & 
                                (disasters['Start Year'] <= actual_max_year)].copy()

annual_disasters = disasters_filtered.groupby('Start Year').size().reset_index(name='Disaster_Frequency')
annual_disasters.columns = ['Year', 'Disaster_Frequency']

print(f"Years with disaster data: {len(annual_disasters)}")

# get GDP data for the same years
year_cols = [col for col in gdp_growth.columns if col.isdigit() and 
             actual_min_year <= int(col) <= actual_max_year]

global_gdp_data = []
for year_col in year_cols:
    year = int(year_col)
    gdp_growth_values = gdp_growth[year_col].dropna()
    
    global_gdp_data.append({
        'Year': year,
        'Global_GDP_Growth': gdp_growth_values.mean()
    })

gdp_df = pd.DataFrame(global_gdp_data)

# combine them (only years with both)
merged_data = annual_disasters.merge(gdp_df, on='Year', how='inner')
merged_data = merged_data.sort_values('Year').reset_index(drop=True)
merged_data['Decade'] = (merged_data['Year'] // 10) * 10

print(f"Merged data: {len(merged_data)} years")
print(f"Year range: {merged_data['Year'].min():.0f} - {merged_data['Year'].max():.0f}")
print()

# get the correlation
freq_gdp_data = merged_data[['Disaster_Frequency', 'Global_GDP_Growth']].dropna()
corr_freq, p_freq = pearsonr(freq_gdp_data['Disaster_Frequency'], 
                               freq_gdp_data['Global_GDP_Growth'])

print(f"Correlation: r = {corr_freq:.4f}, p = {p_freq:.4f}")
print()
print("Creating clean frequency analysis figure...")

# make the figure
fig = plt.figure(figsize=(16, 10))
gs = fig.add_gridspec(2, 2, hspace=0.35, wspace=0.30, 
                      left=0.08, right=0.95, top=0.92, bottom=0.08)

# colors
color_disasters = '#E74C3C'  # Red
color_gdp = '#27AE60'  # Green
color_scatter = '#3498DB'  # Blue

# ============================================================================
# panel A: disaster frequency over time
# ============================================================================
ax1 = fig.add_subplot(gs[0, 0])

# plot it
ax1.plot(merged_data['Year'], merged_data['Disaster_Frequency'], 
         linewidth=2.5, color=color_disasters, marker='o', markersize=4, 
         markerfacecolor='white', markeredgewidth=1.5, markeredgecolor=color_disasters,
         alpha=0.8, label='Annual Frequency')

# add trend line
valid_data = merged_data[merged_data['Disaster_Frequency'] > 0]
if len(valid_data) > 3:
    z = np.polyfit(valid_data['Year'], valid_data['Disaster_Frequency'], 2)
    p = np.poly1d(z)
    trend_years = np.linspace(valid_data['Year'].min(), valid_data['Year'].max(), 100)
    ax1.plot(trend_years, p(trend_years), '--', color='#C0392B', 
             linewidth=3, label='Polynomial Trend', alpha=0.7)

ax1.set_xlabel('Year', fontsize=13, fontweight='bold')
ax1.set_ylabel('Number of Disasters', fontsize=13, fontweight='bold')
ax1.set_title('(A) Annual Disaster Frequency', fontsize=14, fontweight='bold', pad=15)
ax1.legend(fontsize=11, framealpha=0.95, loc='upper left')
ax1.grid(True, alpha=0.25, linestyle='--', linewidth=0.8)
ax1.set_xlim(merged_data['Year'].min() - 1, merged_data['Year'].max() + 1)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

# ============================================================================
# panel B: GDP growth over time
# ============================================================================
ax2 = fig.add_subplot(gs[0, 1])

# plot GDP
ax2.plot(merged_data['Year'], merged_data['Global_GDP_Growth'], 
         linewidth=2.5, color=color_gdp, marker='s', markersize=4,
         markerfacecolor='white', markeredgewidth=1.5, markeredgecolor=color_gdp,
         alpha=0.8, label='GDP Growth Rate')

# line at zero
ax2.axhline(y=0, color='#95A5A6', linestyle='-', linewidth=1.5, alpha=0.7, zorder=1)

# shade positive/negative areas
ax2.fill_between(merged_data['Year'], merged_data['Global_GDP_Growth'], 0,
                  where=(merged_data['Global_GDP_Growth'] > 0), 
                  alpha=0.15, color=color_gdp, interpolate=True, label='Positive Growth')
ax2.fill_between(merged_data['Year'], merged_data['Global_GDP_Growth'], 0,
                  where=(merged_data['Global_GDP_Growth'] <= 0), 
                  alpha=0.2, color='#E74C3C', interpolate=True, label='Negative Growth')

ax2.set_xlabel('Year', fontsize=13, fontweight='bold')
ax2.set_ylabel('GDP Growth Rate (%)', fontsize=13, fontweight='bold')
ax2.set_title('(B) Global GDP Growth Rate', fontsize=14, fontweight='bold', pad=15)
ax2.legend(fontsize=11, framealpha=0.95, loc='lower left')
ax2.grid(True, alpha=0.25, linestyle='--', linewidth=0.8)
ax2.set_xlim(merged_data['Year'].min() - 1, merged_data['Year'].max() + 1)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

# ============================================================================
# panel C: scatter with regression line
# ============================================================================
ax3 = fig.add_subplot(gs[1, 0])

# scatter plot colored by year
scatter = ax3.scatter(merged_data['Disaster_Frequency'], 
                      merged_data['Global_GDP_Growth'],
                      c=merged_data['Year'], cmap='viridis', 
                      s=100, alpha=0.7, edgecolors='black', linewidth=1,
                      vmin=1960, vmax=2020)

# add the line
if len(freq_gdp_data) > 2:
    slope, intercept, r_value, p_value, std_err = stats.linregress(
        freq_gdp_data['Disaster_Frequency'], 
        freq_gdp_data['Global_GDP_Growth'])
    
    x_line = np.linspace(0, merged_data['Disaster_Frequency'].max(), 100)
    y_line = intercept + slope * x_line
    ax3.plot(x_line, y_line, 'r--', linewidth=3, alpha=0.8,
             label=f'Linear Fit (R² = {r_value**2:.3f})')
    
    # show the stats
    stats_text = f'r = {corr_freq:.3f}\np = {p_freq:.3f}\n'
    stats_text += f'{"Significant" if p_freq < 0.05 else "Not Significant"} (α=0.05)'
    ax3.text(0.05, 0.95, stats_text, transform=ax3.transAxes,
             fontsize=11, verticalalignment='top',
             bbox=dict(boxstyle='round,pad=0.7', facecolor='white', 
                      edgecolor='gray', alpha=0.9))

ax3.set_xlabel('Disaster Frequency', fontsize=13, fontweight='bold')
ax3.set_ylabel('GDP Growth Rate (%)', fontsize=13, fontweight='bold')
ax3.set_title('(C) Frequency vs GDP Growth', fontsize=14, fontweight='bold', pad=15)
ax3.legend(fontsize=11, framealpha=0.95, loc='lower right')
ax3.grid(True, alpha=0.25, linestyle='--', linewidth=0.8)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)

# add colorbar
cbar = plt.colorbar(scatter, ax=ax3, pad=0.02)
cbar.set_label('Year', fontsize=11, fontweight='bold')
cbar.ax.tick_params(labelsize=10)

# ============================================================================
# panel D: compare decades
# ============================================================================
ax4 = fig.add_subplot(gs[1, 1])

# average by decade
decade_avg = merged_data.groupby('Decade').agg({
    'Disaster_Frequency': 'mean',
    'Global_GDP_Growth': 'mean'
}).reset_index()

x_pos = np.arange(len(decade_avg))
width = 0.35

# make bar chart
bars1 = ax4.bar(x_pos - width/2, decade_avg['Disaster_Frequency'], 
                width, label='Avg Disasters/Year', 
                color=color_disasters, alpha=0.8, edgecolor='black', linewidth=1.2)

# second y-axis for GDP
ax4_twin = ax4.twinx()
bars2 = ax4_twin.bar(x_pos + width/2, decade_avg['Global_GDP_Growth'], 
                     width, label='Avg GDP Growth (%)', 
                     color=color_gdp, alpha=0.8, edgecolor='black', linewidth=1.2)

# labels
ax4.set_xlabel('Decade', fontsize=13, fontweight='bold')
ax4.set_ylabel('Avg Disaster Frequency', fontsize=12, fontweight='bold', color=color_disasters)
ax4_twin.set_ylabel('Avg GDP Growth (%)', fontsize=12, fontweight='bold', color=color_gdp)
ax4.set_title('(D) Decade-Level Comparison', fontsize=14, fontweight='bold', pad=15)

ax4.set_xticks(x_pos)
ax4.set_xticklabels([f"{int(d)}s" for d in decade_avg['Decade']], fontsize=11)
ax4.tick_params(axis='y', labelcolor=color_disasters, labelsize=10)
ax4_twin.tick_params(axis='y', labelcolor=color_gdp, labelsize=10)

# clean up the look
ax4.grid(True, alpha=0.25, linestyle='--', linewidth=0.8, axis='y')
ax4.spines['top'].set_visible(False)
ax4_twin.spines['top'].set_visible(False)

# combine the legends
lines1, labels1 = ax4.get_legend_handles_labels()
lines2, labels2 = ax4_twin.get_legend_handles_labels()
ax4.legend(lines1 + lines2, labels1 + labels2, loc='upper left', 
          fontsize=11, framealpha=0.95)

# label the bars
for bar in bars1:
    height = bar.get_height()
    if height > 0:
        ax4.text(bar.get_x() + bar.get_width()/2., height * 0.5,
                f'{height:.0f}', ha='center', va='center', 
                fontsize=9, fontweight='bold', color='white')

for bar in bars2:
    height = bar.get_height()
    # put label inside if positive, outside if negative
    if height > 0:
        y_pos = height * 0.5
        va_align = 'center'
    else:
        y_pos = height - 0.3
        va_align = 'top'
    ax4_twin.text(bar.get_x() + bar.get_width()/2., y_pos,
                 f'{height:.1f}%', ha='center', va=va_align, 
                 fontsize=9, fontweight='bold', color='white')

# ============================================================================
# add title
# ============================================================================
year_range = f"{int(merged_data['Year'].min())}-{int(merged_data['Year'].max())}"
fig.suptitle(f'Question 1: Impact of Disaster Frequency on Global GDP Growth ({year_range})', 
             fontsize=16, fontweight='bold', y=0.98)

# save it
plt.savefig('figure1_frequency_analysis.png', dpi=300, bbox_inches='tight', 
            facecolor='white', edgecolor='none')
print("✓ Clean frequency analysis figure saved: figure1_frequency_analysis.png")
print()

# save high-res PDF too
plt.savefig('figure1_frequency_analysis_highres.pdf', dpi=600, bbox_inches='tight', 
            facecolor='white', edgecolor='none')
print("✓ High-resolution PDF saved: figure1_frequency_analysis_highres.pdf")

plt.show()