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

Explorando a base de dados "What Causes Haze Pollution? An Empirical Study of PM2.5 Concentrations in Chinese Cities"

+1 voto
21 visitas
perguntada Mai 24 em Economia por Rodrigo Fernandes (61 pontos)  

"What Causes Haze Pollution? An Empirical Study of PM2.5 Concentrations in Chinese Cities" de Wu et al, Sustainability 2016, 8(2), 132; https://doi.org/10.3390/su8020132

Pelo journal ser aberto, é possível baixar diretamente de:
https://www.mdpi.com/2071-1050/8/2/132

Compartilhe

2 Respostas

+1 voto
respondida Mai 24 por Rodrigo Fernandes (61 pontos)  
 
Melhor resposta

Motivação

China afetada por hazes (ou seja, neblinas de poluição) e um grande causador desse evento são as partículas pequenas (PM 2.5) que causam grande impacto pra saúde, devido sua capacidade de afetar o pulmão por ser diminuta, além da poluição visual. Objetivo dos autores é entender, então, o porquê de algumas regiões terem maior concentração de PM 2.5 que outras. Para isso, observou principalmente fatores socioeconômicos, controlando para variáveis ambientais, que serão mostrados mais a frente.

A importância do artigo é clara se levando em consideração os caminhos soturnos que a humanidade traça com relação a preservação ambiental, seu impacto na saúde humana e segurança alimentar/habitacional, por exemplo. Assim, para os anos de 2013 e 2014, os autores tentam achar os motivos que determinam o grau de concentração de PM 2.5 na atmosfera, baseado nas observações de 74 cidades.

Modelo, dados e metodologia do artigo

O modelo utilizado pelos autores é o que se segue

\(P M_{2.5}(i, t)=\beta_{1} \times\) Industry \(_{i, t-1}+\beta_{2} \times\) Housing \(_{i, t-1}+\beta_{3} \times\) Vehicle \(_{i, t-1}+\beta_{4} \times H_{-}\) Consumption \(_{i, t-1}+\)
\(\beta_{5} \times\) Expenditure \(_{i, t-1}+\beta_{5} \times\) Control \(_{i, t-1}+\beta_{6} \times\) Year \(+\varepsilon\)

Onde i e t são os subscritos de localidade e tempo respectivamente, Industry é o grau de atividade industrial na cidade, Housing é a quantidade de obras na localidade, Vehicle o número de veículos, Consumption é o gasto com gás de cozinha (gás natural), Expenditure é o gasto do governo com preservação ambiental e Control são as variáveis de controle, dentre elas variáveis ambientais e de PIB. Por último, os autores adicionam uma dummy temporal.

Os autores fazem, então, 3 regressões. Primeiro uma OLS simples por ano, só para comparação. Depois, fazem um teste de Hausman para determinar se devem fazer uma regressão com efeito fixo ou random effects e acham que a segunda é melhor. Dessa forma, acabam por optar pela mesma.

Os autores acham que concentração de indústria parece ter o maior efeito em PM 2.5. Paralelamente, gás de cozinha e quantidade de veículos parecem ter efeito mais dissipado. Ainda, os autores acha que os gastos com conservação parecem não afetar a concentração de PM2.5 e considera as políticas públicas ineficientes. Para resumir, a China parece ter um problema de alto crescimento baseado na alta poluição e é essa premissa que nós buscamos testar para o caso brasileiro.

Modelo para nosso caso e diferenças principais

Alto crescimento baseado na alta poluição? Tentei adaptar isso para o caso brasileiro. Existe grande dificuldade de pegar dados do estado chinês, principalmente por grande parte dos dados estar em mandarim. Assim, resolvi por modificar a variável dependente para monóxido de carbono (CO) pela facilidade de se encontrar. Ainda, é importante levar em consideração as diferenças no corpo humano e no meio ambiente entre os tipos de poluição. O trabalho foi feito com dados entre 2009 e 2017 para todos os municípios brasileiros.

Estatística descritiva
A imagem será apresentada aqui.

Modelo
O modelo usado foi o abaixo

\(CO(i, t)=\beta_{1} \times\) Gado \(_{i, t-1}+\beta_{2} \times\) Carros \(_{i, t-1}+\beta_{3} \times\) LogPib\(_{-}\)per\(_{-}\)capita\(_{i, t-1} +\) \(\beta_{4} \times\) Controle\(_{-}\)Ambiental\(_{i, t-1}+\beta_{5} \times\) Ano\( + \varepsilon\)

Onde i e t são os subscritos de localidade e tempo novamente, Gado é o número de gado por município, o mesmo para Carros. Já logpib_ per capita é auto explicativo e ControleAmbiental são as variáveis de controle para meio ambiente. São três: temperatura, precipitação e umidade. A última foi retirada após apresentar um alto grau de correlação observado com um teste de variance inflation factor (VIF). Há ainda, para finalizar, dummy de ano.

É importante mencionar que monóxido de carbono não é um produto direto do gado como é o dióxido de carbono, por exemplo. Não obstante, por ser um setor econômico altamente poluente, ainda é interessante tentar observar seu impacto na concentração.

