Vortrag für die Generalversammlung der AVÖ
17. Mai 2018
# "Standard" Bibliotheken
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import uniform, norm
from pandas_datareader import data
from scipy.fftpack import rfft, irfft, fft, rfft, dct, idct
from scipy.cluster.vq import kmeans2 as kmeans
from mpl_toolkits.mplot3d import Axes3D
from sklearn.preprocessing import scale as scale_this
%matplotlib inline
# !!! DIE FOLGENDEN PFADNAMEN MÃœSSTEN AN LOKALE GEGEBENHEITEN ANGEPASST WERDEN !!!
# Daten vom statistischen Bundesamt Deutschland:
DATAPATH = '/Users/mfulmek/Work/Brotberuf/Vortraege/AVOE2018/destatis/'
# Graphiken:
GRAPHPATH = '/Users/mfulmek/Work/Brotberuf/Vortraege/AVOE2018/graphics/'
Das deutsche Bundesamt für Statistik stellt unter dem Link https://www-genesis.destatis.de/genesis/online/logon?language=de&sequenz=tabelleErgebnis&selectionname=12613-0003&zeitscheiben=5 eine Tabelle der Todesfälle in Deutschland für die Jahre 2011-2015 zur Verfügung, gegliedert nach Geschlecht und Lebensjahr: Ich habe das entsprechende Datenfile (12613-0003.xlsx) heruntergeladen und kann es nun einlesen.
GD = pd.read_excel(DATAPATH+'12613-0003.xlsx')
# Die ersten 6 Datenzeilen
GD[:6]
# Die Datenzeilen 105 bis 109
GD[105:110]
# Die Datenzeilen 208 und 209
GD[208:210]
# Trennen der Kollektive: Todesfälle Männer/Frauen in Deutschland
GMD = GD[5:106] # Zeilen entsprechend Todesfälle Männer
GFD = GD[109:210] # Zeilen entsprechend Todesfälle Frauen
# Überflüssige Informationen entfernen, Spaltenüberschriften = Jahre 2011 bis 2015
COLUMNS = [str(year) for year in range(2011,2016)]
GMD = pd.DataFrame(GMD.values[:,1:], columns=COLUMNS)
GFD = pd.DataFrame(GFD.values[:,1:], columns=COLUMNS)
# "Interne Buchhaltung" für unsere Daten:
DATEN = {
'm' : {
'bez' : 'Männer',
'roh' : GMD
},
'f' : {
'bez' : 'Frauen',
'roh' : GFD
}
}
# Hilfsfunktion zum Plotten der Rohdaten:
def rohdaten_plot(gender):
"""Hilfsfunktion: Plot der Todesfall-Rohdaten für Männer ('m') oder Frauen ('f')"""
myplot = DATEN[gender]['roh'].plot(
title=f'Todesfälle: {DATEN[gender]["bez"]} in Deutschland, 2011-2015', figsize=(12,8)
)
myplot.set_xlabel('Alter')
myplot.set_ylabel('Anzahl Todesfälle')
myplot.grid(b=True)
plt.xticks(np.arange(0, 101, 5))
plt.savefig(GRAPHPATH+f'roh_{gender}.eps')
return myplot
rohdaten_plot('m')
rohdaten_plot('f')
Rein "optisch": Bis zum Pensionsalter wenig Unterschiede, danach anscheinend "Verschiebung" der "lokalen Maxima" nach "oben rechts".
Sehr salopp gesagt: "Differenzieren verstärkt Schwankungen, Integrieren glättet Schwankungen aus."
# Wir bilden die ersten Differenzen
GMD_DIFF = pd.DataFrame(GMD[1:].values - GMD[:-1].values, columns=COLUMNS)
GFD_DIFF = pd.DataFrame(GFD[1:].values - GFD[:-1].values, columns=COLUMNS)
# Wir berücksichtigen diese neuen Daten in unserer "Daten-Buchhaltung"
DATEN['m']['delta'] = GMD_DIFF
DATEN['f']['delta'] = GFD_DIFF
# Eine Hilfsfunktion zum Plotten der Differenzen:
def diff_rohdaten_plot(gender):
"""Hilfsfunktion: Plot der ersten Differenzen der Todesfall-Rohdaten für Männer ('m') oder Frauen ('f')"""
myplot = DATEN[gender]['delta'].plot(
title=f'1. Differenzen Todesfälle: {DATEN[gender]["bez"]} in Deutschland, 2011-2015', figsize=(12,8)
)
myplot.set_xlabel('Alter')
myplot.set_ylabel('Delta Todesfälle')
plt.xticks(np.arange(0, 101, 5))
myplot.grid(b=True)
plt.savefig(GRAPHPATH+f'delta_roh_{gender}.eps')
return myplot
diff_rohdaten_plot('m')
diff_rohdaten_plot('f')
wieder "rein optisch": Daten zeigen ab Alter etwa 65 deutlich anderes "qualitatives Verhalten".
# Eine Hilfsfunktion zum Plotten der Differenzen für einzelne Jahre:
def diff_rohdaten_plot_jahr(gender, years):
"""Hilfsfunktion: Plot der ersten Differenzen der Todesfall-Rohdaten für Männer ('m') oder Frauen ('f')"""
fig, ax = plt.subplots(figsize=(12,8))
for year in years:
ax.plot(DATEN[gender]['delta'][year], label=f'{DATEN[gender]["bez"]} {year}')
ax.legend(loc='lower left', shadow=True, fontsize='x-large')
ax.set_title(f'1. Differenzen Todesfälle: {DATEN[gender]["bez"]} in Deutschland')
ax.set_xlabel('Alter')
ax.set_ylabel('Delta Todesfälle')
plt.xticks(np.arange(0, 101, 5))
ax.grid(b=True)
plt.savefig(GRAPHPATH+f'delta_roh_{gender}_{years[0]}-{years[-1]}.eps')
plt.show()
diff_rohdaten_plot_jahr('m', ['2011','2012'])
diff_rohdaten_plot_jahr('m', ['2012','2013'])
diff_rohdaten_plot_jahr('m', ['2011','2015'])
Entwicklung eines Vektors nach einem vollständigen Orthonormalsystem in einem Hilbertraum.
# Eine Hilfsfunktion zum Plotten der Fourier-Anpassung:
def diff_fourier_plot(gender, year, start, k):
"""Hilfsfunktion: Plot der ersten Differenzen der Todesfall-Rohdaten für Männer ('m') oder Frauen ('f')"""
data = DATEN[gender]['delta'][year][start:]
# Gläten der Daten mit Fourier-Transformation in 3 Zeilen:
fcoeffs = rfft(data) # Fourier-Transformation
fcoeffs[k+1:] = 0.0 # Entfernen der hohen Frequenzen
fitted_data = irfft(fcoeffs) # Inverse Fourier-Transformation
# Formatierte Graphik
fig, ax = plt.subplots(figsize=(12,8))
ax.plot(data.values, label='Rohdaten')
ax.set_title(f'1. Differenzen Todesfälle {DATEN[gender]["bez"]} in {year}, ab Alter {start}')
ax.plot(fitted_data, label=f'Fourier-Filter (k = {k})')
ax.set_xlabel('Alter')
ax.set_ylabel('Delta Todesfälle')
legend = ax.legend(loc='lower left', shadow=True, fontsize='x-large')
xticks_range = np.arange(0, 100 - start + 1, 5)
labels = plt.xticks(xticks_range, xticks_range + start)
ax.grid(b=True)
plt.savefig(GRAPHPATH+f'ff_{gender}_{year}_{start}_{k}.eps')
plt.show()
return fcoeffs[:k]
diff_fourier_plot('m', '2011', 50, 22)
diff_fourier_plot('m', '2012', 65, 15)
diff_fourier_plot('m', '2011', 20, 36)
AGGREGIERT = pd.DataFrame(
np.array([list(GMD.mean(axis=1).values), list(GFD.mean(axis=1).values)]).T,
columns=['Männer', 'Frauen']
)
AGGREGIERT.head()
myplot = AGGREGIERT.plot(title='Todesfälle in Deutschland 2011-2015, aggregiert', figsize=(12,8))
myplot.set_xlabel('Alter')
myplot.set_ylabel('Anzahl Todesfälle')
plt.savefig(DATAPATH+'aggregiert.eps')
myplot
Es fehlen hier die Informationen, wieviele x-Jährige durch die Todesfälle im x-ten Lebensjahr dezimiert wurden, es ist also nicht möglich, die Sterbewahrscheinlichkeiten $$q_x\sim\frac{\text{im x-ten Jahr Verstorbene}}{\text{x-Jährige}}$$ ohne weiteres aus den Daten zu schätzen.
Um dennoch weiterzumachen (es geht hier ja um die Illustration von Software-Tools), fasse ich die vorliegenden Zahlen als erste Differenzen einer Ausscheideordnung auf (mit "Startwert" gleich der Summe über alle Todesfälle):
(Bei der sachlich richtigen Vorgangsweise mit Kohortensterbetafeln müßten die Wanderungssalden berücksichtig werden, siehe Kohortensterbetafeln für Deutschland/destatis).
AOM = np.r_[GMD.sum(axis=1).values[::-1].cumsum()[::-1], np.zeros(1)]
AOM = 10**6 / AOM[0] * AOM
AOF = np.r_[GFD.sum(axis=1).values[::-1].cumsum()[::-1], np.zeros(1)]
AOF = 10**6 / AOF[0] * AOF
AO_AGGREGIERT = pd.DataFrame(
np.array([AOM, AOF]).T,
columns=['Männer', 'Frauen']
)
myplot = AO_AGGREGIERT.plot(title='Nur für Illustrationszwecke!\nInkorrekt geschätzte Absterbeordnung Deutschland, aggregiert.', figsize=(12,8))
myplot.set_xlabel('Alter')
myplot.set_ylabel('Ãœberlebende')
plt.xticks(np.arange(0, 101, 5))
myplot.grid(b=True)
plt.savefig(GRAPHPATH+'ao_aggregiert.eps')
myplot
GMD_SUMS = dict(zip(GMD.columns, [GMD[column].sum() for column in GMD.columns]))
GFD_SUMS = dict(zip(GFD.columns, [GFD[column].sum() for column in GFD.columns]))
GMD_AO=pd.DataFrame(
np.array([np.r_[GMD[column].values[::-1].cumsum()[::-1]/GMD_SUMS[column]*(10**6), np.zeros(1)] for column in GMD.columns]).T,
columns = GMD.columns
)
GFD_AO=pd.DataFrame(
np.array([np.r_[GFD[column].values[::-1].cumsum()[::-1]/GFD_SUMS[column]*(10**6), np.zeros(1)] for column in GFD.columns]).T,
columns = GFD.columns
)
DATEN['m']['ao'] = GMD_AO
DATEN['f']['ao'] = GFD_AO
# Eine Hilfsfunktion zum Plotten der - sehr ungenau geschätzten! - Absterbeordnung:
def ao_plot(gender):
"""Hilfsfunktion: Plot der sehr ungenau geschätzten Absterbeordnung für Männer ('m') oder Frauen ('f')"""
myplot = DATEN[gender]['ao'].plot(
title=f'Nur für Illustrationszwecke!\nInkorrekt geschätzte Absterbeordnung: {DATEN[gender]["bez"]} in Deutschland, 2011-2015.', figsize=(12,8)
)
myplot.set_xlabel('Alter')
myplot.set_ylabel('Ãœberlebende')
plt.xticks(np.arange(0, 101, 5))
myplot.grid(b=True)
plt.savefig(GRAPHPATH+f'ao_{gender}.eps')
return myplot
ao_plot('m')
ao_plot('f')
Welche Unterschiede ergeben sich in den versicherungsmathematischen Barwerten, wenn man die 5 verschiedenen (nur für Illustrationszwecke geschätzten!!!) Sterbetafeln 2011-2015 betrachtet?
def actuarial_pv(x, ao, rz):
"""Barwert einer lebenslänglichen nachschüssigen Leibrente ab Alter x"""
survival_probabilities = ao[x+1:]/ao[x]
discount_factors = np.array([(1+rz/100) ** (-t) for t in range(1, len(ao) - x)])
return np.inner(survival_probabilities, discount_factors)
actuarial_pv(20, GMD_AO['2011'], 1.5)
def sensitivity_plot(gender, rz):
"""Sehr einfache Sensitivitätsanalyse"""
ao = DATEN[gender]['ao']
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
for colour, x in zip(['r', 'g', 'b', 'y', 'c', 'magenta'], [20, 30, 40, 50, 60, 70]):
xs = range(5)
ys = np.array([actuarial_pv(x, ao[year], rz) for year in COLUMNS])
ref_value = ys[0]
ys = (ys - ref_value)/ref_value * 100
cs = [colour] * len(xs)
ax.bar(xs, ys, zs=x, zdir='y', color=cs, alpha=0.8)
# ax.set_xlabel('X')
ax.set_title(f'Nur für Illustrationszwecke!\nRelative Barwert-Änderungen (in %) im Vergleich zu 2011,\n{DATEN[gender]["bez"]} bei Rechnungszins {rz}%')
ax.set_ylabel('Alter x')
ax.set_zlabel('Delta Barwerte (%)')
plt.xticks(np.arange(5), COLUMNS)
plt.savefig(GRAPHPATH+f'sens_{gender}_{rz}.pdf')
plt.show()
sensitivity_plot('m', 1.5)
sensitivity_plot('f', 1.5)
sensitivity_plot('m', 3.5)
Die klassische versicherungsmathematische Modellierung von biometrischen Risiken ist eigentlich sehr einfach:
Tatsächlich kann man mit geringem Aufwand etwas mehr Information gewinnen:
def mc_actuarial_pv(x, ao, rz, nof_persons, nof_scen):
"""Monte-Carlo-Simulation: Barwert einer lebenslänglichen nachschüssigen Leibrente ab Alter x"""
# Wahrscheinlichkeiten, das (x+k)-te Lebensjahr _nicht_ zu erreichen, für k = 1,2,3, ...
probabilities = np.array(1.0 - ao[x+1:]/ao[x])
# Rein finanzmathematische Barwerte von k jährlichen Zahlungen der Höhe 1, für k = 1,2,3, ...
pv = np.cumsum(np.array([0.0]+[(1+rz/100) ** (-t) for t in range(1, len(probabilities)+1)]))
# Erzeuge nof_scen gleichverteilte (Pseudo-)Zufallszahlen
scenarios = uniform.rvs(size=nof_scen*nof_persons)
# Simuliere nof_scen-mal eine Risikogemeinschaft von nof_persons x-Jährigen
# und bestimme die mittlere Auszahlung einer Leibrente der Höhe 1 für jedes Szenario
return np.sort(
np.mean(
np.array([pv[np.searchsorted(probabilities, s)] for s in scenarios]).reshape((nof_scen, nof_persons)),
axis=(1,)
)
)
mc_actuarial_pv(20, DATEN['m']['ao']['2011'], 1.5, 10, 100)
# Erwarteter Barwert:
actuarial_pv(20, DATEN['m']['ao']['2011'], 1.5)
# "Versicherungstechnische" Gewinnne/Verluste
def guv(x, ao, rz, nof_persons, nof_scens):
"""Bestimme GuV-Vektor als Differenz von vers.math. Barwert und simulierten Auszahlungen"""
return mc_actuarial_pv(x, ao, rz, nof_persons, nof_scens) - actuarial_pv(x, ao, rz)
def plot_risk_profile(x, ao, rz, nof_persons, nof_scens):
"""Plot von Monte-Carlo-Simulierten versicherungstechnischen Verlusten/Gewinnen"""
# Bestimme GuV-Vektor
guv_vector = guv(x, ao, rz, nof_persons, nof_scens)
# Bestimme (empirische) Wahrscheinlichkeitsverteilung
v_len = len(guv_vector)
cumul_probs = np.array([p/v_len for p in range(v_len)])
risk_profile = pd.DataFrame(cumul_probs, index=guv_vector, columns=['cumul. prob.'])
# Bestimme das 1%-Quantil (Value at Risk)
VaR = guv_vector[nof_scens // 100]
# Plot der Wahrscheinlichkeitsverteilung
myplot = risk_profile.plot()
myplot.set_xlabel('Verluste / Gewinne')
myplot.set_ylabel('Wahrscheinlichkeit')
myplot.set_title(f'GuV für {nof_persons} {x}-Jährige ({nof_scens} MC-Szenarien),\nVaR = {VaR:3.3}.')
plt.savefig(GRAPHPATH+f'mc_{nof_persons}_{x}_{nof_scens}.eps')
# return guv_vector
plot_risk_profile(20, DATEN['m']['ao']['2011'], 1.5, 10, 500)
plot_risk_profile(20, DATEN['m']['ao']['2011'], 1.5, 100, 500)
plot_risk_profile(20, DATEN['m']['ao']['2011'], 1.5, 1000, 500)
Ein bißchen weit hergeholt: Die Unterscheidung bei den Todesfällen in die Gruppen "0-60" Jahre und "61-100 Jahre" (siehe Fourier-Transformation) kann man als "Clustering von Datenwolken" interpretieren.
# xls_writer = pd.ExcelWriter(DATAPATH+'dummydaten.xlsx')
# Erzeuge dreidimensionale Zufalls-Daten
NOF_POINTS = [400, 500, 450]
SHIFT = np.array([65, 7, 9])
MEAN1 = np.array([1,0,0])*2 + SHIFT
MEAN2 = np.array([0,1,0])*2 + SHIFT
MEAN3 = np.array([0,0,1])*2 + SHIFT
COV1 = np.array([
[1.0, 0.5, -0.6],
[0.5, 2.0, -0.2],
[-0.6,-0.2, 3.0]
])/10
COV2 = np.array([
[3.0, 0.5, -0.6],
[0.5, 2.0, -0.2],
[-0.6,-0.2, 1.0]
])/10
COV3 = np.array([
[2.0, 0.5, -0.6],
[0.5, 2.0, -0.2],
[-0.6,-0.2, 2.0]
])/10
dummy = np.concatenate(
((np.random.multivariate_normal(MEAN1, COV1, NOF_POINTS[0])),
(np.random.multivariate_normal(MEAN2, COV2, NOF_POINTS[1])),
(np.random.multivariate_normal(MEAN3, COV3, NOF_POINTS[2]))
),
axis=0
)
dummydf = pd.DataFrame(dummy, columns=['Alter', 'Vertragsdauer', 'Rating'])
dummydf.head()
# dummydf.to_excel(xls_writer, sheet_name='Dummy 3D', index_label="Nr.", freeze_panes=(1,1))
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(dummy.T[0], dummy.T[1], dummy.T[2], s = 0.5)
ax.set_title('Kunden-Punktwolke')
ax.set_xlabel('Alter')
ax.set_ylabel('Vertragsdauer')
ax.set_zlabel('Internes Rating')
plt.savefig(GRAPHPATH+'Punktwolke3D.eps')
plt.show()
dummy = scale_this(dummy)
centroids = kmeans(dummy, 3)[0]
centroids
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(dummy.T[0], dummy.T[1], dummy.T[2], s=0.5)
ax.set_title('Kunden-Punktwolke, reskaliert\nmit Cluster-Zentren')
ax.scatter(centroids.T[0], centroids.T[1], centroids.T[2], s=30.0, c='red')
ax.set_xlabel('Alter (reskaliert)')
ax.set_ylabel('Vertragsdauer (reskaliert)')
ax.set_zlabel('Internes Rating (reskaliert)')
plt.savefig(GRAPHPATH+'Punktwolke3D_rescaled.eps')
plt.show()
Auch in "nur" drei Dimensionen wird die Visualisierung schon schwieriger.
# Erzeuge zweidimensionale Zufalls-Daten
NOF_POINTS = [3000, 5000, 5500, 3000]
SHIFT = np.array([62, 14])
MEAN1 = np.array([3,0])*0.95 + SHIFT
MEAN2 = np.array([0,1])*0.95 + SHIFT
MEAN3 = np.array([0,0]) + SHIFT
MEAN4 = np.array([3,1]) + SHIFT
COV1 = np.array([
[4.0, 0.5],
[0.5, 1.0]
])/6
COV2 = np.array([
[3.0, 0.5],
[0.5, 1.0]
])/6
COV3 = np.array([
[1.0, -0.5],
[-0.5, 1.0]
])/6
COV4 = np.array([
[3.0, -0.2],
[-0.2, 0.7]
])/6
dummy = np.concatenate(
((np.random.multivariate_normal(MEAN1, COV1, NOF_POINTS[0])),
(np.random.multivariate_normal(MEAN2, COV2, NOF_POINTS[1])),
(np.random.multivariate_normal(MEAN3, COV3, NOF_POINTS[2])),
(np.random.multivariate_normal(MEAN4, COV4, NOF_POINTS[3]))
),
axis=0
)
dummydf = pd.DataFrame(dummy, columns=['Alter', 'Rating'])
dummydf.head()
# dummydf.to_excel(xls_writer, sheet_name='Dummy 2D', index_label="Nr.", freeze_panes=(1,1))
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.scatter(dummy.T[0], dummy.T[1], s=0.3)
ax.set_title('Kunden-Punktwolke.')
ax.set_xlabel('Alter')
ax.set_ylabel('Internes Rating')
plt.savefig(GRAPHPATH+'Punktwolke2D.eps')
dummy = scale_this(dummy)
centroids = kmeans(dummy, 3)[0]
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.scatter(dummy.T[0], dummy.T[1], s=0.3)
ax.scatter(centroids.T[0], centroids.T[1], s=30.0, c='red')
ax.set_title('Kunden-Punktwolke, reskaliert\nmit Cluster-Zentren.')
ax.set_xlabel('Alter (reskaliert)')
ax.set_ylabel('Internes Rating (reskaliert)')
plt.savefig(GRAPHPATH+'Punktwolke2D_rescaled_3.eps')
centroids = kmeans(dummy, 4)[0]
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.scatter(dummy.T[0], dummy.T[1], s=0.3)
ax.scatter(centroids.T[0], centroids.T[1], s=30.0, c='red')
ax.set_title('Kunden-Punktwolke, reskaliert\nmit Cluster-Zentren.')
ax.set_xlabel('Alter (reskaliert)')
ax.set_ylabel('Internes Rating (reskaliert)')
plt.savefig(GRAPHPATH+'Punktwolke2D_rescaled_4.eps')
# xls_writer.save()
Machine learning, deep learning ("neural networks"): Zum Teil ausgereifte Software-Bibliotheken vorhanden, mit entsprechendem Hintergrundwissen und etwas Programmierkenntnissen ist es also relativ leicht möglich, Datenanalyse-Projekte umzusetzen.