# Stochastic Calculus, Courant Institute, NYU, Fall 2013 # http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc2013/ # Assignment 7, compute a value function using Monte Carlo # (your name and email here) #---------------------------------------------------------------------------- cosf <- function( x, a){ # Evaluate a function of x defined for -a <= x <= a. # In this case, the cosine function, V(x) = cos(pi*x/(2*a) return( cos(3.1415926535787*x/(2*a) ) ) } #----------------------------------------------------------------------------- # Simulate a path of n steps and return the value function. # Take gaussian steps for time step dt, return 0 if the path goes outside [-a,a] # Return V(X_T) otherwise. Assume T = n*dt. VXT <- function( x0, dt, n, a, V){ X = x0 # X is the location of the path at time t, initialize. Z = rnorm(n) # Standard normals, for time stepping, a vector of length n dW = sqrt(dt)*Z # Get the corresponding dW values using the time step. This is a vector instruction. for ( j in 1:n ) { X = X + dW[j] if ( abs(X) > a ) { return (0) } } return(V(X,a)) } #----------------------------------------------------------------------------- cat("Hello\n") T = 5 # The final time nt = 10 # The number of time steps up to T = t_n = n*dt dt = T/nt # see above a = 4 # width of the interval x0 = 0 # starting point Ns = 100000 # number of samples to evaluate f S = 0 # The sum of the sample values SSQ = 0 # If you are a statistician, this is obvious for ( k in 1:Ns){ # Take Ns independent samples Y = VXT(x0, dt, nt, a, cosf) S = S + Y SSQ = SSQ + Y*Y } EY = S/Ns EYSQ = SSQ/Ns sigsq = EYSQ - EY*EY cat("got mean = ", EY, ", and variance = ", sigsq, ", and error bar = ", sqrt(sigsq/Ns), "\n")