function s=simcev(n,r,sigma,So,T,gam)

% simulates n sample paths of a CEV process on the interval [0,T] all with

% the same starting value So. assume gamma != 1.

Yt=ones(n,1)*(So^(1-gam))/(1-gam); y=Yt;

dt=T/1000; c1=r*(1-gam)/sigma; c2=gam*sigma/(2*(1-gam));

dw=normrnd(0,sqrt(dt),n,1000);

for i=1:1000

v=find(Yt); % selects positive components of Yt for update

Yt=max(0,Yt(v)+(c1.*Yt(v)-c2./Yt(v))*dt+dw(v,i));

y=[y Yt];

end

196 CHAPTER 3. BASIC MONTE CARLO METHODS

s=((1-gam)*max(y,0)).^(1/(1-gam)); %transforms to St

For example when r = .05, σ = .2, ∆t = .00025, T = .25, γ = 0.8 we can

generate 1000 sample paths with the command

s=simcev(1000,.05,.2,10,.25,.8);

In order to estimate the price of a barrier option with a down-and-out barrier

at b and exercise price K, capture the last column of s,

ST=s(:,1001);

then value a European call option based on these sample paths

v=exp(-r*T)*max(ST-K,0);

¬nally setting the values equal to zero for those paths which breached the

lower barrier and then averaging the return from these 1000 replications;

v(min(s™)<=9)=0;

mean(v);

which results in an estimated value for the call option of around $0.86. Al-

though the standard error is still quite large (0.06), we can compare this with the

Black-Scholes price with similar parameters. [CALL,PUT] = BLSPRICE(10,10,.05,.25,.2,0)

which gives a call option price of $0.4615. Why such a considerable di¬erence?

Clearly the down-and-out barrier can only reduce the value of a call option.

Indeed if we remove the down-and-out feature, the European option is valued

closer to $1.28 so the increase must be due to the di¬erences betwen the CEV

process and the geometric Brownian motion. We can con¬rm this by simulating

the value of a barrier option in the Black_Scholes model later on.

PROBLEMS 197

Problems

1. Consider the mixed generator xn = (axn’1 + 1)mod(m) with m = 64.

What values of a results in the maximum possible period. Can you indicate

which generators appears more and less random?

2. Consider the multiplicative generator with m = 64. What values of

a result in the maximum possible period and which generators appears

more random?

3. Consider a shu¬„ed generator described in Section 3.2 with k = 3, m1 =

7, m2 = 11.

Determine the period of the shu¬„ed random number generator above and

compare with the periods of the two constituent generators.

4. Prove: If Xis random variable uniform on the integers {0, . . . , m ’ 1}

and if Y is any integer-valued random variable independent of X, then

the random variable W = (X + Y )mod m is uniform on the integers

{0, . . . , m ’ 1}.

5. Consider the above quadratic residue generator xn+1 = x2 mod m with

n

m = 4783 — 4027. Write a program to generate pseudo-random numbers

from this generator. Use this to determine the period of the generator

starting with seed x0 = 196, and with seed x0 = 400?

6. Consider a multiplicative generator of the form xn+1 = axn (mod m).

Prove that if m = m1 m2 is the product of relative primes (i.e. gcd(m1 , m2 ) =

1), then the period of the generator xn+1 = axn (mod m) is the least com-

mon multiple of the generators with moduli m1 and m2 .

7. Consider a sequence of independent U [0, 1] random variables U1 , ..., Un .

198 CHAPTER 3. BASIC MONTE CARLO METHODS

De¬ne indicator random variables

Si = 1 if Ui’1 < Ui and Ui > Ui+1 for i = 2, 3, ..., n ’ 1, otherwise Si = 0,

Ti = 1 if Ui’1 > Ui and Ui < Ui+1 for i = 2, 3, ..., n ’ 1, otherwise Ti = 0.

Verify the following:

(a)

X

R=1+ (Si + Ti )

(b)

2n ’ 1

1

and E(R) =

E(Ti ) = E(Si ) =

3 3

(c) cov(Ti , Tj ) = cov(Si , Sj ) = ’ 1 .if |i ’ j| = 1 and it equals 0 if

9

|i ’ j| > 1.

§

⎪ 5 1 7

’ |i ’ j| = 1

⎪ := 72 if

⎪ 24 9

⎨

’1

(d) cov(Si , Tj ) = .

