We denote by \(M'\) the transpose of a matrix \(M\) and by writing \(Y\sim N(0,s)\) we mean that \(Y\) is normally distributed with mean 0 ans standard deviation \(s\).
Often the random variable we are interested in depends on more than one normally distributed input variables so that it’s value is strongly dependend on a liear combination of the variables. Consider, for example the random variable \[Y=\frac{e^{0.2(X_1-3X_2)}}{1+e^{-X_1^2}+e^{-X_2^2}},\] where \(X_1\) and \(X_2\) are independent random variables from the distribution \(N(0,2)\). We can not say much about it’s value if we fix only \(X_1\) or \(X_2\), but if we know the value of the cobination \(X_1-3X_2\), we know quite a lot about possible values of \(Y\). This means that if we want to use stratified sampling for computing the expected value of \(Y\), we should try to define events for stratums by restricting the values of this combination to certain intervals of real numbers. For example, we may want to consider stratums corresponding to events of the form \(A_i =\{a_{i-1} \leq X_1-3X_2 <a_{i}\}\), where \(a_0,a_1,\ldots\) are some increasing numbers, but then we have to find a way to generate \(X_1, X_2\) values corresponding to a given stratum.
We know that the random variable \(W=X_1-3X_2\) is also normally distributed.
Exercise Find mean and standard deviation of \(W\). You can check your answer by expanding the link below. Click to check your answer!
We know that it is easy to generate values of \(W\) for given stratum - the general formula for generating the values of W under the condition that \(a\leq W < b\) is \[W=F_W^{-1}(U_{F_W(a),F_W(b)}),\] where \(F_W\) is the cumulative distribution function of \(W\), \(F_W^{-1}\) is it’s inverse (so called quantile function of the distribution) and \(U_{F_W(a),F_W(b)}\) is a random variable with uniform distribution in the interval \((F_W(a),F_W(b))\). We define the conditions by quantiles of the distribution, so \(F_W(a_i)=\frac{i}{k}\), where \(k\) is the number of stratums and \(i\) is the number of a stratum of interest.
Exercise Generate 10 values of \(W\) from the second stratum out of 10 stratums by using the command set.seed(2020) immediatly before it. You can check your answer by expanding the link below. Hint: in the function qnorm you can indicate the parameters of the distribution by mean=...,sd=....
[1] -6.168733 -6.848970 -6.241591 -6.617990 -7.637233 -7.868137 -7.660107
[8] -6.852128 -8.095945 -6.237196
Unfortunately it is not enough to generate \(W\) values, we also need \(X_1\) and \(X_2\) so that the sum is \(W\) ad which are independent with correct distribution. But fortunately we can get \(X_1\) and \(X_2\) with the right properties by generating first independent \(W,Z_1\) and \(Z_2\) (so that \(Z_1\) and \(Z_2\) are from \(N(0,2)\)) and then combining them so that we get two independent vectors from \(N(0,2)\) so that the sum is \(W\). The formulas for this case look like this: \[X_1=\frac{1}{10}W+Z_1-\frac{1}{10}(Z_1-3Z_2),\] \[X_2=-\frac{3}{10}W+Z_2+\frac{3}{10}(Z_1-3Z_2).\] The random variables \(X_1\) and \(X_2\) are normally distributed since they are linear combinations of independent normally distributed variables, moreover they are jointly normally distributed for the same reason. So if one shows that the covariance of \(X_1\) and \(X_2\) is 0, then for jointly normal random variables it means that they are independent. If other properties also hold, we have a way to generate \(X_1\) and \(X_2\) so that \(W=X_1-3X_2\) is in a given stratum - we just have to start with values of W from a given stratu and define \(X_1\) and \(X_2\) according to the formulas above.
Exercise Check that \(X_1-3X_2=W\), that covariance of \(X_1\) and \(X_2\) is 0 and that variance of both \(X_1\) and \(X_2\) is 4.
Let us now see how the the original MC method works and how much we can gain in speed if we use stratified sampling with respect to \(X_1\) or with respect to \(W=X_1-3X_2\).
Exercise Compute the expected value of \[Y=\frac{e^{0.2(X_1-3X_2)}}{1+e^{-X_1^2}+e^{-X_2^2}},\] where \(X_1\) and \(X_2\) are independent random variables from the distribution \(N(0,2)\) with error less than 0.01 with probability 0.90. Click to check your results!
EY is approximately 1.63183736917303 , the number of generated values is: 416000
Click to see the code!
g <- function(X){
return(exp(0.2*(X[,1]-3*X[,2]))/(1+exp(-X[,1]^2)+exp(-X[,2]^2)))
}
gen <- function(n){
return(cbind(rnorm(n,0,2),rnorm(n,0,2)))
}
ans_mc <- MC2(g,gen,0.01,0.1)
paste("EY is approximately",ans_mc[1],", the number of generated values is:",ans_mc[2])Exercise Repeat previous computations in the case using stratified sampling with respect to \(X_1\) variable.
Click to check your results!
EY is approximately 1.64990621413499 , the number of generated values is: 408000
Click to see the code!
gen <- function(n){
return(cbind(qnorm(runif(n,(i-1)/k,i/k),0,2),rnorm(n,0,2)))
}
epsilon<-0.01
alpha<-0.1
k<-10
result<-0
n<-0
for(i in 1:k){
p_i<-1/k
answer<-MC2(g, gen,alpha=alpha,error=epsilon/sqrt(p_i))
result<-result+p_i*answer[1]
n<-n+answer[2]
}
paste("EY is approximately",result,", the number of generated values is:",n)As you see, stratification with respect to \(X_1\) did not give much effect on the number of values generated.
Exercise Repeat previous computations in the case using stratified sampling with respect to \(W=X_1-3X_2\) variable.
Click to check your results!
EY is approximately 1.63861024613585 , the number of generated values is: 242000
Click to see the code!
gen <- function(n){
W <- qnorm(runif(n,(i-1)/k,i/k),0,sqrt(40))
Z1 <-rnorm(n,0,2)
Z2 <- rnorm(n,0,2)
X1 <- 0.1*W+Z1-0.1*(Z1-3*Z2)
X2 <- -0.3*W+Z2+0.3*(Z1-3*Z2)
return(cbind(X1,X2))
}
epsilon<-0.01
alpha<-0.1
k<-10
result<-0
n<-0
for(i in 1:k){
p_i<-1/k
answer<-MC2(g, gen,alpha=alpha,error=epsilon/sqrt(p_i))
result<-result+p_i*answer[1]
n<-n+answer[2]
}
paste("EY is approximately",result,", the number of generated values is:",n)Stratification with respect to \(W\) is better, but the speedup is still not very large. Using larger value of \(k\) and optimal stratified sampling would give better results.
It turns out that the formulas for generating independent normally distributed random variables are not very complicated if they are presented in the language of linear algebra. When \({\bf v}=(v_1,\ldots,v_m)'\) is a non-zero vector (i.e. \(\|{\bf v}\|=\sqrt{v_1^2+\ldots+v_m^2}>0\)), \(W\sim N(0,a\|{\bf v}\|)\) and \(Z_i\sim N(0,a),\ i=1,\ldots,m\) are independent random variables and \(a>0\), then it is easy to check that by defining a vector \(X\) of random variables \(X_i,\ i=1,\dots,m\) as \[{\bf X}=\frac{W}{\|{\bf v}\|^2}{\bf v}+{\bf Z}-\frac{({\bf v}'{\bf Z})}{\|{\bf v}\|^2}{\bf v}\] we have that the components of \(X\) are independent and have distribution \(N(0,a)\), and also \({\bf v}'{\bf X}=W\).
If we want to generate at once more than one vector, each corresponding to different value of \(W\), then the former formula can be written as \[{\bf X}=\frac{1}{\|{\bf v}\|^2}{\bf v}{\bf W}+{\bf Z}-\frac{({\bf v}{\bf v}'{\bf Z})}{\|{\bf v}\|^2}, \] where \({\bf X}\) is now a \(m\times n\) matrix with independent normally \(N(0,a)\) distributed random variables, \({\bf W}\) is a \(1\times n\) matrix (a row vector) of independent \(N(0,a\|{\bf v}\|)\) distributed random variables and \({\bf Z}\) is a \(m\times n\) matrix with independent normally \(N(0,a)\) distributed random variables. The matrix \({\bf X}\) has now the property \({\bf v}'{\bf X}={\bf W}\) (i.e. each column sums with weights \(v_i\) to the value of \(W_i\)). Thus we can generate independent normally distributed random variables so that we first generate the value of a linear combination of the variables and then determine the variables itself.
Exercise. Write a procedure that, given the inputs \(k,m,T\), would draw \(k\) different Brownian motion paths such that every stratum (based on the value of \(B(T)\) and defined as in the previous lab) would include the terminal value of exactly one path. Hints: the increments of Brownian motion over intervals of length \(\Delta t\) are independent and normally distributed with mean 0 and standard deviation \(\sqrt{\Delta t}\). The final value of a trajectory of Brownian motion is the sum of increments, so the vector \(\bf v\) is vector of ones. The matrix multiplication is %*% and traspose of a matrix M can be computed as t(M).
Solution
Click to see the code!
strat_BM <- function(k,m,T){
dt<-T/m #step size
a<-sqrt(dt)
#we want to stratify according to the sum
v<-matrix(1,nrow=m,ncol=1)
v_norm<-sqrt(sum(v^2))
#generate values for W, one from each stratum
i<-1:k
#one value from every stratum
W<-matrix(a*v_norm*qnorm(runif(k,(i-1)/k,i/k)),nrow=1)
#generate increments of Brwonian motion
#so that sum is W
Z<-matrix(rnorm(m*k,sd=a),nrow=m,ncol=k)
dB<-1/v_norm^2*v%*%W+Z-1/v_norm^2*v%*%t(v)%*%Z
#compute the trajectories of Brownian motion
B<-matrix(0,nrow=m+1,ncol=k)
for(i in 1:m){
B[i+1,]<-B[i,]+dB[i,]
}
t<-seq(0,T,length.out=m+1)
matplot(t,B,type="l")
}
strat_BM(10,100,0.5)Exercise
Enhance the stock price generation function that is based on Euler’s method and non-constant volatility so that one could use the optimal stratified sampling. Using optimal stratified sampling find the price of an European call with total error \(0.01\) in the case when \(r=0.06\), \(D=0.03\), \(\sigma(s)=\frac{95}{95+s}\), \(T=0.5\), \(S(0)=105\), \(E=100\), using \(\alpha=0.05\).
Solution. First, define a function for generating a matrix of increments of Brownian motion so that the final values of corresponding Brownian motion trajectory are in a given stratum.
Click to see the code!
dB_strat<-function(n,i,k,T,v1){
m<-length(v1)
a<-sqrt(T/m)
#make v to be a column vector
v<-matrix(v1,ncol=1)
v_norm<-sqrt(sum(v^2))
#generate values for W, one from each trajectory
W<-matrix(a*v_norm*qnorm(runif(n,(i-1)/k,i/k)),nrow=1)
#generate increments of Brwonian motion
#so that sum is W
Z<-matrix(rnorm(m*n,sd=a),nrow=m,ncol=n)
dB<-1/v_norm^2*v%*%W+Z-1/v_norm^2*v%*%t(v)%*%Z
return(dB)
}Next modify the generator S_euler2 from Lab 4 to generate values of stock prices corresponding to Brownian motion trajectories with the final value in a given stratum
Click to see the code!
S_euler2_strat<-function(n,S0,m,T,r,D,sigma,i,k){
#generator for variable volatility sigma(t,s)
#output: n approximate values for S(T)
#computed with Euler's method with m time steps
#uses stratified sampling according to values of B(T)
delta_t <-T/m
S<-matrix(NA,nrow=m+1,ncol=n)
S[1,]<-S0 #set the starting value for each trajectory
t<-seq(0,T,length.out=m+1)
#######
#Here generate the increments of Brownian motion for a given stratum
dB<-dB_strat(n,i,k,T,rep(1,m))
#since i has a different meaning, use j as a variable for the for cycle
#it is not a good habit to assign values to input parameters
for(j in 1:m ){
Xi<-dB[j,] #each row corresponds to a different time moment
S[j+1,]<-S[j,]*(1+(r-D)*delta_t+sigma(t=t[j],s=S[j,])*Xi)
}
return(S[m+1,])
}Next, include the functions used for optimal stratified sampling
Define the functions and parameters needed for solving the problem
Click to see the code!
r<-0.06
D<-0.03
sigma <- function(t,s){
return(95/(95+s))
}
T<-0.5
S0<-105
E<-100
alpha<-0.05
#for stratified sampling generator should be a function of n,i,k
gen_strat<-function(n,i,k){
return(S_euler2_strat(n,S0,m,T,r,D,sigma,i,k))
}
g_call<-function(s){
return(exp(-r*T)*pmax(s-E,0))
}Now use the usual procedure for computing the answer with a given accuracy. Let us use \(k=50\) stratums.
Click to see the final answer
Price is 16.953, number of generations in the last computation: 151900, MC error estimate of the last computation: 0.0046
Click to see the code!
total_error<-0.01
k<-50
m0<-5
eps0 <- 0.01 #come back, if works fast
# compute the approximate option price for m=m0
m<- m0
answer1 <- Stratified_error(g_call,gen_strat,k,alpha,eps0)
Vm0 <- answer1[1]
#second computation
m <- 2*m0
answer2 <- Stratified_error(g_call,gen_strat,k,alpha,eps0)
V2m0 <- answer2[1]
# recommendataion - check if the difference is larger than 2*eps0
# if not, and if the estimate of m we obtain is large,
#we may want to go back and reduce eps0 to get better estimate for m
# q is 1 for Euler's method
q <- 1
#estimate constant C
#since with the optimal stratified sampling the errors of first 2 computations
#are not always less than the given error eps0, let us use
#the sum of the actual error estimates instead of eps0 values
Cbar <- (2*m0)^q/(2^q-1)*(abs(Vm0-V2m0)+answer1[3]+answer2[3])
#find m such that Cbar/m**q<=total_error/2
m <- ceiling((2*Cbar/total_error)**(1/q) )
m
#final computation
MCerror <- total_error/2
answer_final <- Stratified_error(g_call,gen_strat,k,alpha,total_error/2)
answer_final