##########################################################################################################################################################
### Program: SimTrade_mu_versus_sigma_Python_code_2026_01_04_SB                                                                                             ###
### Author: Saral Bindal                                                                                                                               ###
### Contact: saralbindal.24@kgpian.iitkgp.ac.in                                                                                                        ###
### Article on the SimTrade blog: https://www.simtrade.fr/blog_simtrade/modeling-asset-prices-financial-markets-arithmetic-geometric-brownian-motions/ ###
##########################################################################################################################################################


# Python program to simulate market prices using arithmetic Brownian motion and geometric Brownian motion.

# Objectives of the program:
#   1) Simulate market prices using ABM and GBM.
#   2) Calculate upper and lower line for confidence intervals.


####################################################################################################
### Organization of the program:                                                                 ###
###   STEP 1: Environment setup and package loading                                              ###
###   STEP 2: Define parameters                                                                  ###
###   STEP 3: Computation and plotting of price paths for ABM                                    ###
###   STEP 4: Computation and plotting of price paths for GBM                                    ###
####################################################################################################


####################################################################################################
### Parameters you can change for BSM option pricing:                                            ###
###   n_sims   : Number of simulations                                                           ###
###   n_months : Number of time steps                                                            ###
###   dt       : time step                                                                       ###
###   mu       : Annualized drift/mean                                                           ###
###   sigma    : Annualized volatility/standard deviation                                        ###
###   S0       : Initial market price                                                            ###
####################################################################################################


#################
# Documentation #
#################

# Python packages
# https://numpy.org/doc/stable/
# https://matplotlib.org/stable/
# https://docs.scipy.org/doc/scipy/

# Academic articles
# Bachelier, L. (1900). Theory of Speculation. Annals of the Scientific School of the École Normale Supérieure, 3rd series, 17, 21–86.
# Samuelson P. A. (1965). Rational theory of warrant pricing. Industrial Management Review, 6(2), 13–39.
# Wiener N. (1923). Differential-space. Journal of Mathematics and Physics, 2, 131–174.


#################################################
# STEP 1: Environment setup and package loading #       
#################################################

import sys
import subprocess

# Function to install a package if it is not already installed
def install(package):
    subprocess.check_call([sys.executable, "-m", "pip", "install", package])

# Import numpy (numerical computations)
try:
    import numpy as np
except ImportError:
    install("numpy")
    import numpy as np

# Import matplotlib (data visualization)
try:
    import matplotlib.pyplot as plt
except ImportError:
    install("matplotlib")
    import matplotlib.pyplot as plt

# Import scipy.stats (Normal distribution)
try:
    from scipy.stats import norm
except ImportError:
    install("scipy")
    from scipy.stats import norm

##############################
# STEP 2: Define parameters  #
##############################

n_sims = 100000   # Number of simulations
n_months = 120    # Number of time steps
dt = 1/12         # time step

###########################################################
# STEP 3: Computation and plotting of price paths for ABM #
###########################################################

# 1. ARITHMETIC BROWNIAN MOTION
m = 8            # mu (in $)
s = 15           # sigma (in $)
S0 = 100         # Initial market price (ABM can take negative values)

# GENERATE RANDOM SHOCKS (MONTE CARLO)
# We create a 2D matrix: columns = simulations, rows = days
# Each entry represents a Brownian increment over one time step
increments = np.random.normal(m * dt, s * np.sqrt(dt), size=(n_months, n_sims))

# CALCULATE PRICE PATHS AND CONFIDENCE BANDS
# ABM evolves additively via cumulative sums of increments
price_paths = S0 + np.cumsum(increments, axis=0)

# Time vector (in years)
t = np.arange(1, n_months + 1) * dt
months = np.arange(1, n_months + 1)

# Confidence level
alpha = 0.66
z_alpha = norm.ppf((1 + alpha) / 2)

# Theoretical values
mean_path = S0 + m * t
upper_band = S0 + m * t + z_alpha * s * np.sqrt(t)
lower_band = S0 + m * t - z_alpha * s * np.sqrt(t)

# VISUALIZATION
# Plot only the first 1000 paths for visual clarity
plt.figure(figsize=(12, 6))
plt.plot(months, price_paths[:, :1000], color='black', alpha=0.02)

plt.plot(months, mean_path, color='cyan', linestyle='--', linewidth=2, label='Theoretical Mean')

# Add confidence interval as a shaded region
plt.fill_between(months, lower_band, upper_band, color='cyan', alpha=0.2, label='66% Confidence Interval')

plt.title(f"Monte Carlo Simulation: {n_sims} SBM Price Paths")
plt.xlabel("Time (in months)")
plt.ylabel("Asset price ($)")
plt.grid(True, alpha=0.3)
plt.legend()

plt.xlim(0, 120)
plt.xticks(range(0, 121, 10))
plt.ylim(0, 400)
plt.yticks(range(0, 401, 100))

plt.show()



###########################################################
# STEP 4: Computation and plotting of price paths for GBM #
###########################################################

#Parameters
mu = 0.08         
sigma = 0.15     

# 2. GEOMETRIC BROWNIAN MOTION
m = mu - 0.5 * sigma**2    # Drift adjusted for log-normality
s = sigma
S0 = 100                   # Initial market price (GBM remains strictly positive)

# GENERATE RANDOM SHOCKS (MONTE CARLO)
# We create a 2D matrix: columns = simulations, rows = days
# These represent log-returns
log_returns = np.random.normal(m * dt, s * np.sqrt(dt), size=(n_months, n_sims))

# CALCULATE PRICE PATHS
# Cumulative sum of log-returns followed by exponentiation
# Ensures prices stay positive
price_paths = S0 * np.exp(np.cumsum(log_returns, axis=0))

# Time vector (in years)
t = np.arange(1, n_months + 1) * dt
months = np.arange(1, n_months + 1)

# Confidence level
alpha = 0.66
z_alpha = norm.ppf((1 + alpha) / 2)

# Theoretical values
mean_path = S0 * np.exp(mu * t)
upper_band = S0 * np.exp((mu - 0.5 * sigma**2) * t + z_alpha * sigma * np.sqrt(t))
lower_band = S0 * np.exp((mu - 0.5 * sigma**2) * t - z_alpha * sigma * np.sqrt(t))

# VISUALIZATION
# Plot only a subset of paths to avoid overplotting
plt.figure(figsize=(12, 6))
plt.plot(months, price_paths[:, :1000], color='black', alpha=0.02)

plt.plot(months, mean_path, color='cyan', linestyle='--', linewidth=2, label='Theoretical Mean')

# Add confidence interval shading
plt.fill_between(months, lower_band, upper_band, color='cyan', alpha=0.2, label='66% Confidence Interval')

plt.title(f"Monte Carlo Simulation: {n_sims} GBM Price Paths")
plt.xlabel("Time (in months)")
plt.ylabel("Asset price ($)")
plt.grid(True, alpha=0.3)
plt.legend()

plt.xlim(0, 120)
plt.xticks(range(0, 121, 10))
plt.ylim(0, 400)
plt.yticks(range(0, 401, 100))

plt.show()


# End of the code