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

Explorando a base de dados do paper "The Mission: Human Capital Transmission, Economic Persistence, and Culture in South America"

0 votos
13 visitas
perguntada Mai 19 em Economia por gustavobrangel (6 pontos)  
editado Mai 19 por gustavobrangel

Referência e Motivação

A literatura de persistência econômica sugere que algumas características da sociedade moderna são reflexos de eventos naturais ou quase-naturais ocorridos em um passado distante. Por exemplo, o trabalho mais conhecido dessa literatura é o "The Colonial Origins of Comparative Development: An Empirical Investigation" de Acemoglu, Johnson e Robinson. Os autores argumetam que parte da diferença atual observada entre as instituições de antigas colônias europeias pode ser atribuída à diferença na taxa de mortalidade dos europeus no período colonial. Onde os europeus encontraram melhores condições de sobrevivência (América do Norte e Austrália) foi também onde eles formaram instituições com mais qualidade e esse efeito persiste até os dias de hoje.

Esse é um tipo de literatura bem interessante principalmente por lidar com dados muito antigos. Seus pesquisadores fazem um excelente trabalho de coleta de dados e, pensando nisso, resolvi escolher um artigo dessa área para pode explorar sua base dados. Já a escolha do paper especificamente foi por esse ser um dos poucos artigos nessa literatura tratando de dados brasileiros.

A principal hipótese do artigo é que os lugares que hoje se encontram mais próximos das missões jesuítas do período colonial são também os lugares com mais capital humano. A razão para tanto é que as missões jesuítas não se preocupavam apenas com a pregação religiosa. Os missionários também ensinavam os indígenas a ler, escrever, realizar operações aritméticas básicas e a praticar alguns ofícios. Assim, o autor argumenta que esse desenvolvimento de capital humano persistiu e seus efeitos podem ser observados até hoje a depender da proximidade dos lugares atuais em relação à missão jesuíta mais próxima.

Detalhamento da Base de Dados

Segue a tabela com as principais estatísticas descritivas da base de dados:

A imagem será apresentada aqui.

Resultados econométricos

Os principais resultados do artigo estão contidos na Tabela II, que pode ser visualizada abaixo:

A imagem será apresentada aqui.

Ao ter acesso à base de dados, tento replicar esses resultados, e chego aos seguintes valores:

A imagem será apresentada aqui.

Como pode ser observado, não consigo chegar aos valores apresentado no artigo para os coeficientes e desvios-padrão da regressão que contempla apenas o Brasil. Ao ler atentamente as notas abaixo da tabela, percebe-se o autor inclui efeitos fixos de mesorregiões apenas para o caso do Brasil. Então, resolvi iniciar a investigação na diferença desses resultados por aí.

