Scientific Computing I/MAPL660/Fall 2001/Liu

Homework #2, due October 16

Warning: Monte–Carlo burns computer time. This assignment has been scaled to fit into a recent model Pentium computer or better. If it takes too long on your computer, scale it down a bit.

(1) Explore basic Monte–Carlo and central limit theorem.

Let X be the random variable X = max(T_{1}, T_{2}, T_{3}), where the T_{j} are three independent uniformly distributed random variables in the interval [0, 1].

(a) Show that the cumulative distribution function for X is F (x) = x^{3} and the probability density function for X is f(x) = 3·x^{2}. Construct 10^{4} samples of X using three uniformly distributed random numbers and make a histogram with a bin size .01. Compare the empirical histogram with the underlying density agree f(x). (You can use standard random generators such as ran1.c in Numerical Recipes, Section 7.1, rand in Matlab to generate the uniform random numbers and use hist in Matlab to make the histogram)

(b) Let S_{n} be the sum of n independent X samples:

S_{n} = ∑_{k=1}^{n} X_{k} ,

Let the brackets ⟨ ⟩ represent an ensemble average. Show that the mean value m_{n} is given by m_{n} = ⟨S_{n}⟩ = \frac{3}{4}n, and the variance, v_{n} = ⟨(S_{n} − m_{n})^{2}⟩ = \frac{3}{80}n. Average 10^{4} trajectory value of S_{n} at n = 5j, j = 1, …, 20 and plot the empirical mean and variance with the exact values m_{n} and v_{n}, respectively. (error bar can be ploted by errorbar in matlab)

(c) Construct a histogram of the S_{100} values from problem 1(b) above. Choose bin widths such that 100 equally spaced bins capture approximately 99% of data (from the central limit theorem, three standard deviations shall capture approximately 99% of the data). On the same graph, plot a Gaussian density function with the appropriate mean and variance, and compare it with the empirical histogram.

(2) Let T_{k} be independent uniform random variables in [0, 1]. We define a set of random variables {y_{k}} by

y_{0} = 0, y_{k} = \frac{1}{2}y_{k−1} + T_{k}, for k = 1, 2, ⋯  (1)

Let A_{k} represent a trial in which the value of y_{k} is greater than 1. Define S_{n} as follows

S_{n} = # {A_{k} | k ≤ n} .

(a) Conduct 1×10^{4}, 2×10^{4} and 5×10^{4} samples of S_{n} for n = 100, 200, and 500, respectively, and show that S_{n}/n is approximately normally distributed. Turn in a table showing mean and variance of S_{n}/n for n = 100, 200 and 500.

(b) Insight into the outcome of the simulations can be extracted from characterizing the probability distribution of y_{k}. Rewrite y_{k} in (1) in a more compact manner:

y_{k} = ∑_{j=0}^{k−1} \frac{1}{2^{j}} T_{j}

and show that the expected values of y_{k} are

⟨y_{k}⟩ = 1 − 2^{−k}

Therefore, one expects that, since the expected value of y_{k} is close to one, the expected number of s_{n}/n is close to 1/2.

(3) Explore Gaussian deviates and the Box–Muller algorithm.

We want to study whether the resulting random variables (X, Y) for the Box–Muller algorithm are members of the following joint probability density function:

f(x, y) = \frac{1}{2π} e^{−(x^{2}+y^{2})/2} ,  (2)

(a) The first task is to capture random variable pairs (X, Y) form the Box-Muller algorithm and determine how many “fall” within a pre-defined region on the xy-plane. We will consider only the central “bin” centered over the origin of the xy-plane. Let the xy-plane be subdivided into bins:

B_{jk} = \{(x, y) | |x − x_{j}| ≤ \frac{h}{2} and |y − y_{k}| ≤ \frac{h}{2}\} .

centered at x_{j} = jh and y_{k} = kh. If we take n pairs, (X_{t}, Y_{t}), t = 1, …, n, then the bin counts are

N_{jk} = # {(X_{t}, Y_{t}) ∈ B_{jk} | 1 ≤ t ≤ n} .

The expected counts are

< N_{jk} >= n ∫∫_{B_{jk}} f(x, y)dxdy ,  (3)

where f is given by (2). The midpoint rule for the integral in (3) is second order accurate. In terms of local truncation error, this means that

∫∫_{B_{jk}} f(x, y)dxdy = h^{2}f(x_{j}, y_{k}) + O(h^{4}) .

In two dimensions we do not want h so small because that would create too many bins. Therefore we want a more accurate integration formula. Expand f(x, y) in a Taylor series at the bin center (x_{j}, y_{k}) and show that

∫∫_{B_{jk}} f(x, y)dxdy = h^{2}f(x_{j}, y_{k}) + \frac{h^{4}}{24} (\frac{∂^{2}f}{∂x^{2}} + \frac{∂^{2}f}{∂y^{2}} ) (x_{j}, y_{k}) + O(h^{6}) .

Choose h so that for n = 10^{6}, the central bin, B_{00}, has on the order of a thousand points in it.

(b) Show that the moments:

⟨X^{2p}⟩ ≡ \frac{1}{\sqrt{2π}} ∫_{∞}^{∞} x^{2p}e^{−x^{2}/2}dx = 1 · 3 · 5 ⋯ (2p − 1)

