# QuantFinanceToGo, module 16: Monte Carlo Simulation. Exercise 5 (hard).
# Runs on its own with Python 3 and numpy. Fill in the lines marked TODO.
#
# Antithetic pairs: every standard bell draw z is used twice, as z and -z.
# In the model of step 1, W / W0 = exp(1.8 + 0.2 sqrt(30) z) after 30 years.
# a) Write plain(func, z): the average of func(z) and its standard error.
# b) Write antithetic(func, z): the average of the pair values
#    (func(z) + func(-z)) / 2 and their standard error, one value per pair.
# c) Compare both with the same number of paths, 10,000 plain draws against
#    5,000 pairs, for the final value, its log and the chance of a loss.
#    Then repeat both 200 times and measure the real spread of each estimate.
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):
    # TODO: (average of func(z), its standard error)
    return None

def antithetic(func, z):
    # TODO: pair values (func(z) + func(-z)) / 2, then (average, standard error)
    return None

line = tasks[1][2]
if plain(line, np.ones(4)) is None or antithetic(line, np.ones(4)) is None:
    print("Exercise 5 is not solved yet. Fill in the lines marked TODO.")
else:
    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}}")