Como já havia feito a codificação da variável categória \textit{mesorregi} por meio de variáveis dummies (também conhecido como one-hot encoding, resolvo, então, utilizar tal codificação por meio de apenas uma variável (também conhecido como ordinal encoding). Ao fazer isso, chego exatamente aos valores que o autor do paper publicou.

Entretanto, será que essa é realmente a forma correta de se realizar o ajuste por efeitos fixos? Para testar isso, altero a ordem de representação das mesorregiões na codificação ordinal e também ajusto uma regressão com esses novos valores. Vamos comparar, então, esses três tipos de representação da variável categórica de mesorregiões:

A imagem será apresentada aqui.

A imagem será apresentada aqui.

A imagem será apresentada aqui.

Representar a mesorregião em que o município de, por exemplo, Aceguá se encontra de 0 ou de 5 não deveria alterar o resultado. Caso altere, é porque tem algo errado com esse tipo de codificação. Utilizo a estátistica \(t\) para comparar os três diferentes resultados:

A imagem será apresentada aqui.

Observa-se que, de fato, o ordinal encoding não parece ser a técnica ideal para ajuste por efeitos fixos.

Além disso, quando não há o ajuste por variáveis geográficas, a estatística \(t\) cai para abaixo de 2, ou seja, a distância para a missão jesuíta mais próxima deixa de ser estatisticamente significante na determinação do capital humano. A abordagem clássica possui essa desvantagem: as conclusões são binárias. Existe ou não existe efeito. E tudo depende do p-valor resultante da estatística \(t\). E, por acaso, existe alguma outra abordagem de estimação que não faça uso do p-valor?

A resposta é sim: um paradigma do ramo da estatística iniciado por Thomas Bayes. Na estimação bayesiana, utilizam-se algoritmos MCMC (Markov Chain Monte Carlo) para se conhecer melhor o parâmetro de interesse, como uma alternativa ao MQO clássico. Aqui, eu foco apenas na apresentação do resultado final desse tipo de estimação:

A imagem será apresentada aqui.

Percebe-se que, no lugar de apenas o valor do coeficiente e seu desvio-padrão (como se costuma apresentar os resultados na abordagem clássica), o paradigma bayesiano traz uma informação mais completa a respeito do parâmetro de interesse. Pela figura, podemos concluir que a média do coeficiente é 0.016, mas que há 94% de chance de esse coeficiente estar entre 0.003 e 0.03.

A abordagem bayesiana é uma técnica bem diferente em relação à frequentista. Ambas têm suas desvantagens e desvantagens. A bayesiana não é melhor, nem pior, e um bom pesquisador deve conhecer bem ambas para saber identificar quando utilizar cada uma delas.

Compartilhe
comentou Mai 19 por gustavobrangel (6 pontos)  

Código utilizado na análise do ajuste por efeitos fixos:

import pandas

import seaborn

from pathlib import Path

from statsmodels.regression.linear_model import OLS

white_colour = (250/255, 250/255, 250/255)

pandas.set_option('display.max_colum', 100)
seaborn.set_context('notebook')
seaborn.set_style('whitegrid', {'axes.facecolor': white_colour,
'figure.facecolor': white_colour})

images_path = Path('my_path_to_images')

table2 = pandas.read_stata('my_path_to_database')

br_data = table2[table2.country == 'BRA']

br_data = br_data.astype({'mesorregi': int})

br_data['const'] = 1

br_data['MR_seq1'] = pandas.factorize(br_data.mesorregi, sort=True)[0]
br_data['MR_seq2'] = pandas.factorize(br_data.mesorregi, sort=False)[0]

br_data = pandas.get_dummies(br_data, columns=['mesorregi'], drop_first=True, prefix='MR', prefix_sep=':')

mr_cols = [col for col in br_data.columns if 'MR:' in col]

y = ['illiteracy']
X_col_ordinal1 = 'distmiss lati longi MR_seq1 const'.split()
X_col_ordinal2 = 'distmiss lati longi MR_seq2 const'.split()
X_col_onehot = 'distmiss lati longi const'.split() + mr_cols
geo_cols = 'area tempe alti preci rugg river coast'.split()

ols_onehot = OLS(br_data[y], br_data[X_col_onehot], missing='drop').fit(cov_type='HC1')
ols_onehot_geo = OLS(br_data[y], br_data[X_col_onehot + geo_cols], missing='drop').fit(cov_type ='HC1')

ols_ordinal1 = OLS(br_data[y], br_data[X_col_ordinal1], missing='drop').fit(cov_type='HC1')
ols_ordinal1_geo = OLS(br_data[y], br_data[X_col_ordinal1 + geo_cols], missing='drop').fit(cov_type ='HC1')

ols_ordinal2 = OLS(br_data[y], br_data[X_col_ordinal2], missing='drop').fit(cov_type ='HC1')
ols_ordinal2_geo = OLS(br_data[y], br_data[X_col_ordinal2 + geo_cols], missing='drop').fit(cov_type ='HC1')

ols_models = [ols_ordinal1_geo, ols_ordinal2_geo, ols_onehot_geo, ols_ordinal1, ols_ordinal2, ols_onehot]

coef = []
sd = []

for model in ols_models:
coef.append(model.params['distmiss'])
sd.append(model.bse['distmiss'])

results_br_table2 = pandas.DataFrame({
'coef': coef,
'sd': sd,
'encoding': ['ordinal_seq1', 'ordinal_seq2', 'one-hot'] * 2,
'geo_controls': ['yes', 'no'] * 3
})

results_br_table2['t score'] = abs(results_br_table2.coef / results_br_table2.sd)

fig_ordinal_only = seaborn.catplot(data=results_br_table2, x='encoding', y='t score', col='geo_controls', kind='bar',
sharey=False, palette=['tab:blue', 'tab:orange', white_colour])
fig_ordinal_only.set_titles(col_template='Geo controls: {col_name}')
fig_ordinal_only.set_xticklabels(['ordinal_seq1', 'ordinal_seq2', ''])
fig_ordinal_only.set_axis_labels("", "t score")
fig_ordinal_only.tight_layout()
fig_ordinal_only.savefig(images_path / 'encoding_ordinal_only.png')

fig_complete = seaborn.catplot(data=results_br_table2, x='encoding', y='t score', col='geo_controls', kind='bar',
sharey=False, palette=['tab:blue', 'tab:orange', 'tab:green'])
fig_complete.set_titles(col_template='Geo controls: {col_name}')
fig_complete.set_axis_labels("", "t score")
fig_complete.tight_layout()
fig_complete.savefig(images_path / 'encoding_complete.png')

comentou Mai 19 por gustavobrangel (6 pontos)  
Código utilizado na regressão bayesiana:

import pandas

import pymc3

import seaborn

from pathlib import Path

from matplotlib import pyplot

pandas.set_option('display.max_colum', 100)
seaborn.set_context('notebook')
seaborn.set_style('white')
seaborn.set(rc={'axes.facecolor': (250/255, 250/255, 250/255), 'figure.facecolor': (250/255, 250/255, 250/255)})

__DBPATH__ = Path('my_path_to_db')

__FIGPATH__ = Path('my_path_to_figure')

table2 = pandas.read_stata(__DBPATH__ / 'Literacy Argentina Brazil Paraguay.dta')

br_data = table2[table2.country == 'BRA']

# br_data.loc[:, 'illiteracy'] = br_data.loc[:, 'illiteracy'] / 100

br_data = br_data[['illiteracy', 'distmiss', 'lati', 'longi', 'mesorregi']].dropna()

br_data['mesorregi'] = br_data['mesorregi'].astype(int)

br_data['mesorregi'] = br_data.mesorregi - min(br_data.mesorregi)

br_data = pandas.get_dummies(br_data, columns=['mesorregi'], drop_first=True, prefix='MR', prefix_sep='_')

mr_cols = [col for col in br_data.columns if 'MR_' in col]

regressors = ['distmiss', 'lati', 'longi'] + mr_cols

all_vars = ['illiteracy'] + regressors

formula = 'illiteracy ~ ' + ' + '.join(regressors)

draws = 10_000
tune = 3_000
target_accept = 0.9

with pymc3.Model() as model:
    priors = {
        "illiteracy": pymc3.Beta.dist(alpha=1, beta=1),
    }
    pymc3.GLM.from_formula(formula, data=br_data, priors=priors)
    trace = pymc3.sample(draws=draws, tune=tune, target_accept=target_accept)

    results = pandas.read_csv(__FIGPATH__ / 'results_info.csv')

    if not len(results):
        index = 0
    else:
        index = results.index[-1] + 1

    fig_posterior, ax = pyplot.subplots()

    pymc3.plot_posterior(trace['distmiss'], ax=ax)

    ax.set_title(r'Posterior density of $\beta$')

    fig_posterior.tight_layout()

    fig_posterior.savefig(__FIGPATH__ / f'{index:02}_posterior_plot')

    pymc3.traceplot(trace)

    fig_traceplot = pyplot.gcf()  # get the current figure

    # fig_traceplot.tight_layout()

    fig_traceplot.savefig(__FIGPATH__ / f'{index:02}_trace_plot')

    # pymc3.forestplot(trace)

    # pymc3.autocorrplot(trace)

    # pymc3.energyplot(trace)

    diverging = trace['diverging']
    diverging_pct = diverging.nonzero()[0].size / len(trace) * 100

    with open(__FIGPATH__ / f'{index:02}_trace_summary', 'w') as f:
        print(f'Number of Divergent Chains: {diverging.nonzero()[0].size}', file=f)
        print(f'Percentage of Divergent Chains: {diverging_pct:.1f}\n', file=f)
        print(pymc3.summary(trace), file=f)

    info_dict = {
        'results_index': index,
        'draws': len(trace),
        'tune': tune,
        'chains': len(trace.chains),
        'target_accept': target_accept
    }

    priors_info = {key: priors.pop(key, 'default') for key in all_vars}

    priors_info = {key: priors_info[key] if value == 'default'
                   else value.__str__().split('object')[0].split('.')[-1].strip()
                   for key, value in priors_info.items()
                   }

    priors_info['MR_0'] = 'drop'

    info_dict = {**info_dict, **priors_info}

    results_update = results.append(info_dict, ignore_index=True)

    results_update.to_csv(__FIGPATH__ / 'results_info.csv', index=False)

Entre ou cadastre-se para responder esta pergunta.

...