# Inverting the cumulative probability function to find roots of stochastic function

Given a function:

$$f[x]=a\, \Phi \left[-x+\sigma \sqrt{\tau}\right]-\left(b+c\, e^{-d \tau}\right)\Phi \left[-x\right]$$

where $\Phi$ is the cumulative density function of the standard normal distribution: $$\Phi\left[z\right] = \frac{1}{\sqrt{2 \pi}}\int^{z}_{-\infty}e^{\frac{-u^2}{2}} \, \mathbb{d}u$$

...how can I find $x$ which which satisifies the conditon $f[x]=0$? Suppose that a, b, c, $\sigma$, t are known quantities.

I am stuck trying to use inverse identities since approximations to inverse functions seem to not hold when inverting a probability function multiplied by a constant.

Also, although there is an algorithm to find $x$ through recursion:

$$x \to-\Phi^{-1}[\frac{a}{b+c\, e^{-d \tau}}\, {\Phi\left[-x+\sigma \sqrt{\tau}\right]}] \,\,\, \forall \,\,\, \tau \in \, [t,T]$$

where $\Phi^{-1}$ is the probit function.

...this does not satisfy a closed form requirement.

Acceptable answers may include closed form solutions as well as numerical approximations provided that approximations converge $\forall_{\left|x\right|\lt ~5} \in \mathbb{R}$. I also appreciate any direction or references.

Update:

Following Synia's comments, it is apparent that $g'[x] = 0$ has at least one solution (with $g(x) = a \, \Phi [-x+s] - \Phi[-x]$), i.e.

$$g'[x] = -a\, \Phi'[-x+s] + \Phi'[-x] = \frac{1}{\sqrt{2\pi}} \left( -a e^{(x-s)^2/2 } + e^{ x^2/2 } \right)$$ vanishes in $$x^* : (a, s) \mapsto \frac{\frac{s^2}{2}-\ln\left[a\right]}{s} \quad \| \, a>0$$

Thus, if $g[x^*(a, s)] > 0$, there is no zero, if $g[x^*(a, s)] = 0$, this is the only zero, and if $g[x^*(a, s)] < 0$, there are two zeros, hence, one has to choose which one to approximate.

• It seems to me that what you are looking for is a general formula for the inverse bijection of a function $f$. You can then try the Lagrange reversion theorem en.wikipedia.org/wiki/Lagrange_reversion_theorem. Of course, you have to be able to find a formula for the $k$-th derivative of a power of your function... – Synia Aug 22 '17 at 10:21
• This will be clearer if you state up front that you are looking for good numerical techniques rather than expecting a closed form, and also replace G, B, F, X with a, b, c, x. The use of F as a constant is especially confusing in a probabilistic context. – Matt F. Aug 22 '17 at 15:37
• In fact, there are lots of constants for what you want ; you could instead consider the function $x \mapsto a \Phi(-x + s) - \Phi(-x)$ and change your variables accordingly afterwards (by scaling the equation, there are only two parameters). – Synia Aug 22 '17 at 17:02
• Nevertheless, a study of the equivalent function $g : x \mapsto a \Phi(x + s) - \Phi(x)$ shows that the derivative $g'$ can vanish in $x^* := \frac{ s^2/2 - \log(a) }{s}$ if $a > 0$ (what you certainly suppose). Hence, there are two solutions to your equations (or 1 if $g(x^*) = 0$ or 0 if $g(x^*) > 0$). You must thus specify the interval on which you want a solution in the case $g(x^*) < 0$. – Synia Aug 22 '17 at 17:17
• @Synia. Thank you for your input. I was not familiar with Lagrange reversion. I will check into that. Also, I see how your logic defines that a solution exists, but do not understand how to take it any further. – David Addison Aug 23 '17 at 4:36

Since you are considering a numerical approach: If you have access to Mathematica, you can use the "Reduce" functionality to find all roots in a given interval.

In[1]:= f[x_] := (a*Erf[-Infinity, (-x + sigma*Sqrt[tau])*Sqrt[2]]
- (b + c*Exp[-d*tau])*Erf[-Infinity, -x*Sqrt[2]])/2
In[2]:= a = 1; b = 1; c = 1; d = 2.5; sigma = 4; tau = 2.5;
In[3]:= Reduce[f[x] == 0 && -5 < x < 5]
Out[3]: x == -1.44496

I checked that it also works on Wolfram Alpha (so a Mathematica license is not really needed).

This documentation explains how Reduce tries to find all roots:

Given rough locations for roots, we use Newton-like iterative methods to home in on the roots. And finally, we use Mathematica‘s interval arithmetic to prove that there can only be one root inside each small region we’ve identified. There is some subtlety here. Proving zero equivalence is in general undecidable—and the way this shows up is that in pathological cases it can in principle require unboundedly much precision to distinguish multiple roots from closely spaced single roots.

• Thank you for your response. I really hadn't used Reduce numerically before, but it looks a lot more efficient than "FindRoot". However, I wasn't really looking for a numerical approach. I intended to mean that a rational (i.e., polynomial) approximation would be accepted. – David Addison Aug 28 '17 at 21:42