solve two equations with 2 unknowns(pbar,rho)

Hello everyone!

I need some help to solve a system of two equations to find pbar and rho:

Eq1 = 600000*1000*normcdf((sqrt(rho)*norminv(0.95,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) == 147400000;
Eq2 = 600000*1000*normcdf((sqrt(rho)*norminv(0.99,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) == 293200000;

I tried the following code but it didnt work:

syms pbar rho
  eqns = [600000*1000*normcdf((sqrt(rho)*norminv(0.95,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) == 147400000,600000*1000*normcdf((sqrt(rho)*norminv(0.99,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) == 293200000];
  S = solve(eqns,[pbar,rho])

Thank you in advance!!!

Answer 1

[Stephan]

Try numeric:

format long
S0 = [0.1, 0];
S = fsolve(@eqns,S0)
TestSolution = eqns(S)  % Should be near zero for a good solution
format short

function Eq = eqns(x)
rho = x(1);
pbar = x(2);
Eq(1) = 600000*1000*normcdf((sqrt(rho)*norminv(0.95,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) - 147400000;
Eq(2) = 600000*1000*normcdf((sqrt(rho)*norminv(0.99,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) - 293200000;
end

Answer 2

[John D'Errico]

The symbolic toolbox does not always work well with the stats toolbox. And stats is where normcdf and norminv come from.

So do it with mostly pencil and paper. Start with Eq1.

Eq1 = 600000*1000*normcdf((sqrt(rho)*norminv(0.95,0,1) + norminv(pbar,0,1))/sqrt(1-rho),0,1) == 147400000;

Divide by those constants at the beginning, and take the inverse normcdf of both sides.

(sqrt(rho)*norminv(0.95,0,1) + norminv(pbar,0,1))/sqrt(1-rho) == norminv(147400000/(600000*1000),0,1)

We can reduce some of those numeric terms as...

format long g
norminv(0.95,0,1) 

ans = 
          1.64485362695147

norminv(147400000/(600000*1000))

ans = 
        -0.688189693580052

This leaves Eq1 as:

(sqrt(rho)*1.64485362695147 + norminv(pbar,0,1))/sqrt(1-rho) == -0.688189693580052

Do the same for Eq2.

norminv(0.99,0,1)

ans = 
          2.32634787404084

norminv(293200000/(600000*1000))

ans = 
       -0.0284122759864473

(sqrt(rho)*2.32634787404084 + norminv(pbar,0,1))/sqrt(1-rho) == -0.0284122759864473

We still see pbar in there, inside norminv, in both equations as norminv(pbar,0,1). Just call that number Z. So in the end, we have two equations in two unknowns, totally divorced from the stats toolbox.

syms rho Z
Eq1 = (sqrt(rho)*1.64485362695147 + Z)/sqrt(1-rho) == -0.688189693580052;
Eq2 = (sqrt(rho)*2.32634787404084 + Z)/sqrt(1-rho) == -0.0284122759864473;
parms = solve(Eq1,Eq2,[rho,Z]);
rho = double(parms.rho)

rho = 
         0.483813019861367

pbar = normcdf(double(parms.Z),0,1)

pbar = 
        0.0506542131962006

In fact, I could have taken the pencil and paper solution further, eliminating Z from the problem. Then I could probably have solved for rho by hand too, from the one equation that results. I was too lazy though. But suppose we did that... (as long as rho is not 1)

(EQ1 - EQ2)*sqrt(1-rho)

will look something like this:

12276708948987292*sqrt(rho) == 11885493328088105*sqrt((1 - rho)

Squaring both sides,

rho = (1-rho)*0.937282493509289

We can trivially now solve for rho, as 0.483813019861367. All by pencil and paper and a hand calculator. To no surprise, it agrees with what solve told me before.