# 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

n, t = 10000, 200
mu, sigma = 0.05, 0.3
dn, cn = n // 10, n // 100

R = np.random.normal(mu, sigma, (n, t))
X = R.cumsum(axis=1)
Y = np.exp(X)

S = np.ones((n, 1))
K = np.concatenate((S, S*Y), axis=1)
K.sort(axis=0)

W_Total = K[0:n].sum(axis=0)
Q_Top10 = K[n-dn:n].sum(axis=0) / W_Total
Q_Top100 = K[n-cn:n].sum(axis=0) / W_Total

T = np.linspace(0, t, num=t+1)
P_Top10 = sts.norm.cdf(sigma*np.sqrt(T)-sts.norm.ppf(1-0.1))
P_Top100 = sts.norm.cdf(sigma*np.sqrt(T)-sts.norm.ppf(1-0.01))

fig, ax = plt.subplots(figsize=(10.0, 5.0))
fig.suptitle('Vermögenskonzentration im Zeitverlauf')
ax.set(xlabel='Zeit (Jahre)', ylabel='Vermögensanteil (100% = 1.0)')

ax.plot(Q_Top10, 'bo', markersize=4.0, label='Die oberen 10% der Population (Simulation)')
ax.plot(T, P_Top10, 'c', linewidth=2.0, label='Die oberen 10% der Population (Prognose)')
ax.plot(Q_Top100, 'bo', markersize=2.0, label='Die oberen 1% der Population (Simulation)')
ax.plot(T, P_Top100, 'c', linewidth=1.0, label='Die oberen 1% der Population (Prognose)')

ax.grid()
ax.legend()
plt.savefig("Vermoegenskonzentration_im_Zeitverlauf.png", dpi=150)