#  Stochastic Calculus, Courant Institute, NYU, Fall 2012
#   http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc2012/

#  Assignment 3, simulate the distribution of the max of a Brownian motion.  See HW for line 62 correction.

#  (your name and email here)

#  Set parameters

T    = 20.                   #  The length of a path (needs a decimal point)
N    = 100                   #  The number of time steps
L    = 10000                 #  The number of paths to generate
B    = 3.                    #  Record paths that don't hit this level.

dt   = T/N                   #  Will be wrong if T doesn't have a decimal point.

Mh = array( 0, L)            #  The maxima for the paths, for a histogram
Wh = array( 0, L)            #  The maxima for the paths, for a histogram
Sh = array( 0, L)            #  Record the path ends that never touch B
NS = 0                       #  The number of paths recorded in Sh

b = sqrt(dt)

#        Loop over paths

for ( k in (1:L)) {


  W  = 0.                    #  Standard Brownian motion starts at 0.
  M  = W                     #  M = largest value for this path so far
  
  Z = rnorm(N)

  for (i in (1:N) ) {        #  The loop that creates the path

    
    W = W + sqrt(dt)*Z[i]
    M = max( M, W)
   }
  Mh[k] = M
  Wh[k] = W
  if ( M <= B ) {           #  If this path never hit B
    NS     = NS + 1         #  Record that you found a new path
    Sh[NS] = W              #  Record its end value
   }

 }


#------------------plotting----------------------
  
Nb   = 60                 # number of histogram bins
db   = .25                # bin width
bs   = 0.                 # first bin edge
counts = array(0, Nb)     # bin counts
bc     = array(0, Nb)     # bin centers, for plotting
Nr     = 0                # number of "rejects", values outside the hist range
for ( k in (1:L)){
  bin = as.integer( ( Mh[k] - bs)/db )  #  An integer, the bin number
  bi  = bin + 1                 #  The bin index = (bin #) + 1
  if ( ( bi >0) && ( bi <= Nb ) ) {
    counts[bi] = counts[bi] + 1
   }
  else {
    Nr = Nr + 1
   }
 }

Wd = array( 0., Nb)      #  The probabiliy desnity of W, extimated.

for ( bi in (1:Nb)){
  bc[bi] = bs + (bi - .5)*db
  Wd[bi] = counts[bi]/((L-Nr)*db)
 }
  
plot(bc, Wd)
grid()

setwd("/Users/jg/desktop/notes/StochCalc2012/codes")  #  set the directory
pdf( "Assignment3.pdf")                   # open and name the plot file
 
title = sprintf('Probability density for the maximum of a BM path')
PlotInfo  = sprintf('%8d paths, time step = %8.2f, T = %8.2f',
                      L,                      dt,        T)
xlabel = "Max"
ylabel = "estimated density"
plot(  bc, Wd,           #  plot estimated density vs bin centers
       main = title,     #  overall plot title
       sub  = PlotInfo,  #  subtitle has more info on the plot
       xlab = xlabel,    #  label for the x axis
       ylab = ylabel,    #  label for the y axis
       type = 'p',       #  draw symbols
       pch  = 12,        #  a kind of box
       col  = 'blue' )   #  blue 
       

dev.off( which = dev.cur() ) #  This closes the current "device", which ...
                             #  ... is the .pdf file.  Otherwise the file ..
                             #  ... is not saved (or something).