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

8 vues (au cours des 30 derniers jours)
Tommaso Belluzzo le 7 Mar 2020
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 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

John D'Errico le 8 Mar 2020
Modifié(e) : John D'Errico le 8 Mar 2020
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 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Plus de réponses (1)

Tommaso Belluzzo le 8 Mar 2020
Absolutely amazing, and thank you very much for your great explanation!
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Calculus dans Help Center et File Exchange

R2018a

### Community Treasure Hunt

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

Start Hunting!

Translated by