if i=j

⎪9

⎪

⎪

©0 |i ’ j| > 1

if

(e) var(R) = 2(n’2) 1 ( 2 )+4(n’3)(’ 1 )+4(n’3)( 72 )+2(n’2)(’ 1 ) =

7

33 9 9

3n’5

18 .

(f) Con¬rm these formulae for mean and variance of R in the case n =

3, 4.

8. Verify that for the serial correlation statistic Cj,

§

⎪ 4/45n

⎪ for j=0

⎪

⎨

7

var(Cj ) = j ≥ 2, even

for

⎪ 72n

⎪

⎪

© 13

j ≥ 1, odd

for

144n

9. Consider the turbo-pascal generator xn+1 = (134775813xn + 1) mod 232 .

Generate a sequence of length 5000 and apply the serial correlation test.

Is there evidence of dependence?

PROBLEMS 199

10. Generate 1000 daily “returns” Xi , i = 1, 2, ..., 1000 from each of the two

distributions, the Cauchy and the logistic. In each case, assume that

a = 0, b = 0.1. Graph the total return over an n day period versus n. Is

there a qualitative di¬erence in the two graphs? Repeat with a graph of

the average daily return.

11. Brie¬‚y indicate an e¬cient algorithm for generating one random vari-

able from each of the following distributions and generate such a ran-

dom variable using one or more of the uniform[0,1] random numbers.

Ui : 0.794 0.603 0.412 0.874 0.268 0.990 0.059 0.112 0.395

(a) X∼ U [’1, 2].

3 1/2

(b) a random variable X with probability density function f (x) = 16 x ,0 <

x<4

(c) A discrete random number X having probability function P [X =

x] = (1 ’ p)x p, x = 0, 1, 2, ..., p = 0.3.

(d) A random variable X with the normal distribution, mean 1 and

variance 4.

(e) A random variable X with probability density function

f (x) = cx2 e’x , 0 · x < 1

for constant c = 1/(2 ’ 5e’1 ).

(f) A random variable X with the following probability function:

0 1 2 3

x

0.1 0.2 0.3 0.4

P [X = x]

12. Consider the multiplicative pseudo-random number generator

xn+1 = axn mod 150

200 CHAPTER 3. BASIC MONTE CARLO METHODS

starting with seed x0 = 7. Try various values of the multiplier a and

determine for which values the period of the generator appears to be max-

imal.

13. Consider the linear congruential generator

xn+1 = (axn + c) mod 28

What is the maximal period that this generator can achieve when c = 1

and for what values of a does this seem to be achieved? Repeat when

c = 0.

14. Evaluate the following integral by simulation:

Z1

(1 ’ x2 )3/2 dx.

0

15. Let U be a uniform random variable on the interval [0,1]. Find a function

of U which is uniformly distributed on the interval [0,2]. The interval

[a, b]?

16. Evaluate the following integral by simulation:

Z2

x3/2 (4 ’ x)1/2 dx.

0

17. Evaluate the following integral by simulation:

Z∞

2

e’x dx.

’∞

R∞ 2

(Hint: Rewrite this integral in the form 2 0 e’x dx and then change

variables to y = x/(1 + x))

18. Evaluate the following integral by simulation:

Z 1Z 1

2

e(x+y) dxdy.

0 0

(Hint: Note that if U1 , U2 are independent Uniform[0,1] random variables,

R1R1

E[g(U1 , U2 )] = 0 0 g(x, y)dxdy for any function g).

PROBLEMS 201

√ √

1 ’ 2x2 < y < 1 ’ 2x4 }

