# QuantFinanceToGo, module 16: Monte Carlo Simulation. The lesson program.
# Runs on its own with Python 3 and numpy. Every path in it is simulated.
# For learning only, nothing here is investment advice.
import math
from statistics import NormalDist
import numpy as np

bell = NormalDist()                     # the standard bell curve of module 7
start, years, paths = 10000, 30, 10000
mu, sigma = 0.08, 0.20                  # assumed average and volatility a year
g = mu - sigma ** 2 / 2                 # growth rate of module 10
t_sd = sigma * math.sqrt(years)         # sd of ln(W/W0) after 30 years

# 1) The model of one path
print("1) The model: yearly log returns from a bell with mean g and sd sigma")
print(f"   g = {mu} - {sigma}^2 / 2 = {g:.2f}   after {years} years: mean of "
      f"ln(W/W0) {g * years:.1f}, sd {t_sd:.4f}")

# 2) The simulation: 10,000 paths of 30 years
rng = np.random.default_rng(17)         # fixed seed, same paths on every run
log_ret = rng.normal(g, sigma, (paths, years))
value = start * np.exp(log_ret.cumsum(axis=1))
final = value[:, -1]
low, mid, high = np.percentile(final, [10, 50, 90])
print("2) The simulation: 10,000 simulated futures of $10,000")
print(f"   median ${mid:,.0f}   10th percentile ${low:,.0f}   90th ${high:,.0f}")
print(f"   chance of a loss {(final < start).mean():.2%} "
      f"({int((final < start).sum())} paths)   average ${final.mean():,.0f}")
print(f"   paths below the average: {(final < final.mean()).mean():.2%}")

# 3) The exact answers of the same model
z90 = bell.inv_cdf(0.9)
exact_median = start * math.exp(g * years)
exact_mean = start * math.exp(mu * years)
print("3) Exact answers of the model (no random numbers)")
exact_low = exact_median * math.exp(-z90 * t_sd)
exact_high = exact_median * math.exp(z90 * t_sd)
print(f"   median ${exact_median:,.0f}   10th ${exact_low:,.0f}   "
      f"90th ${exact_high:,.0f}")
print(f"   chance of a loss {bell.cdf(-g * years / t_sd):.2%}   average "
      f"${exact_mean:,.0f}   below the average {bell.cdf(t_sd / 2):.2%}")

# 4) Standard error and 95 % interval of the average
s = final.std(ddof=1)
se = s / math.sqrt(paths)
sd_exact = exact_mean * math.sqrt(math.exp(t_sd ** 2) - 1)
print("4) The average and its standard error SE = s / sqrt(M)")
print(f"   s ${s:,.0f}   SE ${se:,.0f}   95 % interval "
      f"${final.mean() - 1.96 * se:,.0f} to ${final.mean() + 1.96 * se:,.0f}")
print(f"   exact sd ${sd_exact:,.0f}   exact SE ${sd_exact / math.sqrt(paths):,.0f}"
      f"   estimate minus exact: {(final.mean() - exact_mean) / se:+.2f} SE")
for m in (100, 40000, 1000000):
    print(f"   with {m:>9,} paths the exact SE is ${sd_exact / math.sqrt(m):,.0f}")

# 5) The chance of a loss and its standard error
p_loss = (final < start).mean()
se_p = math.sqrt(p_loss * (1 - p_loss) / paths)
print("5) The chance of a loss as an average of 0/1 flags")
print(f"   {p_loss:.2%} with SE {se_p:.3%}: 95 % interval "
      f"{p_loss - 1.96 * se_p:.2%} to {p_loss + 1.96 * se_p:.2%}")

# 6) The error of the median: groups of the same 10,000 paths
print("6) Medians of groups: sd of their logs, theory sqrt(pi/2) * sd / sqrt(n)")
for n in (100, 1000):
    med = np.median(final.reshape(paths // n, n), axis=1)
    spread = np.log(med).std(ddof=1)
    theory = math.sqrt(math.pi / 2) * t_sd / math.sqrt(n)
    print(f"   {paths // n:>3} groups of {n:>5,}: spread {spread:6.2%}   theory "
          f"{theory:6.2%}   from ${med.min():,.0f} to ${med.max():,.0f}")

se_med =math.sqrt(math.pi / 2) * t_sd / math.sqrt(paths)
print(f"   all 10,000: {math.log(mid / exact_median):+.2%} off the exact median, "
      f"theory {se_med:.2%}")

# 7) Seeds: the same seed repeats, other seeds scatter
print("7) The same seed twice, then the seeds 1 to 5")
again = np.random.default_rng(17)
again_log = again.normal(g, sigma, (paths, years))
again_final = start * np.exp(again_log.cumsum(axis=1))[:, -1]
print(f"   seed 17 again: median ${np.median(again_final):,.0f}, every path the "
      f"same: {bool(np.array_equal(again_final, final))}")
seeds = [np.random.default_rng(1), np.random.default_rng(2), np.random.default_rng(3),
         np.random.default_rng(4), np.random.default_rng(5)]
for number, r in enumerate(seeds, start=1):
    f = start * np.exp(r.normal(g, sigma, (paths, years)).cumsum(axis=1))[:, -1]
    off = math.log(np.median(f) / exact_median) / se_med
    print(f"   seed {number}: median ${np.median(f):,.0f} ({off:+.2f} SE)   "
          f"average ${f.mean():,.0f}   loss {(f < start).mean():.2%}")

# 8) Antithetic pairs: every draw z is used again as -z
z = rng.standard_normal((paths // 2, years))
up = start * np.exp((g + sigma * z).sum(axis=1))
down = start * np.exp((g - sigma * z).sum(axis=1))
pair = (up + down) / 2
rho = np.corrcoef(up, down)[0, 1]
rho_exact = -math.exp(-t_sd ** 2)        # correlation of e^(a + bZ) and e^(a - bZ)
anti_exact = sd_exact * math.sqrt((1 + rho_exact) / 2) / math.sqrt(paths // 2)
print("8) 5,000 antithetic pairs, 10,000 paths in all")
print(f"   average ${pair.mean():,.0f}   SE "
      f"${pair.std(ddof=1) / math.sqrt(paths // 2):,.0f}   correlation in a pair "
      f"{rho:+.3f}, exact {rho_exact:+.3f}")
print(f"   exact SE: plain ${sd_exact / math.sqrt(paths):,.0f}, antithetic "
      f"${anti_exact:,.0f}")

# 9) The chance of a loss after 1, 10 and 30 years
print("9) Below the start after t years: simulated and exact")
for t in (1, 10, 30):
    print(f"   t = {t:>2}: {(value[:, t - 1] < start).mean():.2%}   exact "
          f"{bell.cdf(-g * t / (sigma * math.sqrt(t))):.2%}")