Aula 10: Painel de Dados na Inferência Causal

# Implementação de modelos de painel em Python
import pandas as pd
import numpy as np
from linearmodels import PanelOLS, RandomEffects
from linearmodels.panel import compare
import statsmodels.api as sm

# Carregar dados de painel
df = pd.read_csv('dados_empresas_painel.csv')
df = df.set_index(['empresa', 'trimestre'])

# Preparar variáveis
y = df['roe']
X = df[['alavancagem', 'tamanho', 'liquidez', 'crescimento']]
X = sm.add_constant(X)

# 1. Pooled OLS
pooled_model = PanelOLS(y, X)
pooled_results = pooled_model.fit()

# 2. Efeitos Fixos
fe_model = PanelOLS(y, X, entity_effects=True)
fe_results = fe_model.fit()

# 3. Efeitos Aleatórios
re_model = RandomEffects(y, X)
re_results = re_model.fit()

# 4. Teste de Hausman
hausman_test = compare({'Efeitos Fixos': fe_results, 
                       'Efeitos Aleatórios': re_results})

print("Resumo dos Resultados:")
print(f"Pooled OLS - Alavancagem: {pooled_results.params['alavancagem']:.4f}")
print(f"Efeitos Fixos - Alavancagem: {fe_results.params['alavancagem']:.4f}")
print(f"Efeitos Aleatórios - Alavancagem: {re_results.params['alavancagem']:.4f}")
print(f"\nTeste de Hausman p-valor: {hausman_test.pval:.4f}")

# Implementação de DiD com tratamento escalonado
import pandas as pd
import numpy as np
from linearmodels import PanelOLS
import matplotlib.pyplot as plt

# Exemplo: Adoção de práticas ESG por empresas brasileiras
def did_staggered_adoption(df):
    """
    Implementa DiD com adoção escalonada usando a abordagem de 
    Callaway & Sant'Anna (2020)
    """
    # Criar variável de coorte (ano de adoção)
    df['cohort'] = df.groupby('empresa')['esg_adotado'].transform(
        lambda x: x.idxmax() if x.any() else np.nan
    )
    
    # Análise por coorte
    results = []
    cohorts = df['cohort'].dropna().unique()
    
    for cohort in cohorts:
        # Empresas tratadas nesta coorte
        treated = df[df['cohort'] == cohort]['empresa'].unique()
        
        # Empresas nunca tratadas (controle puro)
        never_treated = df[df['cohort'].isna()]['empresa'].unique()
        
        # Subset para análise
        subset = df[df['empresa'].isin(np.concatenate([treated, never_treated]))]
        
        # DiD para esta coorte
        subset['post'] = (subset.index.get_level_values('ano') >= cohort)
        subset['treated'] = subset['empresa'].isin(treated)
        subset['did'] = subset['post'] * subset['treated']
        
        # Estimação
        model = PanelOLS(subset['performance'], 
                        subset[['did', 'controles']], 
                        entity_effects=True, 
                        time_effects=True)
        
        result = model.fit()
        results.append({
            'cohort': cohort,
            'effect': result.params['did'],
            'se': result.std_errors['did'],
            'n_treated': len(treated)
        })
    
    return pd.DataFrame(results)

# Agregação dos efeitos por coorte
def aggregate_effects(cohort_results):
    """Agrega efeitos usando pesos por tamanho da coorte"""
    weights = cohort_results['n_treated'] / cohort_results['n_treated'].sum()
    
    avg_effect = np.sum(cohort_results['effect'] * weights)
    avg_se = np.sqrt(np.sum((cohort_results['se'] ** 2) * (weights ** 2)))
    
    return avg_effect, avg_se

# Event study para validar tendências paralelas
def event_study(df, treatment_var, outcome_var, time_to_treat_var):
    """Implementa event study para verificar pré-tendências"""
    
    # Criar dummies para cada período relativo ao tratamento
    periods = range(-5, 8)  # 5 períodos antes até 7 depois
    
    for period in periods:
        if period != -1:  # Omitir período -1 como referência
            df[f'period_{period}'] = (df[time_to_treat_var] == period).astype(int)
    
    # Modelo de event study
    period_vars = [f'period_{period}' for period in periods if period != -1]
    
    model = PanelOLS(df[outcome_var], 
                    df[period_vars + ['controles']], 
                    entity_effects=True, 
                    time_effects=True)
    
    results = model.fit()
    
    # Plotar resultados
    coefs = [0]  # período -1 (referência)
    ses = [0]
    
    for period in periods:
        if period != -1:
            coefs.append(results.params[f'period_{period}'])
            ses.append(results.std_errors[f'period_{period}'])
        else:
            coefs.append(0)
            ses.append(0)
    
    plt.figure(figsize=(10, 6))
    plt.plot(periods, coefs, 'o-', color='blue', label='Efeito Estimado')
    plt.fill_between(periods, 
                     np.array(coefs) - 1.96*np.array(ses),
                     np.array(coefs) + 1.96*np.array(ses),
                     alpha=0.2, color='blue', label='IC 95%')
    
    plt.axvline(x=-0.5, color='red', linestyle='--', label='Tratamento')
    plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    
    plt.xlabel('Períodos Relativos ao Tratamento')
    plt.ylabel('Efeito no Resultado')
    plt.title('Event Study: Teste de Tendências Paralelas')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    return results

