Das deutsche Python-Forum

Guten Abend,
brauche dringend Hilfe. Ich werde die Aufgabenstellung am ende des Textes einfügen.
Ich muss eine MCS durchführen, das standardverfahren habe ich hinbekommen. Nur muss ich es noch mit der Black-Scholes-Dynamik und lokalen Volatilität jeweils implementieren -.-
Da fehlt mir die Erfahrung :/ Und ich weiß nicht wie ich am besten zum Beispiel eine MCS mit 10.000 durchführen am besten analysieren und interpretieren kann mit einem geeignetem Plot.
Falls mir jemand helfen kann wäre es so nett ist wirklich dringend. Ich versuche auch mein Standard MCS hochzuladen, dass Sie sehen können das ich auch selber versuche was zu machen.

Schritt:
Erzeugung von mindestens 1.000, eher 2.000, 5.000 oder 10.000 zukünftigen Szenarien für den Aktienkurs und [ggf.] auch zugehörigen Volatilitätsszenarien mithilfe einer Monte Carlo Simulation (für den "morgigen" Tag) ausgehend vom heutigen Markt:
Dieser Teil ist eher stochastisch, aber kann auch elementar näherungsweise als eine Irrfahrt [1- oder 2-dimensionaler diskretisierter Wiener Prozesse] verstanden werden: Er bereitet insofern die empirischen Untersuchungen des zweiten Teils vor].

Vorab: Für alle Varianten schlage ich die Wahl von $\mu(t,S_t)= \mu_0 \cdot S_t$ mit einem festen Wert $\mu_0 \in \mathbb{R}$ vor, den man beispielsweise aus dem historischen Mittelwert der Log-Returns der Aktie bestimmen könnte.

(a) Variante 1: [klassische BSM-Dynamik]
Es wird $\sigma(t,S_t) \equiv \sigma_{t_0}$ konstant mit bspw. der historischen Volatilität des betrachteten Aktienkurses gewählt und die klassische Diskretisierung
$$
S_{(j+1) \cdot \Delta} = S_{j \cdot \Delta t} + \mu(t_j,S_{j \cdot \Delta t}) S_{j \cdot \Delta t}) \cdot \Delta t + \sigma_{t_0} S_{j \cdot \Delta t} \sqrt{\Delta t} \cdot Z_t
$$
mit einem pro Zeitschritt und Szenario unabhängig gezogenen standardnormalverteilten $Z \sim \mathcal{N}(0,1)}$ generiert: Die Volatilitätsoberfläche ist dann eine zu $(t,S)$-Ebene parallele Ebene auf der Höhe $S_{t_0}$.

(b) Variante 2: [lokale Volatilität]
Es wird $\sigma(t,S_t) := a_0 + a_{1,0} t + a_{0,1} S_t + \frac{1}{2} \left( a_{2,0} t^2 + 2 a_{1,1} t S_t + a_{0,2} S_t^2 \right) $ oder ganz einfach $\sigma(t,S_t) = \sigma_0 + sigma_1 t + \sigma_2 S_t + \sigma_3 S_t^2$ in der Simulation in jedem jeweiligen Zeitschritt bzw. Szenario verwendet, wobei die Aktienkursszenarien sich damit gemäß
$$
S_{(j+1) \cdot \Delta} = S_{j \cdot \Delta t} + \mu(t_j,S_{j \cdot \Delta t}) S_{j \cdot \Delta t}) \cdot \Delta t + \sigma\left( j \cdot \Delta t,S_{j \cdot \Delta t}\right) \sqrt{\Delta t} \cdot Z_t
$$
sich berechnen lassen. Die Parameter wie $\sigma_j$, $j \in \{0,1,2,3\}$, könnte man beispielsweise versuchen, an impliziten Volatilitätsdaten zu fitten.

import numpy as np
import pandas as pd
from pandas_datareader import data as wb
import matplotlib as mlp
import matplotlib.pyplot as plt
from scipy.stats import norm
import statsmodels.api as sm
import pylab

ticker = 'LHA.DE'
data = pd.DataFrame()
data[ticker] = wb.DataReader(ticker, data_source='yahoo', start='2018-4-1', end='2019-4-1')['Adj Close']

log_returns = np.log(1 + data.pct_change())
log_returns.tail()

u = log_returns.mean()
u

var = log_returns.var()
var

drift = u - (0.5 * var)
drift

stdev = log_returns.std()
stdev

type(drift)
type(stdev)

np.array(drift)

drift.values
stdev.values

norm.ppf(0.95)

x = np.random.rand(10, 2)
x
norm.ppf(x)

Z = norm.ppf(np.random.rand(10,2))
Z

## m=1000
t_intervals = 252
iterations = 1000

daily_returns = np.exp(drift.values + stdev.values * norm.ppf(np.random.rand(t_intervals, iterations)))
daily_returns

S0 = data.iloc[-1]
S0

price_list = np.zeros_like(daily_returns)
price_list
price_list[0]
price_list[0] = S0
price_list

for t in range(1, t_intervals):
price_list[t] = price_list[t - 1] * daily_returns[t]

price_list