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

start, years, g, sigma = 10000, 30, 0.06, 0.2
rng_a = np.random.default_rng(2)                # two generators, same seed
rng_b = np.random.default_rng(2)

def median_of_run(rng, paths):
    # every call takes the next random numbers of rng: a new, independent run
    log_ret = rng.normal(g, sigma, (paths, years))
    final = start * np.exp(log_ret.sum(axis=1))
    return float(np.median(final))

first, twin = median_of_run(rng_a, 1000), median_of_run(rng_b, 1000)
print(f"seed 2 in two generators: ${first:,.2f} and ${twin:,.2f}")
print(f"the first generator again: ${median_of_run(rng_a, 1000):,.2f}")
logs = np.log([median_of_run(rng_a, 1000) for _ in range(100)])
theory = math.sqrt(math.pi / 2) * sigma * math.sqrt(years) / math.sqrt(1000)
print(f"100 runs of 1,000 paths: medians from ${math.exp(logs.min()):,.0f} "
      f"to ${math.exp(logs.max()):,.0f}")
print(f"spread of their logs {logs.std(ddof=1):.2%}   theory {theory:.2%}")