# Implementação de painel dinâmico com GMM
import pandas as pd
import numpy as np
from linearmodels.iv import IVGMM
from linearmodels.panel import PanelData

def dynamic_panel_gmm(df, y_var, x_vars, entity_var, time_var):
    """
    Estima painel dinâmico usando GMM de Arellano-Bond
    """
    # Preparar dados
    panel_data = PanelData(df.set_index([entity_var, time_var]))
    
    # Criar variável dependente defasada
    y = panel_data.dataframe[y_var]
    y_lag = y.groupby(level=0).shift(1)
    
    # Diferenças primeiras
    y_diff = y.groupby(level=0).diff()
    y_lag_diff = y_lag.groupby(level=0).diff()
    
    # Preparar regressores
    X_diff = panel_data.dataframe[x_vars].groupby(level=0).diff()
    X_diff['y_lag'] = y_lag_diff
    
    # Instrumentos: valores defasados em níveis
    instruments = []
    for lag in range(2, 5):  # Usar lags 2-4 como instrumentos
        y_lag_instrument = y.groupby(level=0).shift(lag)
        instruments.append(y_lag_instrument)
    
    # Combinar instrumentos
    Z = pd.concat(instruments, axis=1)
    Z.columns = [f'y_lag{i}' for i in range(2, 5)]
    
    # Estimar GMM
    model = IVGMM(y_diff, X_diff, Z)
    results = model.fit(cov_type='robust')
    
    return results

# Teste de Sargan para validade dos instrumentos
def sargan_test(residuals, instruments, n_params):
    """
    Teste de sobreidentificação de Sargan
    H0: Instrumentos são válidos
    """
    n = len(residuals)
    Z = instruments.values
    u = residuals.values.reshape(-1, 1)
    
    # Estatística J
    J = n * (u.T @ Z @ np.linalg.inv(Z.T @ Z) @ Z.T @ u)[0, 0]
    
    # Graus de liberdade
    df = Z.shape[1] - n_params
    
    # P-valor
    from scipy.stats import chi2
    p_value = 1 - chi2.cdf(J, df)
    
    return {'J_stat': J, 'df': df, 'p_value': p_value}

# Análise completa de dados em painel - Setor Bancário Brasileiro
import pandas as pd
import numpy as np
from linearmodels import PanelOLS, RandomEffects
from linearmodels.panel import compare
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 1. Preparação e exploração dos dados
def prepare_bank_panel_data():
    """
    Simula dados de painel para bancos brasileiros
    """
    np.random.seed(42)
    
    # Parâmetros
    n_banks = 50
    n_periods = 40  # 10 anos trimestrais
    
    # Gerar identificadores
    banks = [f'Bank_{i}' for i in range(n_banks)]
    periods = pd.date_range('2014-01-01', periods=n_periods, freq='Q')
    
    # Criar painel
    panel_index = pd.MultiIndex.from_product([banks, periods], 
                                           names=['bank', 'quarter'])
    df = pd.DataFrame(index=panel_index).reset_index()
    
    # Características dos bancos (efeitos fixos)
    bank_effects = pd.DataFrame({
        'bank': banks,
        'efficiency': np.random.normal(0, 0.5, n_banks),
        'size_factor': np.random.uniform(0.5, 2.0, n_banks),
        'risk_profile': np.random.choice(['Conservative', 'Moderate', 'Aggressive'], 
                                      n_banks, p=[0.3, 0.5, 0.2])
    })
    
    # Merge com painel
    df = df.merge(bank_effects, on='bank')
    
    # Variáveis temporais
    df['time_trend'] = df.groupby('bank').cumcount()
    df['post_regulation'] = (df['quarter'] >= '2018-01-01').astype(int)
    
    # Variáveis explicativas
    df['size'] = np.log(df['size_factor'] * np.exp(df['time_trend'] * 0.02 + 
                                                  np.random.normal(0, 0.2, len(df))))
    df['capital_ratio'] = 0.12 + df['efficiency'] * 0.02 + \
                         np.random.normal(0, 0.03, len(df))
    df['loan_deposit_ratio'] = 0.7 + np.random.normal(0, 0.1, len(df))
    df['npl_ratio'] = 0.03 + (df['risk_profile'] == 'Aggressive') * 0.02 + \
                     np.random.normal(0, 0.01, len(df))
    
    # Variável dependente (ROA)
    df['roa'] = (0.015 + 
                0.5 * df['efficiency'] + 
                0.002 * df['size'] + 
                0.1 * df['capital_ratio'] - 
                0.5 * df['npl_ratio'] +
                df['post_regulation'] * (-0.003 + 0.005 * (df['size'] > df['size'].median())) +
                np.random.normal(0, 0.003, len(df)))
    
    return df

