# Stochastic Calculus, Courant Institute, NYU, Fall 2012
# http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc2012/
# Assignment 2, simulations of an urn process
# (your name and email here)
# Do simulations
m = 50 # The number of balls in the urn
T = 50 # The number of Morkov chain steps
L = 4000 # The number of paths to generate
x0 = m/2 # The starting state
x0 = as.integer(round(x0)) # make it an integer
p = .3 # Prob of replacement with a red ball
h = array(0, (m+1)) # Histogram of X_n
ia = 1:(m+1) # Array of i values, for plotting
a = 1:(m+1) # Arrays to hold transition probabilities
b = 1:(m+1) # a = down prob, c = up prob, b = don't move prob
c = 1:(m+1)
for ( i in ( 1: (m+1) )) {
a[i] = .4 # Put in the correct formulas here
c[i] = .2
b[i] = 1. - a[i] - c[i]
}
b[1] = b[1] + a[1] # Cannot go down from i=0 ...
a[1] = 0
b[m+1] = b[m+1] + c[m+1] # ... or up from m.
c[m+1] = 0
one = as.integer(1) # So that X is always an integer
for ( k in (1:L) ){ # generate path k
X = x0 # start at x0
for (n in (1:T)) { # simulate X_n
U = runif( 1, 0., 1.) # U is uniform in [0,1]
if ( U < a[X+1] && X > 0) X = X - one # Move left with probability = a
if ( U > 1.-c[X+1] && X < m) X = X + one # Move right with probability = c
} # end of path generation loop
h[X+1] = h[X+1]+1 # h[i+1] records that X_T = i
}
# Calculate the occupation probabilities
un = array(0, (m+1)) # The probabilities at time n
unp1 = array(0, (m+1)) # Probabilities at the next time
un[x0+1] = 1 # Start at x0, un[i]=0, all other i
for (n in (1:T) ){ # compute unp1
for ( i in ( 2:m)) { # compute the interior values first
unp1[i] = c[i-1]*un[i-1] + b[i]*un[i] + a[i+1]*un[i+1]
}
unp1[1] = b[1]*un[1] + a[2]*un[2] # the boundary cases
unp1[m+1] = c[m]*un[m] + b[m+1]*un[m+1]
un = unp1 # done with the time step
}
Ecounts = un*L # Expected bin count = prob*L
#-------------graphics---------------------------------
setwd("/Users/jg/desktop/notes/StochCalc2012/codes") # set the directory
pdf( "Assignment2.pdf") # open and name the plot file
title = sprintf('Histogram, urn process, %8d paths, %6d steps, p = %8.2f', L, T, p)
PlotInfo = "more info here"
xlabel = "i"
ylabel = "counts"
plot( ia, Ecounts, # plot Expected bin counts vs. i
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 = 'l', # draw a line (not just symbols)
lty = 2, # draw a dashed line
col = 'blue' ) # a blue line
lines( ia, h, # plot the bin counts vs i on the same plot
type = 'p', # draw symbols
pch = 12, # symbol 12 is a + in a box
col = 'black' ) # black symbols
grid() # dashed grid lines make the plot easier to read
# Create a legend saying what is in the graph
labels = c("Expected bin counts", "bin counts") # a list of two strings
# Where to put it
mh = max(h)
legend( .6*m, .9*mh, # Where to put the legend box
labels, # list of labesls for the lines
col = c("blue", "black"), # corresponding colors
lty = c( 2 , 0 ), # line types. 2 = dashed, 0 = no line
pch = c( NA , 12 ) ) # characters. NA means "none"
dev.off( which = dev.cur() ) # This closes the current "device", which ...
# ... is the .pdf file. Otherwise the file ..
# ... is not saved (or something).