# Entrepreneurs, Chance, and the Deterministic Concentration of Wealth
# Joseph E. Fargione, Clarence Lehman, Stephen Polasky
# https://doi.org/10.1371/journal.pone.0020728

import numpy as np
import scipy.stats as sts
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick

n, t = 100000, 200  # Größe der Population und Zeitraum in Jahren
p = -0.2, 0.3  # Auflistung der möglichen jährlichen Wachstumsraten (Auswahl per Münzwurf)

P = np.random.choice(p, (n, t))  # Matrix mit den zufälligen individuellen Wachstumsraten p_i_k
R = np.log(1 + P)  # Matrix mit den entsprechenden zufälligen individuellen Raten r_i_k = ln(1 + p_i_k)
S = np.ones((n, 1))  # Vektor für das einheitliche Startkapital S_i = K_0 = 1 (z.B. 1 Mio. Euro)
T = np.linspace(0, t, num=t+1)  # Liste für die Jahre 0, 1, 2, ..., t zur Berechnung der Prognosen

fig, ax = plt.subplots(figsize=(12.0, 6.0))
fig.suptitle('Entwicklung der Vermögensanteile der reichsten zehn Prozent der Population mit und ohne Zinseszinsen')
ax.set(xlabel='Zeit (in Jahren)', ylabel='Anteil am Gesamtvermögen (in Prozent)')

X = R.cumsum(axis=1)  # Matrix mit den Summen x_i(t) = r_i_1 + r_i_2 + ... + r_i_t
Y = np.exp(X)  # Matrix mit den Potenzen y_i(t) = e^(x_i(t)) = e^(r_i_1 + r_i_2 + ... + r_i_t)
K = np.concatenate((S, S * Y), axis=1)  # Matrix mit den Kapitalvermögen K_i(t) = S_i * y_i(t) = K_0 * e^(x_i(t))
K.sort(axis=0)  # Sortierung der Kapitalvermögen innerhalb der einzelnen Jahre 0, 1, 2, ..., t

W_Total = K[0:n].sum(axis=0)  # Liste mit den Gesamtvermögen für die Jahre 0, 1, 2, ..., t
W_Top = K[(n - n // 10):n].sum(axis=0)  # Liste mit den Teilvermögen der reichsten zehn Prozent für diese Jahre
Q_Top = W_Top / W_Total  # Liste mit den Vermögensanteilen der reichsten zehn Prozent für diese Jahre

sigma = np.std(np.log(1 + np.array(p)))  # Standardabweichung der Raten r = ln(1 + p)
P_Top = sts.norm.cdf(sigma * np.sqrt(T) - sts.norm.ppf(1 - 0.1))  # Liste mit den prognostizierten Vermögensanteilen

ax.plot(Q_Top, 'bo', markersize=4.0, label='Vermögen im Besitz der reichsten zehn Prozent der Population (Simulation mit Zinseszinsen)')
ax.plot(T, P_Top, 'c', linewidth=2.0, label='Vermögen im Besitz der reichsten zehn Prozent der Population (Prognose mit Zinseszinsen)')

X = P.cumsum(axis=1)  # Matrix mit den Summen x_i(t) = p_i_1 + p_i_2 + ... + p_i_t
Y = 1 + X  # Matrix mit den Summen y_i(t) = 1 + x_i(t) = 1 + p_i_1 + p_i_2 + ... + p_i_t
K = np.concatenate((S, S * Y), axis=1)  # Matrix mit den Kapitalvermögen K_i(t) = S_i * y_i(t) = K_0 * (1 + x_i(t))
K.sort(axis=0)  # Sortierung der Kapitalvermögen innerhalb der einzelnen Jahre 0, 1, 2, ..., t

W_Total = K[0:n].sum(axis=0)  # Liste mit den Gesamtvermögen für die Jahre 0, 1, 2, ..., t
W_Top = K[(n - n // 10):n].sum(axis=0)  # Liste mit den Teilvermögen der reichsten zehn Prozent für diese Jahre
Q_Top = W_Top / W_Total  # Liste mit den Vermögensanteilen der reichsten zehn Prozent für diese Jahre

mu, sigma = np.mean(p), np.std(p)  # Mittelwert und Standardabweichung der Wachstumsraten p
P_Top = 0.1 + sigma * np.exp(- sts.norm.ppf(1 - 0.1) ** 2 / 2) * np.sqrt(T / np.pi / 2) / (mu * T + 1)  # Liste mit den Prognosen

ax.plot(Q_Top, 'bo', markersize=2.0, label='Vermögen im Besitz der reichsten zehn Prozent der Population (Simulation ohne Zinseszinsen)')
ax.plot(T, P_Top, 'c', linewidth=1.0, label='Vermögen im Besitz der reichsten zehn Prozent der Population (Prognose ohne Zinseszinsen)')

ax.legend(loc='lower right', fontsize='small')
ax.set_xticks(np.linspace(0, t, num=21))
ax.set_yticks(np.linspace(0, 1, num=11))
ax.set_ylim([-0.05, 1.05])
ax.grid()
plt.gca().yaxis.set_major_formatter(mtick.PercentFormatter(xmax=1.0))
plt.savefig("Vermoegenskonzentration_Simulation_2.svg")