# 2. Análise exploratória
def exploratory_analysis(df):
    """
    Análise exploratória do painel
    """
    print("=== Estrutura do Painel ===")
    print(f"Número de bancos: {df['bank'].nunique()}")
    print(f"Número de períodos: {df['quarter'].nunique()}")
    print(f"Total de observações: {len(df)}")
    
    # Verificar balanceamento
    balance_check = df.groupby('bank')['quarter'].count()
    print(f"\nPainel balanceado: {balance_check.nunique() == 1}")
    
    # Estatísticas descritivas
    print("\n=== Estatísticas Descritivas ===")
    print(df[['roa', 'size', 'capital_ratio', 'npl_ratio']].describe())
    
    # Visualização da evolução temporal
    plt.figure(figsize=(12, 6))
    avg_roa = df.groupby('quarter')['roa'].mean()
    plt.plot(avg_roa.index, avg_roa.values, linewidth=2)
    plt.axvline(x=pd.Timestamp('2018-01-01'), color='red', linestyle='--', 
                label='Regulação')
    plt.xlabel('Trimestre')
    plt.ylabel('ROA Médio')
    plt.title('Evolução do ROA no Setor Bancário')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    # Matriz de correlação
    plt.figure(figsize=(10, 8))
    corr_matrix = df[['roa', 'size', 'capital_ratio', 'loan_deposit_ratio', 'npl_ratio']].corr()
    sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0)
    plt.title('Matriz de Correlação')
    plt.show()

# 3. Estimação de modelos
def estimate_panel_models(df):
    """
    Estima diferentes especificações de modelos de painel
    """
    # Preparar dados para linearmodels
    df_panel = df.set_index(['bank', 'quarter'])
    
    # Variáveis
    y = df_panel['roa']
    X_base = df_panel[['size', 'capital_ratio', 'npl_ratio', 'post_regulation']]
    
    # Modelo 1: Pooled OLS
    pooled_model = PanelOLS(y, X_base)
    pooled_results = pooled_model.fit()
    
    # Modelo 2: Efeitos Fixos (entidade)
    fe_model = PanelOLS(y, X_base, entity_effects=True)
    fe_results = fe_model.fit()
    
    # Modelo 3: Efeitos Fixos Two-way
    twfe_model = PanelOLS(y, X_base, entity_effects=True, time_effects=True)
    twfe_results = twfe_model.fit()
    
    # Modelo 4: Efeitos Aleatórios
    re_model = RandomEffects(y, X_base)
    re_results = re_model.fit()
    
    # Teste de Hausman
    hausman = compare(fe_results, re_results)
    print(f"Teste de Hausman - Chi2: {hausman.stat:.4f}, p-valor: {hausman.pval:.4f}")
    
    # Comparar resultados
    results_df = pd.DataFrame({
        'Pooled OLS': pooled_results.params,
        'FE': fe_results.params,
        'TWFE': twfe_results.params,
        'RE': re_results.params
    })
    
    print("\n=== Comparação de Coeficientes ===")
    print(results_df)
    
    return pooled_results, fe_results, twfe_results, re_results

# 4. Análise de heterogeneidade
def heterogeneity_analysis(df):
    """
    Analisa efeitos heterogêneos da regulação
    """
    df_panel = df.set_index(['bank', 'quarter'])
    
    # Criar interações
    df_panel['large_bank'] = (df_panel['size'] > df_panel['size'].median()).astype(int)
    df_panel['reg_x_large'] = df_panel['post_regulation'] * df_panel['large_bank']
    
    # Modelo com interação
    y = df_panel['roa']
    X = df_panel[['size', 'capital_ratio', 'npl_ratio', 
                  'post_regulation', 'large_bank', 'reg_x_large']]
    
    het_model = PanelOLS(y, X, entity_effects=True, time_effects=True)
    het_results = het_model.fit()
    
    print("\n=== Análise de Heterogeneidade ===")
    print(het_results.summary.tables[1])
    
    # Visualizar efeitos heterogêneos
    effect_small = het_results.params['post_regulation']
    effect_large = het_results.params['post_regulation'] + het_results.params['reg_x_large']
    
    plt.figure(figsize=(10, 6))
    plt.bar(['Bancos Pequenos', 'Bancos Grandes'], 
            [effect_small, effect_large], 
            color=['#3498db', '#e74c3c'])
    plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    plt.ylabel('Efeito sobre ROA')
    plt.title('Efeito Heterogêneo da Regulação por Tamanho do Banco')
    plt.grid(True, alpha=0.3, axis='y')
    plt.show()
    
    return het_results

