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

#  Computing the Ito integral that represents a proportionate trading
#  rule in which the investor bets a fraction r of her assets at each
#  time on a standard Brownian motion.
#  See Exercises 4 and 5 of Week 2 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( dt, T, r, rbg):
    """Compute an Ito integral process up to time T with time step dt.
       dt = time step for the Ito integral
       r  = the fraction of assets invested at a given time
       T  = the end of the trading simulation time
       return values: X, W, isq
       X   = assets at time T
       W   = Brownian motion value at T
       isq = int from 0 to T of X_t^2
    """
    
    X = 1.                # start with a unit asset
    W = 0.                # start the Brownian motion path at 0
    t = 0.                # starting time
    sdt = np.sqrt(dt)     # avoid recomputing this
    isq = 0.              # accumulate the integral of the square
    while (1):
        dW  = sdt*rg.normal()     #  Generater the next Brownian motion increment
        W  += dW                  #  Add the increment to get W at time t+dt
        dX  = r*X*dW              #  The increment to the Ito integral
        X  += dX                  #  The Ito integral at time t+dt
        isq += (r*X)**2           #  the square of the integrand
        t += dt                   #  Done with work at this time
        if ( t > T ):             #  integrate (add up) terms while t<T
            isq = dt*isq          #  return values, get the integral
            return X, W, isq
        
#------------------ The main program ---------------------------------------------------

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

#    Collect all constants and parameters in one place

T        = 1.5           # final time for the strategy simulation
r        = .5            # the investment proportion
dtList   = [.4, .2, .1]  # values of Delta t to compare
Nsamp    = 10000         # number of Brownian motion paths for each dt
Nvals    = 500000        # number of bins for the cdf plots
dVal     = .00001        # bin size for the cdf plots

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


#        Setop for the cdf plots

fig, ax = plt.subplots()                       # Create a figure containing a single axes.
Xvals = np.linspace(0., dVal*Nvals, Nvals+1)   # X values, for the plot
for dt in dtList:                              # separate cdf for each dt

    counts = np.zeros(Nvals)                   # allocate arrays for the sim results
    cdf = np.zeros(Nvals+1)                    # the first bin of the cdf has zero
    for samp in range(Nsamp):                  # do the simulations
        X, W, isq = sim( dt, T, r, rg)         # catch the three return values
        bin = int( X/dVal)                     # a cdf is like the integral of the histogram
        if ( ( bin < Nvals ) and ( bin >= 0 ) ):   # ignore values out of range
            counts[bin] += 1                   # bin counts to make a histogram
    for bin in range(Nvals):                   # get the cdf by integrating the histogram
        cdf[bin+1] = cdf[bin] + counts[bin]
    cdf = cdf / Nsamp
    
    label = "Del t = {dt:10.2e}"               # label the cdf curve for this dt
    label = label.format(dt = dt)
    ax.plot( Xvals, cdf, label = label)        # add this cdf to the figure
    
ax.grid()                  #  details of the plot, the grid for readibility
ax.set_ylabel('cdf')       # all plots have axis labels, even if they're simple
ax.set_xlabel('X')
ax.legend()                #  multiple curves in one figure require a legend

title = "Wealth outcome, T = {T:6.2f}, r = {r:7.3f}"
title = title.format( T= T, r = r)
ax.set_title(title)
plt.savefig(PlotFile)       # Save a good copy of the plot
plt.show()