Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

É fácil rodar modelos de painel de dados (dados longitudinais) em Python?

0 votos
417 visitas
perguntada Mar 5, 2016 em Estatística por danielcajueiro (5,251 pontos)  
Compartilhe

1 Resposta

+2 votos
respondida Mar 5, 2016 por danielcajueiro (5,251 pontos)  

Essa é uma pergunta muito ampla, pois existem muitos modelos usados para esse fim. Enquanto alguns deles não são trivialmente implementados, outros estão disponíveis no pacote pandas. Entretanto, se considerarmos os modelos mais comuns, eles podem ser facilmente rodados em Python, mesmo que os pacotes específicos estejam ainda em fase de desenvolvimento e provavelmente serão adicionados ao pacote statsmodels (que é uma implementação de modelos econométricos que copia o estilo do R).

Que tal um exemplo?

Vamos usar a amostra de dados "RENTAL.raw" que é bem conhecida e discutida no livro "Introductory Econometrics - Jeffrey Wooldridge", que pode ser baixada no site da Cengage.

Essa base foi coletada nos anos 1980 e 1990 e inclui preços de aluguel de imóveis em cidades universitárias (cidades cuja a população é dominada pela população universitária). A pergunta de interesse é se a forte presença de estudantes afeta os preços dos imóveis.

Usando o modelo de Painel de dados conhecido como Efeitos Fixos (within estimator), vamos estimar o modelo

\(log(rent_{it})=\beta_1 y90_{t} + \beta_2 log(pop_{it}) +\beta_3 log(avinc_{it})\) \( + \beta_4 pctstu_{it} + a_i + u_{it},\)

onde \(rent\) é o preço médio do aluguel, \(y90\) é uma variável dummy para o ano 1990, \(pop\) é a população da cidade, \(avinc\) é a renda per capita e \(pctstu\) é a porcentagem de estudantes na cidade.

Os modelos de painel de dados de efeitos fixos supõem que desejamos estimar um modelo do tipo

\[y_{it}=\beta_1 x_{it} + a_i + u_{it}, \; (1)\]

onde \(a_i\) é um efeito associado a cada indivíduo que não é observável. Uma forma simples de estimar esse modelo é tomar a média desse modelo dada por

\[\overline y_{i}=\beta_1 \overline x_{i} + a_i + \overline u_{i}\; (2),\]

onde \(\overline y_{i}=\sum_{t=1}^{T}y_{it}\).

Subtraindo (2) de (1) chegamos a

\[y_{it}-\overline y_{i}=\beta_1 (x_{it}-\overline x_{i}) + (u_{it}-\overline u_{i})\]

e podemos usar diretamente OLS para estimar o vetor de parâmetros \(\beta_1\).

Em Python, começamos lendo a base de dados

df = pd.read_csv('RENTAL.raw',delim_whitespace=True,header=None,na_values=".")
df.columns=['city','year','pop','enroll','rent','rnthsg','tothsg','avginc',\
'lenroll','lpop','lrent','ltothsg','lrnthsg','lavginc','clenroll','clpop',\
'clrent','cltothsg','clrnthsg','clavginc','pctstu','cpctstu','y90']

df=df[['year','city','lrent','y90','lpop','lavginc','pctstu']]    
df = df.dropna() # Drop missing        

Depois precisamos criar as variáveis acima que foram subtraídas de sua média:

df[['meanlrent','meany90','meanlpop','meanlavginc','meanpctstu']] = df.groupby('city')[['lrent','y90','lpop','lavginc','pctstu']].transform('mean')

df['difflrent']=df['lrent']-df['meanlrent']
df['diffy90']=df['y90']-df['meany90']
df['difflpop']=df['lpop']-df['meanlpop']
df['difflavginc']=df['lavginc']-df['meanlavginc']
df['diffpctstu']=df['pctstu']-df['meanpctstu']    

Finalmente, usamos OLS para estimar o modelo como já mostramos aqui:

y,X = ps.dmatrices('difflrent ~ diffy90 + difflpop + difflavginc + diffpctstu -1',data=df, return_type='dataframe')
modelFE = sm.OLS(y,X) # Describe Model
resultsFE = modelFE.fit() # Fit model
print resultsFE.summary()

Chegamos então ao seguinte sumário:

                            OLS Regression Results                            
==============================================================================
Dep. Variable:              difflrent   R-squared:                       0.977
Model:                            OLS   Adj. R-squared:                  0.976
Method:                 Least Squares   F-statistic:                     1290.
Date:                Sat, 05 Mar 2016   Prob (F-statistic):          5.75e-100
Time:                        08:54:11   Log-Likelihood:                 219.27
No. Observations:                 128   AIC:                            -430.5
Df Residuals:                     124   BIC:                            -419.1
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
===============================================================================
                  coef    std err          t      P>|t|      [95.0% Conf. Int.]
-------------------------------------------------------------------------------
diffy90         0.3855      0.026     15.050      0.000         0.335     0.436
difflpop        0.0722      0.061      1.176      0.242        -0.049     0.194
difflavginc     0.3100      0.046      6.703      0.000         0.218     0.401
diffpctstu      0.0112      0.003      3.898      0.000         0.006     0.017
==============================================================================
Omnibus:                        0.015   Durbin-Watson:                   3.159
Prob(Omnibus):                  0.992   Jarque-Bera (JB):                0.023
Skew:                          -0.000   Prob(JB):                        0.989
Kurtosis:                       2.934   Cond. No.                         23.0
==============================================================================

