# 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
mu, sigma = 0.05, 0.3  # Mittelwert und Standardabweichung (Normalverteilung)

R = np.random.normal(mu, sigma, (n, t))  # Matrix mit den zufälligen individuellen Raten r_i_k
S = np.ones((n, 1))  # Vektor mit gleichverteiltem Startkapital s_i = 1 für i = 1, 2, ..., n
T = np.linspace(0, t, num=t+1)  # Liste für die Jahre k = 0, 1, 2, ..., t zur Berechnung der Prognosen

X = R.cumsum(axis=1)  # Matrix mit den Summen x_i_k = r_i_1 + r_i_2 + ... + r_i_k für k = 1, 2, ..., t
Y = np.exp(X)  # Matrix mit den Potenzen y_i_k = e^(x_i_k) = e^(r_i_1 + r_i_2 + ... + r_i_k) für k = 1, 2, ..., t
W = np.concatenate((S, S * Y), axis=1)  # Matrix mit den Kapitalvermögen w_i_k = s_i * y_i_k für k = 0, 1, 2, ..., t
W.sort(axis=0)  # Sortierung der Kapitalvermögen (wealth) innerhalb der einzelnen Jahre k = 0, 1, 2, ..., t

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

W_Total = W[0:n].sum(axis=0)  # Liste mit den Gesamtvermögen für die Jahre k = 0, 1, 2, ..., t
W_Top = W[(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
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)')
ax.plot(T, P_Top, 'c', linewidth=2.0, label='Vermögen im Besitz der reichsten zehn Prozent der Population (Prognose)')

W_Total = W[0:n].sum(axis=0)  # Liste mit den Gesamtvermögen für die Jahre k = 0, 1, 2, ..., t
W_Top = W[(n - n // 100):n].sum(axis=0)  # Liste mit den Teilvermögen des reichsten Prozents für diese Jahre
Q_Top = W_Top / W_Total  # Liste mit den Vermögensanteilen des reichsten Prozents für diese Jahre
P_Top = sts.norm.cdf(sigma * np.sqrt(T) - sts.norm.ppf(1 - 0.01))  # Liste mit den prognostizierten Vermögensanteilen

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

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_1.svg")