p_estimate_model.py

#!/usr/bin/env python

import os
import copy
import string

import numpy as np
import pandas as pd

import seaborn as sns
import matplotlib.pyplot as plt

import geopandas as gpd
import shapely

import statsmodels
import statsmodels.api as sm

## =============================================

def aggregate_bins(df, bins, first_bin_index=None, last_bin_index=None, factor=1):
	if first_bin_index != None:
		df[bins[first_bin_index]] = df.loc[:, bins[:first_bin_index+1]].sum(axis=1)
	if last_bin_index != None:
		df[bins[last_bin_index]] = df.loc[:, bins[last_bin_index:]].sum(axis=1)
	if first_bin_index != None:
		df = df.drop(columns=bins[:first_bin_index])
	if last_bin_index != None:
		df = df.drop(columns=bins[last_bin_index+1:])
	bins = bins[first_bin_index:last_bin_index+1]
	if factor != 1:
		for c in range(0, len(bins), factor):
			#print(bins[c], bins[c:c+factor])
			df.loc[:, bins[c]] = df.loc[:, bins[c:c+factor]].sum(axis=1)
		columns_drop = [b for b in bins if b not in bins[::factor]]
		df = df.drop(columns=columns_drop)
	return df

## =============================================

DATAPATH = '../data_additional'

## =============================================

ifile = 'climate_panel_all.csv'
df_clim = pd.read_csv(os.path.join(DATAPATH, ifile))

df_clim = df_clim.sort_values(by=['unit', 'year'], ascending=True)\
			.rename(columns={'unit': 'GID_1'})

## =============================================

ifile = 'DOSE/DOSE_V2.csv'
df_econ = pd.read_csv(os.path.join(DATAPATH, ifile))

## =============================================

df = df_clim.merge(df_econ, on=['GID_1', 'year'], how='left')
df = df.rename(columns={'GID_1': 'unit'})

df = df.loc[df['unit'].isin(df.loc[df['grp_lcu'].notnull(), 'unit'].unique()), :]
df = df.loc[df['year'].between(1970, 2020), :]

## =============================================

temperature_variables = [c for c in df.columns if (c[:3] == 'bin')]

## ==========================================================================================

## default specification

frequency = 0.005

dfg = df.loc[df['grp_pc_usd'].notnull(), temperature_variables].sum(axis=0)
dfg = dfg.cumsum(axis=0)
dfg = dfg / dfg.iloc[-1]

bins = dfg.loc[(dfg >= frequency) & (dfg <= (1. - frequency))].index.values

## =============================================

first_bin_index = int(bins[0].split('bin')[-1])
last_bin_index = int(bins[-1].split('bin')[-1])

if first_bin_index % 2 == 0:
	first_bin_index = first_bin_index + 1
else:
	pass # Odd

dfx = aggregate_bins(df, bins=temperature_variables,
			first_bin_index=first_bin_index, last_bin_index=last_bin_index,
			factor=2)

#v = [b for b in dfx.columns if (b[:3] == 'bin')]
#dfx.loc[:, v].sum(axis=1)

explanatory_variables_all = [b for b in dfx.columns if (b[:3] == 'bin')] + ['tp_total']
explanatory_variables = [b for b in explanatory_variables_all if (b != 'bin61')]
dfx['year_str'] = dfx['year'].astype(str)
dfx['log_grp_pc_usd'] = np.log(dfx['grp_pc_usd']).replace([-np.inf, np.inf], np.nan)
dfx = dfx.sort_values(by=['unit', 'year'], ascending=True)
dfx['y'] = dfx.groupby('unit')['log_grp_pc_usd'].diff()

formula = 'y ~ year_str + unit + {0:s}'.format(' + '.join(explanatory_variables))
res = sm.OLS.from_formula(formula, dfx).fit(missing='drop').get_robustcov_results()
dx = res.summary2(float_format="%.5f").tables[1].iloc[:, [0, 1, 3]]

dfxr = dfx.dropna(subset=explanatory_variables + ['y'])
res = sm.OLS.from_formula(formula, dfxr).fit(cov_type='cluster', cov_kwds={'groups': dfxr['unit']})
dxc1 = res.summary2(float_format="%.5f").tables[1].iloc[:, [1]].rename(columns={'Std.Err.': 'Std.Err.Cl.Unit'})

dfxr = dfx.dropna(subset=explanatory_variables + ['y'])
res = sm.OLS.from_formula(formula, dfxr).fit(cov_type='cluster', cov_kwds={'groups': dfxr['country']})
dxc2 = res.summary2(float_format="%.5f").tables[1].iloc[:, [1]].rename(columns={'Std.Err.': 'Std.Err.Cl.Country'})

dx = dx.merge(dxc1, left_index=True, right_index=True, how='left')
dx = dx.merge(dxc2, left_index=True, right_index=True, how='left')

dx.loc[explanatory_variables].to_csv('../damage_frompython/coefficients_model_global_01.csv')

## =============================================

def extract_errorbar(dx, variable):
	if variable in dx.index.values:
		y = dx.loc[variable, 'Coef.']
		y_error = 1.96*dx.loc[variable, 'Std.Err.Cl.Country']
	else:
		y = 0.
		y_error = 0.
	return y, y_error

df_limits = pd.read_csv(os.path.join(DATAPATH, 'bins.csv'), sep=',')
df_limits = df_limits.set_index('Name')

