Finding Multiple Roots In Matlab (Equivalent of R rootSolve/multiroot)
11 views (last 30 days)
Show older comments
Tommaso Belluzzo
on 7 Mar 2020
Answered: Tommaso Belluzzo
on 8 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 Comments
Accepted Answer
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.
0 Comments
More Answers (1)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!