What is the value of the variance for X^{2p}? For p = 1, 2, and 5, find out the values of n required to achieve an accuracy of 1% for the sampling moments

\frac{1}{n} ∑_{t=1}^{n} X_{t}^{2p}

approximate ⟨X^{2p}⟩. (In Monte Carlo simulation, a 66% confidence intervals (one standard derivation) is usually used in the error estimation.) Now call the independent Gaussian given by Box–Muller and estimate moments for p = 1, 2, 5. Compare with the above accuracy predication. Explain why more samples are required for large p?

(4) Explore Brownian motion simulation.

A standard Brownian motion can be simulated by taking

X_{0} = 0, X_{k} = X_{k−1} + \sqrt{Δt}Z_{k}, k = 1, 2, ⋯

where the Z_{k} are i.i.d. standard normals, X_{k} ≈ X(kΔt). The simulation will be performed for 10 steps, each time step is Δt = .1. We seek to calculate the probability P of a trajectory extending to X_{k} ≥ 3 for some k:

P = Pr (X_{k} ≥ 3 for some t_{k} ≤ 1) .

(a) What precision can be achieved using 10^{6} trajectories?

(b) Express P as an integral form and show how it can be computed using importance sampling. Let a trajectory be represented by the vector Z⃗:

Z⃗ = {z_{1}, z_{2}, ⋯, z_{10}}

The corresponding differential elements is dz⃗ = dz_{1}dz_{2} ⋯ dz_{10}. The local of the individual steps within the trajectory are represented by x⃗:

x⃗ = {x_{0}, x_{1}, ⋯, x_{10}}

Each element of x⃗ is defined as follows:

x_{0} = 0, x_{k} = x_{k−1} + \sqrt{Δt}z_{k}, k = 1, 2, ⋯

Let the decision function be ϕ(z⃗):

ϕ(z⃗) = \{ \begin{matrix}1 & if x_{k} ≥ 3 for some k ≤ 10 \\ 0 & otherwise. \end{matrix}

Show that

P = ∫ ϕ(z⃗)f(z⃗)dz⃗ = ∫ ⋯ ∫ ϕ(z_{1}, …, z_{10})f(z_{1}, …, z_{10})dz_{1} ⋯ dz_{10} ,

where

f(z⃗) = f(z_{1}, …, z_{10}) = \frac{1}{(2π)^{10/2}} exp \{−(z_{1}^{2} + ⋯ + z_{10}^{2})/2\} ; .

To implement importance sampling, we can choose a new distribution for the Z_{k} that selects random deviates with a greater chance of giving ϕ = 1 for a given trajectory. Let these new random variables be represented by Y_{k}. A simple choice is to take Y_{k} to be i.i.d. Gaussians with mean a and variance 1. Show that P can be rewritten as

P = ∫ ϕ(y⃗)h(y⃗)g(z⃗)dy⃗, h(y⃗) = \frac{f(y⃗)}{g(y⃗)}  (4)

where g represents the joint density for ten independent biased Gaussian random variables. Show that

h(y⃗) = e^{5a^{2}−a(y_{1}+y_{2}+···+y_{10})} ,

and the (4) can be implemented by

P ≈ \frac{1}{n} ∑_{j=1}^{n} ϕ(Y^{j})h(Y^{j}) ,

where the Y^{j} are each ten dimensional vectors.

(c) Choose a in g(y⃗) to achieve a precision of 1% using far fewer trajectories than required without importance sampling. There is no simple way to choose a. You will vary the value of a from a = 0, .2, .4, ⋯, 2 and perform the Monte Carlo simulation. For each a, list in a table the number of trajectories, empirical value of P, standard derivation, relative error. Find the optimum value of a and number of trajectories need to achieve 1% accuracy.

(5) Application of Brownian motion simulation to finance. We try to simulate the stochastic differential equation

dS = rSdt + σSdZ  (5)

for 0 ≤ t ≤ 1 with initial condition S(0) = 1 and with the parameters r = .08, σ = .08. This stochastic differential equation is not our model of actual movements of the stock price (that has μ instead of r). Rather, it is the “risk neutral martingale measure” that we use in valuation problems. To solve (5) by the forward Euler method is to generate (approximate) sample paths starting with S_{0} = S(0) = 1 with M = 1/Δt time steps using

S_{k+1} = S_{k} + rS_{k}Δt + σS_{k}\sqrt{Δt}Z_{k} ,  (6)

where S_{k} ≈ S(kΔt). The pay-out function, ϕ(S), for a European call option at a struck price S = 1, is defined as

ϕ(S) = max(S(1) − 1, 0) .

Use (6) to generate N realizations of S and compute its pay-out value (denoted as X_{j}). The excepted value of the option, after discounting the interest, in a risk-neutral world then can be computed by

c = e^{−r}⟨ϕ(S)⟩ ≈ e^{−r} \frac{1}{N} ∑_{j=1}^{N} X_{j} ,

The exact value for the option can be computed by the following Black-Scholes pricing formula (See J.C. Hull, Options, Futures, and Other Derivative Securities, Prentice Hall, 1989, Section 10.9).

c = F (d_{1}) − e^{−r}F (d_{2}), d_{1} = r/σ + σ/2, d_{2} = r/σ − σ/2

where F (x) is the cumulative density function for the standard normal N(0, 1). Compare your simulation with the above exact solution.