import numpy as np
import pandas as pd
import plotly.figure_factory as ff
from scipy.stats import entropy, mielke
pd.options.plotting.backend = "plotly"
Intuition für Verteilungsmaße#
Um ein Gefühl für die Verteilungsmaße zu bekommen, illustrieren wir Sie auf Basis von geschätzten Verteilungen. Wir nehmen geschätzte Parameter für eine Dagum-Verteilung (in Scipy implementiert als Mielke-Verteilung nach vorheriger Wiederentdeckung der Burr-Typ-III-Verteilung) für Deutschland, Frankreich, Niederlande und Polen im Jahr 2011, um die Daten zu simulieren.
Die Quelle der geschätzten Parameter ist: Koen Decancq und Philippe Van Kerm (Hrsg.), What Drives Inequality?, Emerald Publishing Group 2019.
NB:
Die Daten sind auf Basis des EU-SILC geschätzt und die parametrische Verteilung deckt den rechten Rand nicht gut ab.
Alle Daten sind für das kaufkraftbereinigte Äquivalenzeinkommen; das genaue Verfahren für die Kaufkraftbereinigung konnten wir nicht nachvollziehen, in jedem Fall schauen Sie nicht zu genau auf die Vergleiche über Länder hinweg. Das heißt, behalten Sie bitte keine Aussagen wie “Die Niederlande haben eine geringere Einkommensungleichheit als Deutschland” im Hinterkopf.
Ziehung von je 100.000 Beobachtungen#
Zunächst setzen wir den Zustand des Zufallszahlengenerators, damit wir auf allen Rechnern die gleichen Werte bekommen. Bei echten Anwendungen sollte die Zahl der Ziehungen immer so groß sein, dass die Ergebnisse nicht von dieser Zahl abhängen!
np.random.seed(seed=233423) # noqa: NPY002
n = 100_000
daten = pd.DataFrame(index=pd.RangeIndex(n))
Das Sortieren macht unten die Berechnung der Lorenzkurven einfacher.
scale = 20_394
s = 3.8206
k = s * 0.756
daten["DE"] = sorted(mielke.rvs(k=k, s=s, scale=scale, size=n))
scale = 19_274
s = 4.0893
k = s * 0.9489
daten["NL"] = sorted(mielke.rvs(k=k, s=s, scale=scale, size=n))
scale = 17_390
s = 3.3272
k = s * 1.1053
daten["FR"] = sorted(mielke.rvs(k=k, s=s, scale=scale, size=n))
scale = 9_073
s = 3.4433
k = s * 0.7843
daten["PL"] = sorted(mielke.rvs(k=k, s=s, scale=scale, size=n))
Deskriptive Statistiken#
daten.describe().round(-2)
Plots der Verteilungen#
Um Einkommensverteilungen zu schätzen, nutzt man in der Regel eine nichtparametrische Kerndichteschätzung (kernel density estimation, KDE).
Am einfachsten können Sich sich diese als eine stetige Verallgemeinerung eines Histograms vorstellen, wie sie an der nächsten Graphik sehen. Die Zahl der Histogrammintervalle ist absichtlich so gewählt, dass es offensichtlich wird. Wir nutzen hier aufgrund des schweren rechten Rands der Verteilung, den wir zunächst abschneiden, sehr viele Intervalle.
max_income = 75_000
fig = ff.create_distplot(
[daten["DE"][daten["DE"] < max_income]],
["DE"],
bin_size=5000,
show_rug=False,
)
fig.update_layout(
title="Geschätzte Verteilung der Einkommen in Deutschland",
xaxis_title="Einkommen in Euro",
yaxis_title="Dichte",
)
fig.show()
max_income = 75_000
fig = ff.create_distplot(
[daten["DE"][daten["DE"] < max_income]],
["DE"],
bin_size=5000,
show_rug=False,
)
fig.update_layout(
title="Geschätzte Verteilung der Einkommen in Deutschland",
xaxis_title="Einkommen in Euro",
yaxis_title="Dichte",
)
fig.show()
fig = ff.create_distplot(
[daten[c] for c in daten.columns],
daten.columns,
bin_size=5000,
show_rug=False,
)
fig.update_layout(
title="Geschätzte Verteilung der Einkommen",
xaxis_title="Einkommen in Euro",
yaxis_title="Dichte",
)
fig.show()
max_income = 75_000
fig = ff.create_distplot(
[daten[c].loc[daten[c] < max_income] for c in daten.columns],
daten.columns,
bin_size=5000,
show_rug=False,
)
fig.update_layout(
title="Geschätzte Verteilung der Einkommen",
xaxis_title="Einkommen in Euro",
yaxis_title="Dichte",
)
fig.show()
Perzentile#
def perzentile(daten):
"""Berechnet die Perzentile der Daten.
Args:
daten (pd.DataFrame): Daten
Returns:
pd.Series: Perzentile
"""
p = np.arange(1, 100)
ausgabe = daten.quantile(p / 100)
ausgabe.index = p
return ausgabe
perz = perzentile(daten)
perz.round(-2)
perz.plot.scatter(x=perz.index, y=perz.columns, title="Perzentile der Einkommen")
def mittelwert_pro_perzentil(daten):
"""Berechnet den Mittelwert pro Perzentil.
Args:
daten (pd.DataFrame): Daten
Returns:
pd.DataFrame: Mittelwert pro Perzentil
"""
ausgabe = pd.DataFrame(index=range(1, 101))
_tmp = daten.copy()
for c in daten.columns:
_tmp[c + "_perzentil"] = pd.qcut(_tmp[c], q=100, labels=ausgabe.index)
ausgabe[c] = _tmp.groupby(c + "_perzentil")[c].mean()
return ausgabe
mw_pro_perzentil = mittelwert_pro_perzentil(daten)
mw_pro_perzentil.round(-2)
mw_pro_perzentil.plot.scatter(
x=mw_pro_perzentil.index,
y=mw_pro_perzentil.columns,
title="Mittleres Einkommen pro Perzentil",
)
Palma-Ratio#
(mw_pro_perzentil.loc[91:].sum() / mw_pro_perzentil.loc[:40].sum()).to_frame().rename(
columns={0: "Palma-Ratio"},
).round(2)
Armutsrisikoquote#
def armutsrisikoquote(einkommen):
"""Berechnet die Armutsrisikoquote.
Args:
einkommen (pd.Series): Einkommen
Returns:
float: Armutsrisikoquote
"""
n = einkommen.size
schwelle = einkommen.quantile(0.5) * 0.6
n_unter_schwelle = einkommen[einkommen < schwelle].size
return n_unter_schwelle / n
daten.apply(armutsrisikoquote).to_frame().rename(
columns={0: "Armutsrisikoquote"},
).round(3)
Theil-Index#
Der Theil-Index ist allgemeines Entropiemaß mit $\alpha = 1$.
daten.apply(entropy).to_frame().rename(columns={0: "Theil"}).round(3).sort_values(
"Theil",
)
Lorenzkurve#
def lorenzkurve(daten):
"""Berechnet die Lorenzkurve.
Args:
daten (pd.DataFrame): Daten
Returns:
pd.DataFrame: Lorenzkurve
"""
lorenz = daten.copy()
for c in lorenz.columns:
lorenz[c] = lorenz[c].cumsum() / lorenz[c].sum()
lorenz["Gleichverteilung"] = lorenz.index / lorenz.index.max()
lorenz.index /= lorenz.index.max()
return lorenz
lorenz = lorenzkurve(daten)
lorenz
lorenz.plot.line(x=lorenz.index, y=lorenz.columns, title="Lorenzkurve")
Gini Koeffizient#
def gini(einkommen_sortiert):
"""Berechnet den Gini-Koeffizienten.
Args:
einkommen_sortiert (pd.Series): Einkommen sortiert
Returns:
float: Gini-Koeffizient
"""
n = einkommen_sortiert.size
coef = 2 / n
const = (n + 1) / n
weighted_sum = sum([(i + 1) * yi for i, yi in enumerate(einkommen_sortiert)])
return coef * weighted_sum / (einkommen_sortiert.sum()) - const
daten.apply(gini).to_frame().rename(columns={0: "Gini"}).round(3)
Aufgabe: Stellen wir uns vor, die Londoner City migriert komplett nach Amsterdam…#
Nun geht es darum, zu sehen, inwieweit die Maße auf Änderungen am rechten Rand der Verteilung reagieren. Wir setzen daher in den Niederlanden das oberste Perzentil der Einkommen auf 250_000 €.
daten_verändert = daten[["DE", "NL"]].copy()
daten_verändert["NL + Banker"] = daten["NL"]
_idx = daten_verändert.iloc[-1_000:].index
daten_verändert.loc[_idx, "NL + Banker"] = 250_000
daten_verändert.describe().round(-2)
Plot der Dichtefunktionen#
fig = ff.create_distplot(
[daten_verändert[c] for c in daten_verändert.columns],
daten_verändert.columns,
bin_size=1000,
show_rug=False,
)
fig.update_layout(
title="Geschätzte Verteilung der Einkommen",
xaxis_title="Einkommen in Euro",
yaxis_title="Dichte",
)
fig.show()
Plot der Perzentile#
Plot der Mittelwerte pro Perzentil#
Lorenzkurven mit Beispiel, dass diese sich schneiden können#
lorenz_verändert = lorenzkurve(daten_verändert)
lorenz_verändert
lorenz_verändert.plot.line(
x=lorenz_verändert.index,
y=lorenz_verändert.columns,
title="Lorenzkurve",
)
fig = lorenz_verändert.plot.line(
x=lorenz_verändert.index,
y=lorenz_verändert.columns,
title="Lorenzkurve",
)
fig.update_xaxes(range=[0.2, 0.7])
fig.update_yaxes(range=[0.05, 0.55])
fig.show()
Gini-Koeffizient#
Palma-Ratio#
Armutsrisikoquote#
Wie ranken Sie die Einkommensverteilungen in daten_verändert
hinsichtlich der Ungleichheit? Begründen Sie Ihre Antwort!#
Hier Platz für Ihre Antwort, bitte laden Sie diese auch in eCampus hoch