integration - MATLAB: Pricing a digital option, Monte Carlo vs. explicit integral formula?

Question

I am stuck with the following problem using MATLAB:

Let Z be lognormally distributed such that ln Z has mean m and variance w. Let eta be a negative number and c a positive constant.

I am trying to compute the expected value (let I(Z<=c) denote the indicator function of the set (Z<=c))

E[Z^(eta+1) I(Z<=c)] = (1/sqrt(w)) integral_0^c x^(eta) phi((ln x - m)/sqrt(w)) dx,

where phi() denotes the probability distribution function of a standard normal random variable.

First thing I did was to simulate 10.000 trials of Z, set the entries of the vector with value >c to 0, raised to the power of (eta+1) and then calculated the mean. This should give me the MC estimate of the expected value.

ST = random('Lognormal', m,w_sq,10000,1);
hlp = zeros(10000,1);
hlp(ST<=2) = ST(ST<=2);
hlp(hlp>0) = hlp(hlp>0).^(eta1+1); % 0^(eta1+1) gives infinity
mean(hlp)

For the integral I used the following code

tmpp = integral(@(x) x.^(eta1) .* normpdf((log(x)-m)/sqrt(w_sq),0,c);
tmpp / sqrt(w_sq(1))

Unfortunately the procedures lead to totally different results, although mathematically they should be they same.

This whole thing is part of a bigger code and it would be a lot more convenient for me to use the integral version. Originally I tried to double check using the MC simulation and then saw that something must be wrong...

Can someone help?

Answer 1

For the first piece of code,

I guess the reason why you have unexpected result is when you are doing the calculations on hlp, you try to avoid 0 values, as 0^eta will blow up - and that is not wanted result, so you simply drop it. But in the last step, mean(hlp) will take all values in array hlp, including those 0's. Try this:

mean(hlp(hlp>0)).

My results fall in roughly 2.x with 10,000-point simulation, 2.3x ~ 2.4x with 1,000,000 points.

I was wrong. Your problem is you are using too few points. try 10,000,000 points, and you will be satisfied :)

For the second, I have a problem understanding the way you define your variables. (I don't have enough reputations to add a comment, so I put it here. ) Does the "sq" in w_sq mean the square root of w? Because according to the documentation of random, the argument should be "sigma", which is standard deviation. It's natural to define SD as the square root of the variance, which, I guess, is w. Then you are doing right in the first piece.

If so, in the second piece of your code, why do you put sqrt() on w_sq? Do you mean taking a 4th root of the variance? From your definition of expectation, I don't think it's correct. Please take a look at it.

On the other hand, is w_sq a single number, or an array?

• If it is a number -

tmpp / sqrt(w_sq(1)) should be tmpp / sqrt(w_sq), although they don't make difference in fact.

• If it's an array -

You may want to put all codes in a for loop. The loop iterates over w_sq, each time it picks out one of the elements in the array (name the variable as say w_sq_elem), and let the rest of the code do things as if it is a single number.

Anyway, (log(x)-m)/sqrt(w_sq) and tmpp / sqrt(w_sq(1)) are telling different information on w_sq. The first one assumes it's a number, so the division can be simply a /, rather than ./. The second one indicates it's an array so you are picking its first element. But an array does not make sense here because, in my understanding, you are not dividing 10,000 points in x by 10,000 different variances.

m = 3;
w_sq = 2;
eta1 = -2;
ST = random('Lognormal',m,w_sq,10000000,1);
hlp = zeros(10000000,1);
hlp(ST<=2) = ST(ST<=2);
hlp(hlp>0) = hlp(hlp>0).^(eta1+1); % 0^(eta1+1) gives infinity
mean(hlp)

tmpp = integral(@(x) x.^(eta1) .* normpdf((log(x)-m)/w_sq),0,2) ;
tmpp / w_sq

>> untitled

ans =

    0.2944


ans =

    0.2948

>>