Regressões

Foram então, realizadas três regressões. A primeira foi uma OLS simples, como os autores. Após isso, realizei o Teste de Hausman para definir se seria melhor efeito fixo ou random effects. Com isso, rejeitei hipótese nula de random effects melhor, ou seja, recomendou-se efeito fixo para tempo e município.

Resultados

As três regressões são apresentadas abaixo onde os pontos de interesse estão destacadas. Para todas utilizados constantes e cluster para nível estadual com objetivo de captar parte da heterogeneidade.
A imagem será apresentada aqui.

De forma geral, observamos principalmente impactos de gados e número de veículos no grau de concentração de monóxido de carbono. Não obstante, quando usamos efeito fixo (o recomendado de acordo com o teste de Hausman), o efeito de gado não é observado. Ainda, o PIB não necessariamente está ligado com a quantidade de emissões de CO (o que demonstra uma possibilidade de um crescimento sustentável, pelo menos no que se refere a concentração desse gás específico) mas é significante quando analisamos OLS e seu p-valor se aproxima do crítico (5%) no Random Effects e é significante para 10%.

Alguns problemas podem ser levantados como a pouca disponibilidade de dados do Estado/oficiais, assim a maioria são de organizações (no que se refere ao meio ambiente). Ainda, alguns dados que não tive tempo de encontrar como gasto com conservação ambiental e PIB industrial.

0 votos
respondida Mai 24 por Rodrigo Fernandes (61 pontos)  

Código do python

import os
import pandas as pd
import pickle
path=os.getcwd()
os.chdir()
files = os.listdir(path)
files = [f for f in files if  f.lower().endswith(".dta")]

var = dict()

for file in files: 
    var[file] = pd.read_stata(file)    
a2008 = var['2008.dta']
filename = 'mypickle.pk'

with open(filename, 'wb') as fi:
    pickle.dump(var, fi)

for key in var:
    print(key)


import pickle
filename = 'mypickle.pk'
with open(filename, 'rb') as fi:
    var = pickle.load(fi)

ambiente2008 = var['2008.dta']
ambiente2009 = var['2009.dta']
ambiente2010 = var['2010.dta']
ambiente2011 = var['2011.dta']
ambiente2012 = var['2012.dta']
ambiente2013 = var['2013.dta']
ambiente2014 = var['2014.dta']
ambiente2015 = var['2015.dta']
ambiente2016 = var['2016.dta']
ambiente2017 = var['2017.dta']

dataframes = [ambiente2008, ambiente2009, ambiente2010, ambiente2011,ambiente2012,ambiente2013,ambiente2014,ambiente2015,ambiente2016,ambiente2017]

for i in dataframes:
    i.drop(i.columns[[0,2,3,4,5,6,10,11]], axis=1, inplace=True)    
dataframes = [ambiente2008, ambiente2009, ambiente2010, ambiente2011,ambiente2012,ambiente2013,ambiente2014,ambiente2015,ambiente2016,ambiente2017]

filename2 = 'mypickle2.pk'

with open(filename2, 'wb') as fi:
    pickle.dump(dataframes, fi)

import os
import pandas as pd
import pickle
import numpy as np
import linearmodels as lm
import numpy.linalg as la
from statsmodels.stats import outliers_influence
import statsmodels.api as sm
from scipy import stats

path=os.getcwd()
os.chdir("")

filename2 = 'mypickle2.pk'
with open(filename2, 'rb') as fi:
    dataframes = pickle.load(fi)


### importando as variáveis independentes
pib = pd.read_csv('pib por municipio.csv',sep=';')
pib['pib'] = pd.to_numeric(pib['pib'])
pib.set_index(['cod','ano'],inplace = True)

pop = pd.read_csv('população.csv',sep=';')
pop['populacao'] = pd.to_numeric(pop['populacao'])
pop.set_index(['cod','ano'],inplace = True)

pibpercapita = pib[['pib']].div(pop[['populacao']].values, axis=0)
logpibpercapita = np.log(pibpercapita)
logpibpercapita.rename(columns={'pib':'logpibpercapita'},inplace=True)

poluicao = pd.read_csv('poluição.csv',sep=';', decimal=',')
poluicao['poluicao'] = pd.to_numeric(poluicao['poluicao'])
poluicao.set_index(['cod','ano'],inplace = True)

gado = pd.read_csv('gado.csv',sep=';')
gado['gado'] = pd.to_numeric(gado['gado'])
gado.set_index(['cod','ano'],inplace = True)

frota = pd.read_csv('carrosf.csv',sep=';')
frota['carros'] = pd.to_numeric(frota['carros'])
frota.set_index(['cod','ano'],inplace = True)

ambiental = pd.concat(dataframes)
ambiental.head()
ambientalmean = ambiental.groupby(['municipio_codigo','ano']).mean()