Falta algo? De fato, sim! Note que para calcularmos modelo sem média acima reduzimos os graus de liberdade do modelo. Enquanto o modelo calculado pelo OLS convencional tem \(NT-K\) graus de liberdade, o nosso modelo tem \(NT-N-K\) graus de liberdade (pois calculamos \(N\) médias adicionais). Logo, precisamos corrigir os desvios padrões, as estatísticas \(t\) e os \(p\)-valores. É muito simples de fazer isso em Python e incluímos essa função para fazer isso:

def FixedEffectStats(coefficient,results,significanceLevel,N,T,):
    df=results.df_model
    coefficientValue=results.params[coefficient]
    degreesOfFreedom=N*T-N-df
    standardErrorEstimator=results.bse[coefficient]*(np.sqrt(N*T-df)/np.sqrt(N*T-N-df))

    tstat=coefficientValue/standardErrorEstimator

    pvalue=2*(1-scipy.stats.t.cdf(np.abs(tstat),degreesOfFreedom))

    print coefficient,' ',"{0:.4f}".format(coefficientValue),' ',  "{0:.4f}".format(standardErrorEstimator),' ',"{0:.4f}".format(tstat) ,' ',"{0:.4f}".format(pvalue)

Essa função basicamente corrige o desvio padrão do estimador multiplicando aquele calculado pelo OLS por \(\sqrt(N*T-K)/\sqrt(N*T-N-K)\) e recalcula as novas estatísticas \(t\) e \(p\)-valores considerando essa correção (o programa completo é apresentado abaixo) que nos dá os seguintes resultados corrigidos:

          coef     std err       t       pvalue
diffy90   0.3855   0.0368   10.4692   0.0000
difflpop   0.0722   0.0883   0.8178   0.4167
difflavginc   0.3100   0.0665   4.6627   0.0000
diffpctstu   0.0112   0.0041   2.7114   0.0087

De acordo com esse modelo podemos concluir que aumentando a porcentagem de estudantes na cidade, também aumenta-se o preço dos alugueis.

Se você tem interesse em Python e Econometria, que tal dar uma olhada aqui:

Como o python pode ser usado para lidar com dados e, particularmente, organização de dados?

Eu sei que com o Python, você pode fazer praticamente tudo... Mas uma pessoa que não tem muito traquejo com computação, consegue rodar modelos de regressões em Python?

Exemplos do livro Introductory Econometrics - Wooldridge em Python!

É possível trabalhar com análise de séries temporais usando Python, mesmo sem dominar muito computação?

É possível usar o python como framework para trabalhar com problemas em estatística?

O código completo:

import numpy as np
import statsmodels.api as sm
import pandas as pd
import patsy as ps
import scipy


def FixedEffectStats(coefficient,results,significanceLevel,N,T,):
    df=results.df_model
    coefficientValue=results.params[coefficient]
    degreesOfFreedom=N*T-N-df
    standardErrorEstimator=results.bse[coefficient]*(np.sqrt(N*T-df)/np.sqrt(N*T-N-df))

    tstat=coefficientValue/standardErrorEstimator

    pvalue=2*(1-scipy.stats.t.cdf(np.abs(tstat),degreesOfFreedom))

    print coefficient,' ',"{0:.4f}".format(coefficientValue),' ',  "{0:.4f}".format(standardErrorEstimator),' ',"{0:.4f}".format(tstat) ,' ',"{0:.4f}".format(pvalue)

if __name__ == '__main__':

    df = pd.read_csv('/home/daniel/Documentos/Projetos/Prorum/Python For Econometrics/DataSets/Txt/RENTAL.raw',delim_whitespace=True,header=None,na_values=".")
    df.columns=['city','year','pop','enroll','rent','rnthsg','tothsg','avginc',\
    'lenroll','lpop','lrent','ltothsg','lrnthsg','lavginc','clenroll','clpop',\
    'clrent','cltothsg','clrnthsg','clavginc','pctstu','cpctstu','y90']

    df=df[['year','city','lrent','y90','lpop','lavginc','pctstu']]    
    df = df.dropna() # Drop missing        

    # Create the demeaned data

    df[['meanlrent','meany90','meanlpop','meanlavginc','meanpctstu']] = df.groupby('city')[['lrent','y90','lpop','lavginc','pctstu']].transform('mean')

    df['difflrent']=df['lrent']-df['meanlrent']
    df['diffy90']=df['y90']-df['meany90']
    df['difflpop']=df['lpop']-df['meanlpop']
    df['difflavginc']=df['lavginc']-df['meanlavginc']
    df['diffpctstu']=df['pctstu']-df['meanpctstu']    

    y,X = ps.dmatrices('difflrent ~ diffy90 + difflpop + difflavginc + diffpctstu -1',data=df, return_type='dataframe')
    modelFE = sm.OLS(y,X) # Describe Model
    resultsFE = modelFE.fit() # Fit model
    print resultsFE.summary()

    print "=================================================================="
    print "          "+"coef"+"     "+"std err"+"     "+"  t  "+"     "+"pvalue"

    listOfParameters=['diffy90','difflpop','difflavginc','diffpctstu']
    for param in listOfParameters:
        FixedEffectStats(param,resultsFE,0.05,64,2)
...