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

rng = np.random.default_rng(1)                  # fixed seed

def share_and_se(flags):
    # a share is the average of 0/1 flags, so its SE is sqrt(p (1 - p) / M)
    flags = np.asarray(flags, dtype=float)
    p = float(flags.mean())
    return p, math.sqrt(p * (1 - p) / flags.size)

for m in (100, 10000, 1000000):
    dice = rng.integers(1, 7, size=(m, 2))      # whole numbers 1 to 6
    share, se = share_and_se(dice.sum(axis=1) >= 10)
    print(f"M = {m:>9,}:  estimate {share:.4%}   SE {se:.4%}   "
          f"off by {(share - 1 / 6) / se:+.2f} SE")