Solving equations when variable is included in mu and sigma of distribution function

Hello

I'm trying to solve the following equation

c = 1.60

q = 0.08

equity = 0.08

p = 3.26

sigma = 0.1

mu = 0.08

    syms r

    vpasolve((c*q-equity)*(1+r)/p==(-p/(1+r)+p*(1-normcdf((c*q-equity)*(1+r)/p,(c*mu-equity)*(1+r)/p,c*(1+r)/p*sigma)))*(c*(1-normcdf(q,mu,sigma))-(c*q-equity)*normpdf(q,mu,sigma))/((p/(1+r))^2*normpdf(q,mu,sigma)-c^2*normpdf((c*q-equity)*(1+r)/p,(c*mu-equity)*(1+r)/p,c*(1+r)/p*sigma)),r,[0 1],'random',true)

However, Matlab returns an error message

symengine

Cannot prove '(1809252434946447*r)/36729926713742350 + 1809252434946447/36729926713742350 == 0 &
0.0000000000000000025671305227897154891565148192726*r + 0.0000000000000000025671305227897154891565148192726 <
- (12505220420059*r)/554413988131960 - 12505220420059/554413988131960' literally. To test the statement
mathematically, use isAlways.

   sym/subsindex (line 762)
                  X = find(mupadmex('symobj::logical',A.s,9)) - 1;

   sym/privsubsasgn (line 1031)
                  L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

   sym/subsasgn (line 868)
              C = privsubsasgn(L,R,inds{:});

   normcdf>localnormcdf (line 100)
      p(sigma==0 & x<mu) = 0;

   normcdf (line 46)
  [varargout{1:max(1,nargout)}] = localnormcdf(uflag,x,varargin{:});

My guess is that Matlab is trying to test sigma==0 & x<mu, but it cannot because the symbolic variable r is included in sigma and mu. However, sigma is always >0 because of the range of r I have set.

How can I solve the equation in this case?

Answer 1

[Walter Roberson]

normpdf() and normcdf() are numeric routines. They are coded in a way that cannot succeed if the mu or sigma are symbolic. They internally need to test sigma == 0 & mu < 0 and also mu <= 0, which will fail for symbolic values that have insufficient assertions active (and even then, the assertions might not be enough.) It is okay to use symbolic values for the x for normpdf and normcdf though.

The above is why the vpasolve() is failing; it is failing while looking at the normpdf() expression before even getting into the vpasolve() call.

The solution is to turn it into a numeric system: transform vpasolve(A(r) == B(r),r) into @(r) A(r) - B(r) and apply a numeric solver to that.

Your vpasolve() is for r from 0 to 1. There is no solution to your system in that range. The solution is at approximately r = 1.008849984 .

[ByungChul Min]

Thank you for your answer. I am not a experienced Matlab user, so can you please show me how you did it?

[ByungChul Min]

Instead of using normcdf and normpdf, I directly used the original function form substituted the mu and sigma into it. It yields r of 0.9809. Thank you for your answer

[Walter Roberson]

I copied directly from your code; r = 0.9809 does not satisfy the equations you provided. Possibly you had mis-entered the equations.

f = @(r) ((c*q-equity)*(1+r)/p) - ((-p/(1+r)+p*(1-normcdf((c*q-equity)*(1+r)/p,(c*mu-equity)*(1+r)/p,c*(1+r)/p*sigma)))*(c*(1-normcdf(q,mu,sigma))-(c*q-equity)*normpdf(q,mu,sigma))/((p/(1+r))^2*normpdf(q,mu,sigma)-c^2*normpdf((c*q-equity)*(1+r)/p,(c*mu-equity)*(1+r)/p,c*(1+r)/p*sigma)));

g = @(r) f(r).^2;
r0 = fminsearch(g, 1);

Here I use the old trick of finding a zero crossing by finding the minimum of f^2. I use this rather than use fzero() because the function has a singularity not much further along.

[ByungChul Min]

You are right. I reentered the equations and it yields 1.0088499837318395856408398200666, just as what you said. Thank you very much for your help.