MATLAB:为数字期权定价，蒙特卡洛还是显式积分公式? [英] MATLAB: Pricing a digital option, Monte Carlo vs. explicit integral formula?
问题描述
我在使用MATLAB时遇到以下问题:
令Z为对数正态分布，以使ln Z具有均值m和方差w.令eta为负数，c为正常数.我正在尝试计算期望值(让I(Z <= c)表示集合的指标函数(Z <= c))
p [E [Z ^(eta + 1)I(Z <= c)] =(1/sqrt(w))积分_0 ^ cx ^(eta)phi((ln x-m)/sqrt(w)) dx其中phi()表示标准正态随机变量的概率分布函数.
我要做的第一件事是模拟Z的10.000次试验，将值> c的向量的条目设置为0，提高到(eta + 1)的幂，然后计算平均值.这应该给我MC期望值的估计值.
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)
对于积分，我使用了以下代码
tmpp = integral(@(x) x.^(eta1) .* normpdf((log(x)-m)/sqrt(w_sq),0,c);
tmpp / sqrt(w_sq(1))
不幸的是，尽管从数学上讲它们应该是相同的，但这些程序却导致了完全不同的结果.
这整个都是更大代码的一部分，使用整数版本对我来说将更加方便.最初，我尝试使用MC仿真进行仔细检查，然后发现一定有问题...
有人可以帮忙吗?
对于第一段代码，
我想你之所以有意想不到的结果的原因是当你 在
hlp上进行计算时，您尝试避免使用0值，如0^eta会爆炸-这不是想要的结果，因此您只需将其删除即可. 但在最后一步，mean(hlp)将采用数组hlp中的所有值， 包括那些0.试试这个:
mean(hlp(hlp>0)).通过10,000点仿真，我的结果大约是2.x，是2.3x〜 2.4x 1,000,000点.
我错了.您的问题是您使用的积分太少.尝试10,000,000点，您会满意的:)
第二点，我在理解变量定义方式时遇到问题. (我没有足够的名气来添加评论，所以在这里输入了它.)w_sq中的"sq"是否表示w的平方根?因为根据 random 的文档，该参数应为"sigma" ，这是标准偏差.将SD定义为方差的平方根是很自然的，我猜这是w.然后，您在第一部分中就做对了.
如果是这样，为什么在代码的第二部分中，为什么将sqrt()放在w_sq上?您是说要接受方差的第4个根吗?根据您对期望的定义，我认为这是不正确的.请看一看.
另一方面，w_sq是单个数字还是数组?
- 如果是数字-
tmpp / sqrt(w_sq(1))应该为tmpp / sqrt(w_sq)，尽管实际上并没有什么区别.
- 如果是数组-
您可能希望将所有代码放入for循环中.每次选择数组中的一个元素(将变量命名为w_sq_elem)时，循环都会在w_sq上进行迭代，并让其余代码将其当作单个数字来处理.
无论如何，(log(x)-m)/sqrt(w_sq)和tmpp / sqrt(w_sq(1))告诉w_sq的不同信息.第一个假设它是一个数字，因此除法可以简单地是/，而不是./.第二个表示它是一个数组，因此您要选择第一个元素.但是数组在这里没有意义，因为据我了解，您没有将x中的10,000点除以10,000个不同的方差.
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
>>
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?
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, as0^etawill 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 arrayhlp, 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
>>
这篇关于MATLAB:为数字期权定价，蒙特卡洛还是显式积分公式?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
