# QuantFinanceToGo, module 16: Monte Carlo Simulation. Solution to exercise 3.
# Runs on its own with Python 3 and numpy.
import math
from statistics import NormalDist
import numpy as np

start, years, paths, g, sigma = 10000, 30, 10000, 0.06, 0.2
rng = np.random.default_rng(3)                  # fixed seed
value = start * np.exp(rng.normal(g, sigma, (paths, years)).cumsum(axis=1))

def exact_percentile(t, u):
    # the log grows by g t on average and spreads by sigma sqrt(t)
    z = NormalDist().inv_cdf(u)
    return start * math.exp(g * t + z * sigma * math.sqrt(t))

def exact_loss(t):
    # below the start means ln(W / W0) below 0
    return NormalDist().cdf(-g * t / (sigma * math.sqrt(t)))

print(" t   10th sim / exact      median sim / exact     loss sim / exact")
for t in (1, 5, 10, 20, 30):
    lo, mid = np.percentile(value[:, t - 1], [10, 50])
    print(f"{t:>2}  {lo:>8,.0f} / {exact_percentile(t, 0.1):>8,.0f}"
          f"   {mid:>8,.0f} / {exact_percentile(t, 0.5):>8,.0f}"
          f"   {(value[:, t - 1] < start).mean():6.2%} / {exact_loss(t):6.2%}")

for t in (1, 4, 9, 16, 30):
    hi, lo = np.percentile(value[:, t - 1], [90, 10])
    print(f"t = {t:>2}: 90th / 10th = {hi / lo:5.2f}   exact "
          f"{exact_percentile(t, 0.9) / exact_percentile(t, 0.1):5.2f}")