# 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).