variables = [c for c in explanatory_variables_all if 'tp' not in c]
fig, axes = plt.subplots(nrows=2, height_ratios=[4,1], sharex=True)
ax = axes[0]
for x, v in enumerate(variables):
	y, yerror = extract_errorbar(dx, v)
	ax.errorbar(x, y, yerror, capsize=4., elinewidth=2., marker='o', color='b')
ax.plot([-1.5, len(variables) + 0.5], [0., 0.], 'k-', lw=0.5)
xticks = np.arange(0., len(variables) - 1., 1.) + 0.5
ax.set_xticks(xticks)
xticklabels = ['{0:.0f}'.format(x) for x in df_limits.loc[variables[1:], 'Lower bound'].astype(float)]
ax.set_xticklabels(xticklabels, fontsize='x-small', rotation=90.)
ax.set_xlim([-1.5, len(variables) + 0.5])
ax.set_ylabel('Daily mean temperature')
ax.set_ylabel('Change in growth rate\n from one additional day in bin \nrelative to a day in bin (20, 22]')
#ax.annotate(text='{0:d}-{1:d}'.format(year1, year2), xy=(0.05, 0.9), xycoords='axes fraction', ha='left', va='top')
#sns.despine(ax=ax, offset=1., right=True, top=True)
ax = axes[1]
y = dfx.loc[dfx['y'].notnull(), variables].sum(axis=0)
ax.bar([xticks[0]-0.5] + list(xticks+0.5), y, width=1., color='grey')
ax.set_ylabel('N')
#ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels, fontsize='x-small', rotation=90.)
fig.savefig('../damage_frompython/coefficients_model_global_01.pdf', bbox_inches='tight', transparent=True)

## ==========================================================================================

## first differences of climate variable

for variable in explanatory_variables:
	dfx[variable] = dfx.groupby('unit')[variable].diff()


formula = 'y ~ year_str + unit + {0:s}'.format(' + '.join(explanatory_variables))
res = sm.OLS.from_formula(formula, dfx).fit(missing='drop').get_robustcov_results()
dx = res.summary2(float_format="%.5f").tables[1].iloc[:, [0, 1, 3]]

dfxr = dfx.dropna(subset=explanatory_variables + ['y'])
res = sm.OLS.from_formula(formula, dfxr).fit(cov_type='cluster', cov_kwds={'groups': dfxr['unit']})
dxc1 = res.summary2(float_format="%.5f").tables[1].iloc[:, [1]].rename(columns={'Std.Err.': 'Std.Err.Cl.Unit'})

dfxr = dfx.dropna(subset=explanatory_variables + ['y'])
res = sm.OLS.from_formula(formula, dfxr).fit(cov_type='cluster', cov_kwds={'groups': dfxr['country']})
dxc2 = res.summary2(float_format="%.5f").tables[1].iloc[:, [1]].rename(columns={'Std.Err.': 'Std.Err.Cl.Country'})

dx = dx.merge(dxc1, left_index=True, right_index=True, how='left')
dx = dx.merge(dxc2, left_index=True, right_index=True, how='left')

dx.loc[explanatory_variables].to_csv('../damage_frompython/coefficients_model_global_01fd.csv')

## =============================================

def extract_errorbar(dx, variable):
	if variable in dx.index.values:
		y = dx.loc[variable, 'Coef.']
		y_error = 1.96*dx.loc[variable, 'Std.Err.Cl.Country']
	else:
		y = 0.
		y_error = 0.
	return y, y_error

df_limits = pd.read_csv(os.path.join(DATAPATH, 'bins.csv'), sep=',')
df_limits = df_limits.set_index('Name')

variables = [c for c in explanatory_variables_all if 'tp' not in c]
fig, axes = plt.subplots(nrows=2, height_ratios=[4,1], sharex=True)
ax = axes[0]
for x, v in enumerate(variables):
	y, yerror = extract_errorbar(dx, v)
	ax.errorbar(x, y, yerror, capsize=4., elinewidth=2., marker='o', color='b')
ax.plot([-1.5, len(variables) + 0.5], [0., 0.], 'k-', lw=0.5)
xticks = np.arange(0., len(variables) - 1., 1.) + 0.5
ax.set_xticks(xticks)
xticklabels = ['{0:.0f}'.format(x) for x in df_limits.loc[variables[1:], 'Lower bound'].astype(float)]
ax.set_xticklabels(xticklabels, fontsize='x-small', rotation=90.)
ax.set_xlim([-1.5, len(variables) + 0.5])
ax.set_ylabel('Daily mean temperature')
ax.set_ylabel('Change in growth rate\n from one additional day in bin \nrelative to a day in bin (20, 22]')
#ax.annotate(text='{0:d}-{1:d}'.format(year1, year2), xy=(0.05, 0.9), xycoords='axes fraction', ha='left', va='top')
#sns.despine(ax=ax, offset=1., right=True, top=True)
ax = axes[1]
y = dfx.loc[dfx['y'].notnull(), variables].sum(axis=0)
ax.bar([xticks[0]-0.5] + list(xticks+0.5), y, width=1., color='grey')
ax.set_ylabel('N')
#ax.set_xticks(xticks)
ax.set_xticklabels(xticklabels, fontsize='x-small', rotation=90.)
fig.savefig('../damages_frompython/coefficients_model_global_01fd.pdf', bbox_inches='tight', transparent=True)