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

rng = np.random.default_rng(4)                  # fixed seed
start, mean_log, sd_log = 10000, 1.8, 0.2 * math.sqrt(30)
exact_final = start * math.exp(0.08 * 30)

def interval(values):
    # the central limit theorem makes the average about a bell around the truth
    values = np.asarray(values, dtype=float)
    se = values.std(ddof=1) / math.sqrt(values.size)
    return values.mean() - 1.96 * se, values.mean() + 1.96 * se

for m in (100, 1000):
    hit_final = hit_log = 0
    for _ in range(1000):
        logs = rng.normal(mean_log, sd_log, m)
        low, high = interval(start * np.exp(logs))
        hit_final += low <= exact_final <= high
        low, high = interval(logs)
        hit_log += low <= mean_log <= high
    print(f"M = {m:>5,}: the interval of the final value holds it in "
          f"{hit_final / 1000:.1%}, the interval of the log in "
          f"{hit_log / 1000:.1%}")