Figure 3.28: The region {(x, y); ’1 < x < 1, y > 0,

19. Find the covariance cov(U, eU ) by simulation where U is uniform[0,1]

and compare the simulated value to the true value.

20. Find by simulation the area of the region {(x, y); ’1 < x < 1, y > 0,

√ √

1 ’ 2x2 < y < 1 ’ 2x4 }. The boundaries of the region are graphed in

Figure 3.28.

21. For independent uniform random numbers U1 , U2,.... de¬ne the random

P

variable N = min{n; n Ui > 1}.

i=1

Estimate E(N ) by simulation. Repeat for larger and larger numbers of

simulations. Guess on the basis of these simulations what is the value of

E(N ). Can you prove this?

22. Assume that you have available a Uniform[0,1] random number generator.

Give a precise algorithm for generating observations from a distribution

with probability density function

(x ’ 1)3

f (x) =

4

for 1 · x · 3. Record the time necessary to generate the sample mean of

5,000 random variables with this distribution. .

202 CHAPTER 3. BASIC MONTE CARLO METHODS

23. Assume that you have available a Uniform[0,1] random number generator.

Give a precise algorithm for generating observations from a distribution

(x’20)

for 20 · x · 40. Record the time

with probability density function 200

necessary to generate the sample mean of 5,000 random variables with

this distribution. .

24. Assume that you have available a Uniform[0,1] random number generator.

Give a precise algorithm for generating observations from a distribution

with a density function of the form f (x) = cx3 e’x/2 for x > 0 and ap-

propriate constant c. Record the time necessary to generate the sample

mean of 100,000 random variables with this distribution.

25. Assume that you have available a Uniform[0,1] random number genera-

tor. Give a precise algorithm for generating observations from a discrete

distribution with P [X = j] = (2/3)(1/3)j ; j = 0, 1, ....Record the time

necessary to generate the sample mean of 100,000 random variables with

this distribution. .

26. Assume that you have available a Uniform[0,1] random number generator.

Give a precise algorithm for generating observations from a distribution

with probability density function f (x) = e’x , 0 · x < ∞. Record the

time necessary to generate the sample mean of 100,000 random variables

with this distribution. Compute as well the sample variance and compare

withe the sample mean. How large would the simulation need to be if

we wanted to estimate the mean within 0.01 with a 95% con¬dence

interval?

27. Assume that you have available a Uniform[0,1] random number generator.

Give a precise algorithm for generating observations from a distribution

√

which has probability density function f (x) = x3 , 0 < x < 2. Record the

time necessary to generate the sample mean of 100,000 random variables

PROBLEMS 203

with this distribution. Determine the standard error of the sample mean.

How large would the simulation need to be if we wanted to estimate the

mean within 0.01 with a 95% con¬dence interval?

28. Assume that you have available a Uniform[0,1] random number genera-

tor. Give a precise algorithm for generating observations from a discrete

distribution with probability function

x= 0 1 2 3 4 5

P[X=x]= 0.1 0.2 0.25 0.3 0.1 0.05

Record the time necessary to generate the sample mean of 100,000 random

variables with this distribution. Compare the sample mean and variance

with their theoretical values. How large would the simulation need to be

if we wanted to estimate the mean within 0.01 with a 95% con¬dence

interval? How much di¬erence does it make to the time if we use an

optimal binary search? If we use the alias method?

29. Give an algorithm for generating observations from a distribution which

x+x3 +x5

has cumulative distribution function F (x) = < x < 1. Record

,0

3

the time necessary to generate the sample mean of 100,000 random vari-

ables with this distribution. (Hint: Suppose we generate X1 with cumu-

lative distribution function F1 (x) and X2 with cumulative distribution

function F2 (x) , X3 with cumulative distribution function F3 (x) We then

generate J = 1, 2, or 3 such that P[J = j] = pj and output the value

XJ . What is the cumulative distribution function of the random variable

output?)

30. Consider independent random variables Xi i = 1, 2, 3 with cumulative

distribution function

204 CHAPTER 3. BASIC MONTE CARLO METHODS

§

⎪ x2 ,

⎪ i=1

⎪

⎨

ex ’1

Fi (x) = i=2

⎪ e’1

⎪

⎪

©

xex’1 , i=3

for 0 < x < 1. Explain how to obtain random variables with cumulative

distribution function G(x) = Π3 Fi (x) and G(X) = 1’Π3 (1’Fi (x)).

i=1 i=1

(Hint: consider the cumulative distribution function of the minimum and

maximum).

31. Suppose we wish to estimate a random variable X having cumulative

distribution function F (x) using the inverse transform theorem, but the

exact cumulative distribution function is not available. We do, however,

b b

have an unbiased estimator F (x) of F (x) so that 0 · F (x) · 1 and E

b

F (x) = F (x) for all x. Show that provided the uniform variate U is

b b

independent of F (x), the random variable X = F ’1 (U ) has cumulative

distribution function F (x).

32. Give an algorithm for generating a random variable with probability den-

sity function

f (x) = 30(x2 ’ 2x3 + x4 ), 0<x<1

Discuss the e¬ciency of your approach.

33. The interarrival times between consecutive buses at a certain bus stop are

independent uniform[0, 1] random variables starting at clock time t = 0.

You arrive at the bus stop at time t = 1. Determine by simulation the

expected time that you will have to wait for the next bus. Is it more than

1/2 ? Explain.

√

34. What is the probability density function of X = a(1 ’ U ) where U ∼

U [0, 1]?

PROBLEMS 205

35. Develop an algorithm for generating variates from the density:

√ 2 2 2

f (x) = 2/ πe2a’x ’a /x , x > 0

36. Develop an algorithm for generating variates from the density:

2

, for ’ ∞ < x < ∞

f (x) =

eπx + e’πx

37. Explain how the following algorithm works and what distribution is gen-

erated.

(a) Let I = 0

(b) Generate U ∼ U [0, 1]and set T = U .

(c) Generate U — . IF U · U — return X = I + T .

(d) Generate U . If U · U — go to c.

(e) I = I + 1. Go to b

38. Obtain generators for the following distributions:

(a) Rayleigh

x ’x2 /2σ2

,x ≥ 0 (3.44)

f (x) = e

σ2

(b) Triangular

2 x

(1 ’ ), 0 · x · a (3.45)

f (x) =

a a

√

X2 + Y 2

39. Show that if (X, Y ) are independent standard normal variates, then

has the distribution of the square root of a chi-squared(2) (i.e. exponen-

tial(2)) variable and arctan(Y /X) is uniform on [0, 2π].

40. Generate the pair of random variables (X, Y )

(3.46)

(X, Y ) = R(cos˜, sin˜)

where we use a random number generator with poor lattice properties

such as the generator xn+1 = (383xn + 263) mod(10,000) to generate

206 CHAPTER 3. BASIC MONTE CARLO METHODS

our uniform random numbers. Use this generator together with the Box-

Mueller algorithm to generate 5000 pairs of independent random normal

numbers. Plot the results. Do they appear independent?

41. Assume that a option has payo¬ at expiry one year from now (T = 1)

ST ’20

given by the function g(ST ) = 0, ST < 20, and g(ST ) = ST , ST > 20.

What is the approximate present value of the option assuming that the

risk-neutral interest rate is 5 percent, the current price of the stock is 20,

and the annual volatility is 20 percent. Determine this by simulating 1000

stock prices ST and averaging the discounted return from a corresponding

option. Repeat with 100000 simulations. What can you say about the

precision?

42. (Log-normal generator ) Describe an algorithm for generating log-normal

random variables with probability density function given by

1

√ exp{’(logx ’ log· + σ 2 /2)2 /2σ 2 }. (3.47)

g(x|·, σ) =

xσ 2π

43. (hedging with futures). I need to buy 1000 barrels of heating oil on No-

vember 1 1998. On June 1, I go long a December futures contract which

allows me to purchase 1000 barrels of heating oil on December 1 for $20

per barrel. Suppose we have observed that the price of heating oil is log-

normally distributed with monthly volatility 2 percent. The spot interest

rate is presently 5 percent per annum

(a) What is the value of the oil future on November 1 as a function of

the current price of oil?

(b) Determine by simulation what is the standard deviation of the value

of my portfolio on November 1 assuming I sell the futures contract

at that time.

(c) How much di¬erence would it have made if I had purchased the op-

timal number of futures rather than 1000?

PROBLEMS 207

44. (Multivariate Normal generator ) Suppose we want to generate a mul-

tivariate normal random vector (X1 , X2, ..., XN ) having mean vector

(µ1 , ..., µN ) and covariance matrix the N — N matrix Σ. The usual pro-

cedure involves a decomposition of Σ into factors such that A0 A = Σ. For

example, A could be determined from the Cholesky decomposition, in Mat-

lab, A=chol(sigma), or in R, A= chol(sigma, pivot = FALSE, LINPACK

= pivot) which provides such a matrix A which is also upper triangular,

in the case that Σ is positive de¬nite. Show that if Z = (Z1 , ..., ZN ) is a

vector of independent standard normal random variables then the vector

X = (µ1 , ..., µN ) + ZA has the desired distribution.

45. (Ahrens-Dieter) Show that the rejection algorithm of Ahrens and Dieter

(b = 1) has rejection constant c that is bounded for all ±²(0, 1] and ap-

proaches 1 as ± ’ 0.

46. What distribution is generated by the following algorithm where U is

p p

uniform[0, 1] and V is uniform [’ 2/e, 2/e]?

(a) GENERATE U, V

(b) PUT X = V /U

(c) IF ’ln(U ) < X 2 /4,GO TO a.; ELSE RETURN X.

(Hint: First prove the following result due to Kinderman and Mona-

han, 1977: Suppose (U, V ) are uniformly distributed on the region

pv

{(u, v); 0 · v · 1, 0 · u · f ( u )} for some integrable function

1 ≥ f (x) ≥ 0. Then V /U has probability density function cf (x)

R

with c = 1/ f (x)dx. )

47. (Euler vs. Milstein Approximation) Use the Milstein approximation with

step size .001 to simulate a geometric Brownian motion of the form

dSt = .07St dt + .2St dWt

208 CHAPTER 3. BASIC MONTE CARLO METHODS

Compare both the Euler and the Milstein approximations using di¬erent

step sizes, say ∆t = 0.01, 0.02, 0.05, 0.1 and use each approximation to

price an at-the-money call option assuming S0 = 50 and expiry at T =

0.5. How do the two methods compare both for accurately pricing the call

option and for the amount of computing time required?

48. (Cox, Ingersoll, Ross model for interest rates) Use the Milstein approxi-

mation to simulate paths from a CIR model of the form

√

drt = k(b ’ rt ) + σ rt dWt

and plot a histogram of the distribution of r1 assuming that r0 = .05 for

b = 0.04. What are the e¬ects of the parameters k and b?

49. Simulate independent random variables from the Normal Inverse Gamma

distribution using parameter values so that the expected value, the vari-

ance, the skewness and the kurtosis of daily returns are respectively 0,

0.004, 0.5 and 4 respectively. Evaluate an at the money call option with

time to maturity 250 days using these simulated values and compare the

price of the option with the Black-Scholes price. Repeat with an option

whose strike is 10% over the initial stock price and one 10% under.

50. Suppose interest rates follow the constant elasticity of variance process of

the form

drt = k(b ’ rt ) + σ|rt | γ dWt

for parameters value γ, b, k > 0. For various values of the parameters k, γ

and for b = 0.04 use both Euler and Milsten to generate paths from this

process. Draw conclusions about the following:

(a) When does the marginal distribution of rt appear to approach a

steady state solution. Plot the histogram of this steady state dis-

tribution.

PROBLEMS 209

(b) Are there simulations that result in a negative value of r? How do

you rectify this problem?

(c) What does the parameter σ represent? Is it the annual volatility of

the process?

51. Consider a sequence of independent random numbers X1 , X2 , ...with a

continuous distribution and let M be the ¬rst one that is less than its

predecessor:

M = min{n; X1 · X2 · ... · Xn’1 > Xn }

P∞

(a) Use the identity E(M ) = P [M > n} to show E(M ) = e.

n=0

(b) Use 100,000 simulation runs and part a to estimate e with a 95%

con¬dence interval.

(c) How many simulations are required if you wish to estimate e within

0.005 (using a 95% con¬dence interval)?

52. Daily relative losses from a portfolio are assumed to follow the NIG dis-

tribution with parameters

± = 100, δ = .02, β = 0, µ = 0.

In other words if the portfolio is worth St at the end of day t, then

St ’ St+1

St

has the above NIG distribution. Assume daily relative losses are indepen-

dent of one another. Assume S0 = 100000. Use simulation to determine

the weekly 99% VAR, i.e. the value x such that

P [S5 ’ S0 · x] = 0.99

Compare this result with the VAR if we replace the NIG distribution with

the Normal having the same mean and variance.

210 CHAPTER 3. BASIC MONTE CARLO METHODS

53. The interarrival times between consecutive buses at a certain bus stop

are independent Beta(1/2,1/2) random variables with probability density

function

1

f (x) = p , 0 < x < 1.

π x(1 ’ x)

starting at clock time t = 0. You arrive at the bus stop at time t = 10.

Determine by simulation the expected time that you will have to wait for

the next bus. Is it more than 1/2 ? Repeat when the interarrival times

have the distribution of Y 2 where Y is an exponential with expected value

1

Compare your average wait times with the expected time between buses

2.

and explain how your results are possible.