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

22 views (last 30 days)
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!

Accepted Answer

John D'Errico
John D'Errico on 8 Mar 2020
Edited: John D'Errico on 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.

More Answers (1)

Tommaso Belluzzo
Tommaso Belluzzo on 8 Mar 2020
Absolutely amazing, and thank you very much for your great explanation!

Tags

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!