teste = frota.merge(gado, on=(['cod','ano']), how='left').merge(logpibpercapita, on=(['cod','ano']), how='left').merge(poluicao, on=(['cod','ano']), how='left')#.merge(pop, on=(['cod','ano']), how='left')
teste.rename(columns={'logpib/populacao':'logpibpercapita'},inplace=True) #renomeando pib pra facilitar


# renomear index ambientalmean e dando merge com teste
ambientalmean = ambientalmean.reset_index()
ambientalmean.rename(columns={'municipio_codigo':'cod'}, inplace=True)
ambientalmean.set_index(['cod','ano'],inplace = True)
#ambientalmean.pop('index','level_0')

del ambiental

teste2 = teste.merge(ambientalmean, on=(['cod','ano']), how='left')
pd.set_option('display.float_format', lambda x: '%.2f' % x)



# criando cluster e balanceando [nível estadual]
group = teste2.groupby(level=teste2.index.names)  
teste2 = group.last() 
teste2['poluicao'] = teste2.groupby(level=0)["poluicao"].shift(1) #jogando poluicao 1 pra frente
teste2 = teste2.dropna()
teste2.reset_index(inplace=True)
teste2['coduf'] = teste2['cod']
teste2.set_index(['cod','ano'],inplace = True)
teste2['uf'] = teste2.coduf.astype(str).str[:2].astype(int)
teste2.pop('coduf')

balanceamento = pd.read_csv('balanceamento.csv',sep=';')
balanceamento.set_index(['cod','ano'],inplace = True)
teste2 = balanceamento.merge(teste2, on=(['cod','ano']), how='left')
teste2 = teste2.fillna(0)

vif_data = pd.DataFrame()
vif_data["feature"] = teste2.columns
vif_data["VIF"] = [outliers_influence.variance_inflation_factor(teste2.values, i)  #umidade e as outras ambientais com alto VIF. outra questão é pib
                          for i in range(len(teste2.columns))]

print(vif_data)
teste2.drop(['umidade_relativa_percentual'],axis=1,inplace=True)

#teste2.pop('index')

desc = teste2.describe()
desc.to_excel('descrtivo.xlsx')

# regressões
#aleatorio
random = lm.RandomEffects(teste2[['poluicao']], sm.add_constant(teste2[['gado','carros','logpibpercapita','precipitacao_mmdia','temperatura_c']]))
res = random.fit(cov_type='clustered', clusters=teste2[['uf']], cluster_entity=True)
print(res)

#fixo
fixo = lm.PanelOLS(teste2[['poluicao']], sm.add_constant(teste2[['gado','carros','logpibpercapita','precipitacao_mmdia','temperatura_c']]),entity_effects=True,time_effects = True)

re_res = fixo.fit(cov_type='clustered',clusters=teste2[['uf']],cov_kwds='uf', cluster_entity=True) #cluster?

print(re_res)

# ols 
ols = lm.BetweenOLS(teste2[['poluicao']], sm.add_constant(teste2[['gado','carros','logpibpercapita','precipitacao_mmdia','temperatura_c']]))

resols = ols.fit(cov_type='clustered', cluster_entity=True)
print(resols)

beginningtex = """\\documentclass{report}
\\usepackage{booktabs}
\\begin{document}"""
endtex = "\end{document}"
f = open('myreg.tex', 'w')
f.write(beginningtex)
f.write(resols.summary.as_latex())
f.write(re_res.summary.as_latex())
f.write(res.summary.as_latex())
f.write(endtex)
f.close()

def hausman(re_res, re):
    """
    Compute hausman test for fixed effects/random effects models
    b = beta_fe
    B = beta_re
    From theory we have that b is always consistent, but B is consistent
    under the alternative hypothesis and efficient under the null.
    The test statistic is computed as
    z = (b - B)' [V_b - v_B^{-1}](b - B)
    The statistic is distributed z \sim \chi^2(k), where k is the number
    of regressors in the model.
    Parameters
    ==========
    fe : statsmodels.regression.linear_panel.PanelLMWithinResults
        The results obtained by using sm.PanelLM with the
        method='within' option.
    re : statsmodels.regression.linear_panel.PanelLMRandomResults
        The results obtained by using sm.PanelLM with the
        method='swar' option.
    Returns
    =======
    chi2 : float
        The test statistic
    df : int
        The number of degrees of freedom for the distribution of the
        test statistic
    pval : float
        The p-value associated with the null hypothesis
    Notes
    =====
    The null hypothesis supports the claim that the random effects
    estimator is "better". If we reject this hypothesis it is the same
    as saying we should be using fixed effects because there are
    systematic differences in the coefficients.
    """

    # Pull data out
    b = re_res.params
    B = re.params
    v_b = re_res.cov
    v_B = re.cov

    # NOTE: find df. fe should toss time-invariant variables, but it
    #       doesn't. It does return garbage so we use that to filter
    df = b[np.abs(b) < 1e8].size

    # compute test statistic and associated p-value
    chi2 = np.dot((b - B).T, la.inv(v_b - v_B).dot(b - B))
    pval = stats.chi2.sf(chi2, df)

    return chi2, df, pval

hausman(re_res,res) #rejeitamos H0 -> fe melhor
...