4404
Simulation
Prof. IEOR Jose Blanchet
Solutions to Chapter 8 Ross Book
May 1, 2008
Page 1 of 4
Solutions to Chapter 8 Ross Book
8.3
(a) Dene the indicate function 1( 5=1 iXi 21.6), then E [1( 5=1 iXi 21.6)] =
i
i
P ( 5=1 iXi 21.6). So, 1( 5=1 iXi 21.6) is an estimator. Because logU gives a
i
i
an exponential random variable with mean 1, denote h(U) = 1( 5=1 ilogUi 21.6),
i
where Ui is Uniform r.v for all i = 1, ..., 5.
(b) (b) Let Y1 = h(U) and Y2 = h(1 U), so Y1 +Y2 is an antithetic estimator. Because
2
h(U) is decreasing in each Ui , Y1 and Y2 are negatively correlated.
(c) Yes, it is ecient, because the new estimator has smaller variance.
8.5
(a) Consider the control variate Z , so the estimator is Z 3 +eZ +c Z , where c = Cov (Z 3 +
1
eZ , Z ) and can be approximated by c = n1 n (Xi X )(Zi Z ), where Xi =
i=1
3 + eZi . Let X + c Y be the controlled estimator with variance 2 = V ar (X + c Y ) =
Zi
1
n (V ar (X ) Cov (X, Z ))
(b) The 95% condence interval is (X + c Y 1.96
length is within 0.1 require that
2 1.96
2
M,X
+ c Y + 1.96
2
M ).
That the
2
0.1
M
1 96
i.e., M ( 20..1 )2 2 . 2 can be approximated to 11.7611 by running 104 iterations.
Therefore, we need to at least M = 2.0 104 replications. By running simulation, we
get the 95% condence interval (1.5976, 1.6889).
The following is the Matlab code:
clear all;
N=2*10^4;
X=zeros(N,1);
Z=zeros(N,1);
for i=1:N
Z(i)=randn;
X(i)=Z(i)^3+exp(Z(i));
end
covXZ=cov(X,Z);
cov_XZ=covXZ(1,2);% convariance(X,Z)
c_star=-cov_XZ/var(Z); % c^*
theta_bar=mean(X)+c_star*mean(Z) % sample mean
sigma_square=(var(X)-cov_XZ^2)
Confidence_interval=[theta_bar-1.96*sqrt(sigma_square/N),
theta_bar+1.96*sqrt(sigma_square/N)]
8.10
(a) We compute c as follows:
c =
cov (I, X )
var(X )
2
IEOR 4404, Solutions to Chapter 8 Ross Book
Now, cov (I, X ) = E [IX ] E [I ]E [X ]. And, E [IX ] = E [IX |X a]P [X a] =
1
E [X |X a]P [X a] = a2 /2, while E [I ]E [X ] = a/2. Therefore, cov (I, X ) = 2 a(a
1
1) 0. Recall V ar(X ) = 12 and V ar(I ) = a a2 .
So, the percentage of variance reduction is
[ 1 a(a 1)]2
[Cov (I, X )]2
= 21
= 3a(1 a)
V ar(I )V ar(X )
12 a(1 a)
(b) If X were exponential with mean 1, E [IX ] = E [X |X a]P [X a] = 1 (1 + a)ea ,
E [X ] = 1, V ar(X ) = 1. So, Cov (I, X ) = (1 a) (1 + a)ea .
So, the percentage of variance reduction is
[(1 a) (1 + a)ea ]2
[Cov (I, X )]2
=
V ar(I )V ar(X )
a(1 a)
(c) Here, note that as X increases, I has a tendency to decrease which justies the negative
correlation.
8.17
(a) The raw estimator would be to generate and X Y and average the number of times we
get that X + Y t.
(b) Using indicator function, P (X + Y t) = E [1(X + Y t)]. We can condition on
either X or Y .
(c) Once again, you can use either X or Y here as control variables. For example, E [1(X +
Y t)|X ] + c(Y y ) is new estimator.
8.28 & 8.29
When using importance sampling, we rst need to nd the appropriate tilting that would
give E [X ] 5. To do this recall that, under the tilted distribution, component i has a
Bernoulli distribution with parameter
p,i =
pi
e
pi e
+ 1 pi
where pi = 0.5 + i/50. Therefore, under the tilted distribution, X has mean
20
f () = E [X ] =
20
p,i =
i=1
i=1
pi e
pi e + 1 pi
After plotting the function f () it can be seen that f (2.165) = 4.9990, so we choose
= 2.165. The moment generating function of X is the product of the moment generating
functions of the components, that is,
MX () = 20 (pi e + 1 pi )
i=1
The importance sampling estimator for = P (X 5) is then
Z = I (X 5)e
X
MX ( )
where the components dening X are generated as Bernoulli random variables with parameters p,i , i = 1, 2, ..., 20. For 106 simulations the following values were obtained
Z (n) = 6.7621e 10
S 2 (n) = 3.7427e 14
where Z (n) is the sample mean, and S 2 (n) is the sample variance.
Similarly, dene
Z 1 = I (X == 5)e X MX ( )
IEOR 4404, Solutions to Chapter 8 Ross Book
. By running simulation, we get sample mean and sample variance
Z (n) = 6.6980e 10
S 2 (n) = 3.7386e 14
P (X 5) = 6.7621e 10 0
P (X = 5)
P (X = 5|X 5) =
0.9905
P (X 5)
% function f()
function y=f(x)
p=zeros(20,1);
pp=zeros(20,1);
for i=1:20
p(i)=.5+i/50;
pp(i)=p(i)*exp(x)/(p(i)*exp(x)+1-p(i));
end y=sum(pp);
Then, by plotting to nd
% prob 3
clear all;
for i=1:10/0.001
x(i)=-5+.001*i;
y(i)=f(x(i));
if i>1 && (y(i)-5)*(y(i-1)-5)<=0
y(i-1)
x(i-1)
plot(x,y)
break;
end
end
The code for the importance sampling
% inportance sampling
clear all;
theta=-2.165;
p=zeros(20,1);
p_M=zeros(20,1);
for i=1:20
p(i)=.5+i/50;
p_M(i)=p(i)*exp(theta)+1-p(i);
end
x_Moment=prod(p_M);
n=10^6; % number of simulations
z=zeros(n,1);
z1=zeros(n,1);
x=zeros(20,1);
for j=1:n
3
4
IEOR 4404, Solutions to Chapter 8 Ross Book
for k=1:20
if rand<=p(k)
x(k)=1;
else
x(k)=0;
end
end
if sum(x)<=5
z(j)=exp(-theta*sum(x))*x_Moment;
else
z(j)=0;
end
if sum(x)==5
z1(j)=exp(-theta*sum(x))*x_Moment;
else
z1(j)=0;
end
end
Sample_mean=mean(z) % p(X<=5)
Sample_variance=var(z)
Sample_mean1=mean(z1) %p(x==5)
Sample_variance1=var(z1)
Conditional_prob=Sample_mean1/Sample_mean