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

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!

Answer 1

[John D'Errico]

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.

Answer 2

[Tommaso Belluzzo]

Absolutely amazing, and thank you very much for your great explanation!