MATLAB Answers

Finding Multiple Roots In Matlab (Equivalent of R rootSolve/multiroot)

38 ビュー (過去 30 日間)
Tommaso Belluzzo
Tommaso Belluzzo 2020 年 3 月 7 日
回答済み: Tommaso Belluzzo 2020 年 3 月 8 日
Hi all! I'm trying to replicate in Matlab the following R code:
library(rootSolve)
A <- 560223.07
B <- 358176
r <- 0.0171
Ve <- 0.1276263
T <- 1
D <- exp(-r * T)
F <- function(x)
{
d1 <- x[1]
d2 <- x[2]
V <- (1 - ((B * D * pnorm(d2)) / (A * pnorm(d1)))) * Ve
F1 <- ((log(A / B) + ((r + (0.5 * V^2)) * T)) / (V * sqrt(T))) - d1
F2 <- d1 - (V * sqrt(T)) - d2
c(F1=F1,F2=F2)
}
solution <- multiroot(f = F, start = c(0,0))
solution
Which produces the following output:
$root
[1] 9.818820 9.771408
$f.root
F1 F2
-7.854641e-10 -3.478249e-10
$iter
[1] 3
$estim.precis
[1] 5.666445e-10
Now, let's move to my Matlab implementation:
A = 560223.07;
B = 358176;
r = 0.0171;
Ve = 0.1276263;
T = 1;
D = exp(-r * T);
solution = fsolve(@(x)objective(x,A,B,T,Ve,r,D),[0 0]);
solution
function F = objective(x,A,B,T,Ve,r,D)
d1 = x(1);
d2 = x(2);
V = (1 - ((B * D * normcdf(d2)) / (A * normcdf(d1)))) * Ve;
F1 = ((log(A / B) + ((r + (0.5 * V^2)) * T)) / (V * sqrt(T))) - d1;
F2 = d1 - (V * sqrt(T)) - d2;
F = [F1 F2];
end
The solution is very different from the one proposed by R (which I know to be the correct one):
solution =
1.9522 -0.6945
I also tried the trick of minimizing the square of the joint function outputs as follows:
solution = fminunc(@(x)sum(objective(x,A,B,T,Ve,r,D))^2,[0 0]);
solution
But again, the result is not what I'm expecting:
solution =
5.2336 -0.70497
If I run the objective function with the roots proposed by R, there is what I obtain:
objective([9.818820 9.771408],A,B,T,Ve,r,D)
ans =
-2.1915e-07 -4.7316e-07
Maybe I'm using the wrong Matlab tools to solve this problem? Maybe I have to target a specific algorithm through options? Maybe I need to specify the Matlab objective function in a different way?
Any help will be really appreciated!

  0 件のコメント

サインインしてコメントする。

採用された回答

John D'Errico
John D'Errico 2020 年 3 月 8 日
編集済み: John D'Errico 2020 年 3 月 8 日
First, when you do this:
objective([9.818820 9.771408],A,B,T,Ve,r,D)
you are using a 7 digit approximation to the solution. 9.818820 is surely not the exact number produced by the estimation. So you should expect to get an objective that is roughly only accurate to 7 digits or so.
ans =
-2.1915e-07 -4.7316e-07
NEVER just copy numbers from the screen, and then use them in a calculation.
Anyway, how might I try to solve it? I might take a shot with vpasolve.
A = 560223.07;
B = 358176;
r = 0.0171;
Ve = 0.1276263;
T = 1;
D = exp(-r * T);
syms d1 d2 x
ncdf = (erf(x/sqrt(2))+1)/2;
V = (1 - ((B * D * subs(ncdf,d2)) / (A * subs(ncdf,d1)))) * Ve;
F1 = ((log(A / B) + ((r + (0.5 * V^2)) * T)) / (V * sqrt(T))) - d1;
F2 = d1 - (V * sqrt(T)) - d2;
F = [F1 F2];
Sol = vpasolve(F,d1,d2)
Sol =
struct with fields:
d1: [1×1 sym]
d2: [1×1 sym]
Sol.d1
ans =
9.8188197808502313168391797572718
Sol.d2
ans =
9.7714073076927737876772496940379
How well did it do?
vpa(subs(F,Sol))
ans =
[ 8.8162076311671563097655240291668e-39, -1.1479437019748901445007192746311e-40]
Seems ok to me.
Can I get fsolve to do the same? Well, not with the same accuracy. Not even close. Why not? you have some serious numerical problems.
format long g
normcdf(9.8188197808502313168391797572718)
ans =
1
So in double precision, normcdf has problems that far out. DID YOU SERIOUSLY EXPECT BETTER? THINK ABOUT IT! You are using double precision here. 10 sigma? normcdf will produce 1. Can you do better? This is why I switched to syms, because you need more precision than a double can relly offer, at least not unless you reformulate those equations. perhaps to use erfc. That far out, you should recognize that you need more than 22 digits of precision to get something different from 1 in the call to normcdf.
erfc(9.8/sqrt(2))
ans =
1.12585646227532e-22
The point is, this is not a problem with fsolve. This is a problem of using double precision to try to solve that problem using brute force, and hoping it will survive.

  0 件のコメント

サインインしてコメントする。

その他の回答 (1 件)

Tommaso Belluzzo
Tommaso Belluzzo 2020 年 3 月 8 日
Absolutely amazing, and thank you very much for your great explanation!

  0 件のコメント

サインインしてコメントする。

タグ

製品


リリース

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by