<<

. 4
( 4)






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.

<<

. 4
( 4)