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

rng = np.random.default_rng(5)                  # fixed seed
start, mean_log, sd_log = 10000, 1.8, 0.2 * math.sqrt(30)
tasks = [("final value", ",.0f", lambda z: start * np.exp(mean_log + sd_log * z)),
         ("log of W / W0", ".4f", lambda z: mean_log + sd_log * z),
         ("chance of a loss", ".4f", lambda z: (mean_log + sd_log * z < 0) * 1.0)]

def plain(func, z):
    # one value per draw, SE = s / sqrt(M)
    x = func(z)
    return float(x.mean()), float(x.std(ddof=1) / math.sqrt(x.size))

def antithetic(func, z):
    # the two halves of a pair are not independent, so the pair is one value
    x = (func(z) + func(-z)) / 2
    return float(x.mean()), float(x.std(ddof=1) / math.sqrt(x.size))

for name, spec, func in tasks:
    a, se_a = plain(func, rng.standard_normal(10000))
    b, se_b = antithetic(func, rng.standard_normal(5000))
    print(f"{name:<17} plain {a:{spec}} (SE {se_a:{spec}})   antithetic "
          f"{b:{spec}} (SE {se_b:{spec}})")

for name, spec, func in tasks:
    runs_a = [plain(func, rng.standard_normal(10000))[0] for _ in range(200)]
    runs_b = [antithetic(func, rng.standard_normal(5000))[0] for _ in range(200)]
    print(f"{name:<17} spread of 200 runs: plain {np.std(runs_a, ddof=1):{spec}}"
          f"   antithetic {np.std(runs_b, ddof=1):{spec}}")