# 5. Diagnósticos e testes de robustez
def diagnostic_tests(df, model_results):
    """
    Realiza testes diagnósticos
    """
    # Teste de autocorrelação
    residuals = model_results.resids
    
    # Teste de Wooldridge para autocorrelação em painel
    df_res = pd.DataFrame(residuals).reset_index()
    df_res.columns = ['bank', 'quarter', 'residual']
    df_res['lag_residual'] = df_res.groupby('bank')['residual'].shift(1)
    
    # Regressão dos resíduos no lag
    from sklearn.linear_model import LinearRegression
    
    mask = ~df_res['lag_residual'].isna()
    X = df_res.loc[mask, ['lag_residual']].values
    y = df_res.loc[mask, 'residual'].values
    
    autoreg = LinearRegression().fit(X, y)
    rho = autoreg.coef_[0]
    
    print(f"\nCoeficiente de autocorrelação (ρ): {rho:.4f}")
    
    # Teste de heterocedasticidade
    # Breusch-Pagan simplificado
    squared_resids = residuals ** 2
    het_model = PanelOLS(squared_resids, model_results.model.exog)
    het_results = het_model.fit()
    
    n = len(residuals)
    lm_stat = n * het_results.rsquared
    p_value = 1 - stats.chi2.cdf(lm_stat, het_results.model.exog.shape[1] - 1)
    
    print(f"Teste de Heterocedasticidade - LM: {lm_stat:.4f}, p-valor: {p_value:.4f}")

# 6. Pipeline completo
def run_complete_analysis():
    """
    Executa análise completa de painel
    """
    # Gerar dados
    df = prepare_bank_panel_data()
    
    # Análise exploratória
    exploratory_analysis(df)
    
    # Estimação de modelos
    pooled, fe, twfe, re = estimate_panel_models(df)
    
    # Análise de heterogeneidade
    het_results = heterogeneity_analysis(df)
    
    # Diagnósticos
    diagnostic_tests(df, twfe)
    
    # Event study
    event_study_analysis(df)
    
    return df, twfe, het_results

# 7. Event study para validar identificação
def event_study_analysis(df):
    """
    Implementa event study para verificar tendências paralelas
    """
    df_panel = df.set_index(['bank', 'quarter'])
    
    # Criar variáveis de tempo relativo à regulação
    df_panel['time_to_treat'] = df_panel.index.get_level_values('quarter').map(
        lambda x: (x - pd.Timestamp('2018-01-01')).days // 90
    )
    
    # Criar dummies para cada período
    periods = range(-8, 8)
    for period in periods:
        if period != -1:  # Omitir período -1 como base
            df_panel[f'period_{period}'] = (df_panel['time_to_treat'] == period).astype(int)
    
    # Estimar modelo
    y = df_panel['roa']
    period_vars = [f'period_{p}' for p in periods if p != -1]
    controls = ['size', 'capital_ratio', 'npl_ratio']
    X = df_panel[period_vars + controls]
    
    event_model = PanelOLS(y, X, entity_effects=True)
    event_results = event_model.fit()
    
    # Plotar resultados
    coefs = []
    ses = []
    
    for period in periods:
        if period == -1:
            coefs.append(0)
            ses.append(0)
        else:
            coefs.append(event_results.params[f'period_{period}'])
            ses.append(event_results.std_errors[f'period_{period}'])
    
    plt.figure(figsize=(12, 8))
    plt.plot(periods, coefs, 'o-', color='blue', linewidth=2, markersize=8)
    plt.fill_between(periods,
                     np.array(coefs) - 1.96*np.array(ses),
                     np.array(coefs) + 1.96*np.array(ses),
                     alpha=0.2, color='blue')
    
    plt.axvline(x=-0.5, color='red', linestyle='--', linewidth=2, label='Regulação')
    plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)
    
    plt.xlabel('Trimestres Relativos à Regulação')
    plt.ylabel('Efeito sobre ROA')
    plt.title('Event Study: Impacto da Regulação Bancária')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    return event_results

# Executar análise
if __name__ == "__main__":
    df, twfe_results, het_results = run_complete_analysis()