#  Code for Stochastic Calculus class, Jonathan Goodman
#  http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc20/
#  StockSim.py

#  Use forward Euler Maruyama to simulate SDE paths related to geometric
#  Brownian motion.
#  Make a histogram of the results and compute some expected values
#  See Exercises 4 and 5 of Week 4 notes for details

#  The author gives permission for anyone to use this publicly posted 
#  code for any purpose.  The code was written for teaching, not research
#  or commercial use.  It has not been tested thoroughly and probably has
#  serious bugs.  Results may be inaccurate, incorrect, or just wrong.

# illustrate using the Python random number generator from numpy version 19

from   numpy.random import Generator, PCG64       #  numpy randon number generator routines
import numpy             as np
import matplotlib.pyplot as plt

#-------------------- The strategy simulator, an Ito integral  --------------------

def sim( T, mu, sig, n, S0, rg):
    """Simulate a GBM path up to time T using n time steps, return S_T
       T   = simulate from time t=0 to time t=T
       mu  = expected rate of return
       sig = volatility
       n   = number of time steps
       S0  = starting value of the path
       rg = instantiated random number generator
       return values:
       ST = path value at time T
    """
    
    dt = T/n    # time step
    mudt = mu*dt   # for the drift term in the time step
    sqdt = np.sqrt(dt)  # standard deviation of dW
    
    dW = rg.normal(loc = 0., scale = sqdt, size = n)
    
    S = S0
    for k in range(n):
        S = S + mudt*S + S*sig*dW[k]
        
    return S
    
        
#------------------ The main program ---------------------------------------------------

bg = PCG64(12345)        # instantiate a bit generator with seed 12345
rg = Generator(bg)       # instantiate a random number generator

#    Problem parameters

T   = .05           # final time for the strategy simulation
mu  = .3            # expected rate of retrn
sig = 1.5  # vol
S0  = 10.   # starting price

#    Computing parameters

n   = 200      # number of time steps
Np  = 50000      # number of paths
Nb  = 50      # number of bins for the pdf plots

S = np.zeros([Np])
for path in range(Np):
    S[path] = sim( T, mu, sig, n, S0, rg)
    
Smax = S.max()
Smin = S.min()
dS = ( Smax - Smin)/(Nb-1)
bc = np.zeros( Nb, dtype = int)
print("max is " + str(Smax) + ", min is " + str(Smin))

for path in range(Np):
    bin = int( ( S[path] - Smin )/dS )
    bc[bin] += 1
pdf = np.zeros( Nb)
for bin in range(Nb):
    pdf[bin] = bc[bin]/(dS*Np)
    
Sax = np.linspace( Smin, Smax-dS, Nb)


PlotFile = "StockSim.pdf"   # filename for plot file


#        Setop for the cdf plots

fig, ax = plt.subplots()                       # Create a figure containing a single axes.
pdfLabel = "n = {n:7.2f}"
pdfLabel = pdfLabel.format(n = n)
ax.plot( Sax, pdf, label=pdfLabel)        # add this cdf to the figure
title = "Final price PDF, T = {T:6.2f}, mu = {mu:6.2f}, sig = {sig:6.2f}"
title = title.format( T = T, mu = mu, sig = sig)
ax.set_title(title)
ax.set_ylabel('PDF')
ax.set_xlabel('S')
ax.legend()
    
ax.grid()                  #  details of the plot, the grid for readibility

plt.savefig(PlotFile)       # Save a good copy of the plot
plt.show()