Computing wedge probabilities: finite time horizon case

TL;DR

This paper introduces an alternative infinite series formula for computing the probability that a Brownian motion stays within sloping boundaries over a finite time, offering faster convergence and high precision.

Contribution

It provides a new series representation for wedge probabilities with improved convergence properties and practical rules for choosing the optimal formula based on parameters.

Findings

Convergence rate depends on boundary slopes and time.

Six terms of the series achieve precision of 10^{-16}.

The new formula offers an efficient alternative to Anderson's formula.

Abstract

We present an alternative to the well-known Anderson's formula for the probability that a first exit time from the planar region between two slopping lines -a_1 t -b_1 and a_2 t + b_2 by a standard Brownian motion is greater than T. As the Anderson's formula, our representation is an infinite series from special functions. We show that convergence rate of both formulas depends only on terms (a_1 + a_2)(b_1 + b_2) and (b_1 + b_2)^2 /T and deduce simple rules of appropriate representation's choose. We prove that for any given set of parameters a_1, b_1, a_2, b_2, T the sum of first 6 terms ensures precision 10^{-16}.

Figures

Equations

$$T^s=\inf \{t\ge 0,X_t\notin [-a_1,a_2]\},T^w=\inf \{t\ge 0,X_t\notin [-a_{1}t-b_1,a_{2}t+b_2]\},$$

$${\mathbb{P}}_{x_0}(T^s>T)={\int }_{-a_1}^{a_2}{u}_{x_0}^s(x,T)dx,{\mathbb{P}}_{x_0}(T^w>T)={\int }_{-a_{1}T-b_1}^{a_{2}T+b_2}{u}_{x_0}^w(x,T)dx.$$

$$\left\{\begin{array}[]{l}\left(\mathscr{L^{*}}-\partial/\partial_{t}\right)u_{x_{0}}^{s}=0,\\ u_{x_{0}}^{s}(-a_{1},t)=0,\\ u_{x_{0}}^{s}(a_{2},t)=0,\\ u_{x_{0}}^{s}(x,0)=\delta(x-x_{0}).\end{array}\right.\quad\quad\quad\quad\left\{\begin{array}[]{l}\left(\mathscr{L^{*}}-\partial/\partial_{t}\right)u_{x_{0}}^{w}=0,\\ u_{x_{0}}^{w}(-a_{1}t-b_{1},t)=0,\\ u_{x_{0}}^{w}(a_{2}t+b_{2},t)=0,\\ u_{x_{0}}^{w}(x,0)=\delta(x-x_{0}).\end{array}\right.$$

$${\mathcal{L}}^{\ast }=\mathcal{L}=\frac{1}{2}\frac{{\partial }^2}{\partial x^2}.$$

$$\displaystyle\begin{array}[]{l}u_{x_{0}}^{s}(x,t)=\frac{2}{l}\sum_{n=0}^{\infty}\sin\left(\frac{n\pi}{l}(x+a_{1})\right)sin\left(\frac{n\pi}{l}(x_{0}+a_{1})\right)e^{-\frac{n^{2}\pi^{2}}{2l^{2}}t}\\ u_{x_{0}}^{s}(x,t)=\frac{1}{\sqrt{2\pi t}}\sum_{n=-\infty}^{+\infty}\left[e^{-\frac{(x-x_{0}+2nl)^{2}}{2t}}-e^{-\frac{(x+x_{0}+2a_{1}+2nl)^{2}}{2t}}\right].\end{array}$$

$$u_{x_0}^s(x,t)=\sum _{n=1}^{\infty }{A}_n\sin \left({\omega }_n(x+\alpha )\right)e^{-\frac{{\omega }_n^{2}t}{2}}.$$

$$\displaystyle\left\{\begin{array}[]{l}\frac{1}{2}\frac{d^{2}w}{dx^{2}}-\zeta w=-\delta(x-x_{0}),\\ w(-a_{1};\zeta)=0,\\ w(a_{2};\zeta)=0.\end{array}\right.$$

$$w(x;\zeta )={\int }_0^{\infty }{e}^{-\zeta t}{u}_{x_0}^s(x,t)dt,u_{x_0}^s(x,t)=\frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }{e}^{\zeta t}w(x;\zeta )d\zeta .$$

$$\sum _{n=-\infty }^{+\infty }{e}^{-\frac{{\pi }^2{n}^2}{2u}}\cos \left(\pi nv/u\right)=\sqrt{\frac{2u}{\pi }}{e}^{-\frac{v^2}{2u}}\sum _{n=-\infty }^{+\infty }{e}^{-2n^{2}u}\cosh \left(2nv\right).$$

$$u_{x_0}^w(x,t)=\frac{a_{+}{e}^{\frac{(a_{+}{x}_0+d{)}^2}{2a_{+}{b}_{+}}-\frac{(a_{+}x-d{)}^2}{2a_{+}(a_{+}t+b_{+})}}}{\sqrt{b_{+}}\sqrt{a_{+}t+b_{+}}}{u}_{(a_{+}{x}_0-d)/b_{+}}^s\left(\frac{a_{+}x-d}{a_{+}t+b_{+}},\frac{1}{b_{+}/a_{+}}-\frac{1}{t+b_{+}/a_{+}}\right).$$

$$a_{+}=\frac{a_1+a_2}{2},b_{+}=\frac{b_1+b_2}{2},d=\frac{a_1{b}_2-a_2{b}_1}{2}.$$

$$\frac{a_{+}x-d}{a_{+}t+b_{+}}=-a_1,\frac{a_{+}x-d}{a_{+}t+b_{+}}=a_2.$$

$$u_{x_0}^w(x,0)$$

$$\Psi (\alpha ,\beta )=\frac{1}{\sqrt{2\pi }}{\int }_{\alpha }^{\beta }{e}^{-x^2/2}dx.$$

$$erf(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{-t^2}dt,\Phi (x)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^x{e}^{-t^2/2}dt.$$

$$\Psi (\alpha ,\beta )=\Phi (\beta )-\Phi (\alpha ),\Psi (\alpha ,\beta )=\frac{erf(\beta /\sqrt{2})-erf(\alpha /\sqrt{2})}{2}.$$

$$\Psi (\alpha ,\beta )=\Psi (-\beta ,-\alpha ),\Psi (-\infty ,+\infty )=1.$$

$$k(a_1,b_1;a_2,b_2;T)=\mathbb{P}(T^w>T),k(a_1,b_1;a_2,b_2)=\mathbb{P}(T^w=\infty ),$$

$$\begin{array}[]{l}k(a_{1},b_{1};a_{2},b_{2};T)=\Psi\left(\frac{-a_{1}T-b_{1}}{\sqrt{T}},\frac{a_{2}T+b_{2}}{\sqrt{T}}\right)\\ -\sum_{n=1}^{\infty}e^{-2A_{n}}\Psi\left(\frac{-a_{1}T-b_{1}+2b_{1}-4nb_{+}}{\sqrt{T}},\,\frac{a_{2}T+b_{2}+2b_{1}-4nb_{+}}{\sqrt{T}}\right)\\ -\sum_{n=1}^{\infty}e^{-B_{n}}\Psi\left(\frac{-a_{1}T-b_{1}+2b_{1}+4(n-1)b_{+}}{\sqrt{T}},\,\frac{a_{2}T+b_{2}+2b_{1}+4(n-1)b_{+}}{\sqrt{T}}\right)\\ +\sum_{n=1}^{\infty}e^{-2C_{n}}\Psi\left(\frac{-a_{1}T-b_{1}+4nb_{+}}{\sqrt{T}},\,\frac{a_{2}T+b_{2}+4nb_{+}}{\sqrt{T}}\right)\\ +\sum_{n=1}^{\infty}e^{-2D_{n}}\Psi\left(\frac{-a_{1}T-b_{1}-4nb_{+}}{\sqrt{T}},\,\frac{a_{2}T+b_{2}-4nb_{+}}{\sqrt{T}}\right)\end{array}$$

$$\begin{array}[]{c}A_{n}=n^{2}a_{2}b_{2}+(n-1)^{2}a_{1}b_{1}+n(n-1)(a_{2}b_{1}+a_{1}b_{2}),\\ B_{n}=(n-1)^{2}a_{2}b_{2}+n^{2}a_{1}b_{1}+n(n-1)(a_{2}b_{1}+a_{1}b_{2}),\\ C_{n}=n^{2}(a_{1}b_{1}+a_{2}b_{2})+n(n-1)a_{2}b_{1}+n(n+1)a_{1}b_{2},\\ D_{n}=n^{2}(a_{1}b_{1}+a_{2}b_{2})+n(n+1)a_{2}b_{1}+n(n-1)a_{1}b_{2},\end{array}$$

$$\frac{1}{b_{+}/a_{+}}-\frac{1}{T+b_{+}/a_{+}}=\frac{a_{+}^2}{b_{+}(a_{+}T+b_{+})}$$

$$k(a_1,b_1;a_2,b_2;T)=\frac{e^{\frac{d^2}{2a_{+}{b}_{+}}}}{\sqrt{2\pi T}}{\int }_{-a_{1}T-b_1}^{a_{2}T+b_2}\sum _{n=-\infty }^{+\infty }{e}^{-\frac{(a_{+}x-d{)}^2}{2a_{+}(a_{+}T+b_{+})}}\left[e^{-E_1(n)}-e^{-E_2(n)}\right]dx.$$

$$E_1(n)={\left(\frac{a_{+}x-d}{a_{+}T+b_{+}}+\frac{d}{b_{+}}+4n{a}_{+}\right)}^2,E_2(n)={\left(\frac{a_{+}x-d}{a_{+}T+b_{+}}-\frac{d}{b_{+}}+2a_1+4n{a}_{+}\right)}^2\frac{b_{+}(a_{+}T+b_{+})}{2a_{+}^{2}T};$$

$$E_1(n)+\frac{(a_{+}x-d{)}^2}{2a_{+}(a_{+}T+b_{+})}-\frac{d^2}{2a_{+}{b}_{+}}=\frac{(x+4n{b}_{+}{)}^2}{2T}+4nd+8n^2{a}_{+}{b}_{+}$$

$$E_2(n)+\frac{(a_{+}x-d{)}^2}{2a_{+}(a_{+}T+b_{+})}-\frac{d^2}{2a_{+}{b}_{+}}$$

$$\begin{array}[]{c}k(a_{1},b_{1};a_{2},b_{2};T)=\sum_{n=-\infty}^{+\infty}e^{-2((2n)^{2}a_{+}b_{+}+2nd)}\Psi\left(\frac{-a_{1}T-b_{1}+4nb_{+}}{\sqrt{T}},\frac{a_{2}T+b_{2}+4nb_{+}}{\sqrt{T}}\right)\\ -\sum_{n=-\infty}^{+\infty}e^{-2((2n)^{2}a_{+}b_{+}+n(a_{1}b_{2}+b_{1}a_{2})+(2n+1)a_{1}b_{1})}\Psi\left(\frac{-a_{1}T-b_{1}+4nb_{+}+2b_{1}}{\sqrt{T}},\frac{a_{2}T+b_{2}+4nb_{+}+2b_{1}}{\sqrt{T}}\right).\end{array}$$

$$\sum _{n=-\infty }^{+\infty }F(n)=\sum _{n=1}^{+\infty }F(n)+\sum _{n=1}^{+\infty }F(-n)+F(0),\sum _{n=-\infty }^{+\infty }F(n)=\sum _{n=1}^{+\infty }F(n-1)+\sum _{n=1}^{+\infty }F(-n)$$

$$\begin{array}[]{l}A_{n}=4n^{2}a_{+}b_{+}-n(a_{1}b_{2}+b_{1}a_{2})-(2n-1)a_{1}b_{1},\\ B_{n}=4(n-1)^{2}a_{+}b_{+}+(n-1)(a_{1}b_{2}+b_{1}a_{2})+(2n+1)a_{1}b_{1},\\ C_{n}=4n^{2}a_{+}b_{+}+2nd,\\ D_{n}=4n^{2}a_{+}b_{+}-2nd.\end{array}$$

$$k(a_1,b_1;a_2,b_2)=1-\sum _{n=1}^{\infty }{e}^{-2A_n}+e^{-2B_n}-e^{-2C_n}-e^{-2D_n}.$$

$$\begin{array}[]{c}k(a_{1},b_{1};a_{2},b_{2};T)=i\sqrt{\frac{\pi}{2a_{+}b_{+}}}e^{\frac{d^{2}}{2a_{+}b_{+}}}\sum_{n=-\infty}^{+\infty}e^{-\frac{\pi^{2}n^{2}}{8a_{+}b_{+}}-\frac{in\pi a_{1}}{2a_{+}}}\sin\left(\frac{\pi nb_{1}}{2b_{+}}\right)\times\\ \times\Psi\left(\frac{i\pi n-2a_{1}(a_{+}T+b_{+})}{2a_{+}\sqrt{T+b_{+}/a_{+}}},\frac{i\pi n+2a_{2}(a_{+}T+b_{+})}{2a_{+}\sqrt{T+b_{+}/a_{+}}}\right)\end{array}$$

Full text

1. 1 Introduction

2. 2 PDE and Lie symmetries approach to analysis of stopping times $T^{s}$ and $T^{w}$

3. 3 Explicit formulas for wedge probability $\mathbb{P}(T^{w}>T)$

4. 3.1 Some auxiliary functions and terms

5. 3.2 Explicit formulas for wedge probability

6. 4 Computational aspects

7. 4.1 Convergence rates

8. 4.2 Convergence analysis and implementation

1 Introduction

Let $W=\{W_{t},\,t\geq 0\}$ be a standard Brownian motion. We consider the following exit times

[TABLE]

where process $X$ is defined as $X=\{x_{0}+W_{t},\,t\geq 0\}$. Let’s denote by $\mathbb{P}_{x_{0}}$ probability conditional on the process $X_{t}$ started at $X_{0}=x_{0}$. We omit subscript $x_{0}$ if process $X_{t}$ starts at zero, i.e. if $X_{t}$ is a standard Brownian motion. The probability $\mathbb{P}(T^{w}>T)$ is called wedge probability and have been studied by many authors [2], [3]. Anderson [2] was the first to obtain the explicit representation for the $\mathbb{P}(T^{w}>T)$ in terms of infinite series from normal c.d.f. functions. The explicit formula for $\mathbb{P}(T^{w}=\infty)$ was found by Doob [4]. Based on properties of Jacobi theta functions, Ycart and Drouilhet [6] have found alternative to Doob’s formula. They also computed uniform precision estimates and proposed efficient numerical algorithm. Let us mention that connection of $\mathbb{P}(T^{w}=\infty)$ with theta functions has been pointed out by Salminen and Yor [5].

In this paper we generalize results from [6] to the finite time horizon case. Based on Lie symmetries for the heat equation, we deduce alternative to Anderson’s formula. Our representation is an infinite series from Error functions with complex argument. For numerical computations we derive another representation containing only real functions. Combining our results with Anderson’ formula we propose simple numerical algorithm for computation of $\mathbb{P}(T^{w}>T)$. We show that convergence rate of these two series depends only on two terms: $(a_{1}+a_{2})(b_{1}+b_{2})$ and $(b_{1}+b_{2})^{2}/T$. Anderson’ formula converges fast if at least one term is relatively high, our formula has fast convergence for the opposite case. The algorithmic consequence is that computing at most six terms of the series either in Anderson’s formula or in the new alternative suffices to approximate $\mathbb{P}_{x_{0}}(T^{w}>T)$ with precision smaller than $10^{-16}$.

The rest of paper is organized by the following scheme: in Section 2 we show connections of $T^{s}$ with $T^{w}$ in terms of corresponded PDE boundary-value problems. Based on Lie symmetries for the heat equation we derive the solutions of these PDEs in explicit form. Section 3 contains explicit formulas for probability $\mathbb{P}(T^{w}>T)$. We re-derive famous Anderson’s formula and present two new representations. Also we show that Doob’s and Ycart and Drouilhet’ formulas can be derived as the limiting form $T\rightarrow\infty$. In the last Section 4 we present upper bounds for the remainders of infinite series from Anderson’ formula and its alternative. Based on these results, we deduce rules describing in which cases Anderson’s formula or its alternative should be used and propose simple numerical algorithm for computation of wedge probability $\mathbb{P}(T^{w}>T)$.

2 PDE and Lie symmetries approach to analysis of stopping times $T^{s}$ and $T^{w}$

From standard results in probability theory $\mathbb{P}_{x_{0}}(T^{s}>T)$ and $\mathbb{P}_{x_{0}}(T^{w}>T)$ can be represented in the following form

[TABLE]

Functions $u_{x_{0}}^{s}(x,t)$ and $u_{x_{0}}^{w}(x,t)$ solve the following Cauchy problems with killed boundary conditions for Fokker–Planck–Kolmogorov equation:

[TABLE]

Here $\mathscr{L^{*}}$ is the adjoint of $\mathscr{L}$, which is the infinitesimal generator of the process $X_{t}$. In case of standard Brownian motion operator $\mathscr{L}$ is self-adjoint, i.e.

[TABLE]

Initial condition $\delta(x-x_{0})$ is a Dirac measure at the point $x_{0}$. In the next Proposition we present explicit formulas for $u_{x_{0}}^{s}$.

Proposition 1**.**

Function $u_{x_{0}}^{s}(x,t)$ has two equivalent representations ($l=a_{1}+a_{2}$):

[TABLE]

Proof.

First formula can be derived from the following well-known representation

[TABLE]

The phase shift $\alpha$ and frequencies $\omega_{n}$ are set to satisfy boundary conditions $u_{x_{0}}^{s}(-a_{1},t)=u_{x_{0}}^{s}(a_{2},t)=0$, i.e, $\alpha=a_{1}$, and $\omega_{n}=\pi n/l$. Hence at the moment $t=0$ function $u_{x_{0}}^{s}(x,0)$ is Fourier series on the segment $[-a_{1},a_{2}]$. Calculation of coefficients $A_{n}$ turns out to the main formula (4).

Second formula in (4) can be obtained by using of Laplace transform with respect to the time variable $t$. This integral transform reduces original initial boundary problem to the simple boundary problem for linear ODE:

[TABLE]

Here $w(x;\zeta)$ is an image of transformation, i.e.

[TABLE]

The problem (8) can be easily solved by a standard techniques. Inversion of Laplace transform yields second formula (4). ∎

Remark 1**.**

We can also derive identity between two representations from (4) by using Poisson’ summation formula

[TABLE]

In he next proposition we establish connection of $u_{x_{0}}^{w}$ with $u_{x_{0}}^{s}$.

Proposition 2**.**

Function $u_{x_{0}}^{w}$ can be explicitly represented in terms of $u_{x_{0}}^{s}$:

[TABLE]

here constants $a_{+}$, $b_{+}$ and $d$ are equal

[TABLE]

Proof.

This fact follows from Lie symmetries for the heat equation. One can check by the direct calculations that function $u_{x_{0}}^{w}$ from (9) solves equation $\left(\mathscr{L^{*}}-\partial_{t}\right)u_{x_{0}}^{w}=0$. Note that function $u_{x_{0}}^{s}$ is equal to zero if $x=-a_{1}$ or $x=a_{2}$, therefore $u_{x_{0}}^{w}$ equals zero if $x$ solves one of the following equations:

[TABLE]

Hence the boundary conditions $u_{x_{0}}^{w}(-a_{1}t-b_{1},t)=u_{x_{0}}^{w}(a_{2}t+b_{2},t)=0$ are satisfied. Now we check the initial condition $u_{x_{0}}^{w}(x,0)$:

[TABLE]

∎

3 Explicit formulas for wedge probability $\mathbb{P}(T^{w}>T)$

3.1 Some auxiliary functions and terms

Let’s define the following function $\Psi(\alpha,\beta)$:

[TABLE]

and recall well known Error function $erf(x)$ and normal c.d.f. $\Phi(x)$:

[TABLE]

Function $\Psi(\alpha,\beta)$ can be represented in terms of $erf(x)$ and $\Phi(x)$

[TABLE]

Let us note that function $\Psi(\alpha,\beta)$ also has the following properties:

[TABLE]

3.2 Explicit formulas for wedge probability

We will use the following notations for wedge probabilities $\mathbb{P}(T^{w}>T)$ and $\mathbb{P}(T^{w}=\infty)$:

[TABLE]

Proposition 3** (Anderson).**

Wedge probability $k(a_{1},b_{1};a_{2},b_{2};T)$ has the following representation:

[TABLE]

Here $A_{n},B_{n},C_{n},D_{n}$ are defined by the following formulas:

[TABLE]

and function $\Psi$ is defined in (11).

Proof.

Employ identity

[TABLE]

in the (9) and choose the first formula for $u_{x_{0}}(x,t)$ from (4) with $x_{0}=0$.

[TABLE]

Here

[TABLE]

One can check that

[TABLE]

and

[TABLE]

In the result we have

[TABLE]

Rearrangement of the first and second sum

[TABLE]

yields formula (16). It is easy to check that (see also [6])

[TABLE]

∎

Remark 2**.**

Doob’ formula

[TABLE]

can be easily derived from (16) if we tends $T\rightarrow\infty$. Note that if $a_{1}$ and $a_{2}$ have opposite signs then $k(a_{1},b_{1};a_{2},b_{2};\infty)=0$.

Proposition 4**.**

Formula (16) has the following alternatives:

[TABLE]

or

[TABLE]

Here function $\Psi$ is defined in (11).

Proof.

From Propositions 1 and 2 we have the following formula for $k(a_{1},b_{1};a_{2},b_{2};T)$:

[TABLE]

Change of integration variable $\varphi=a_{1}+(a_{+}x-d)/(a_{+}T+b_{+})$ yields formula (22). Definite integral from (22) is known (see [1]) and can be expressed in terms of Error functions (from complex argument):

[TABLE]

or (see (13))

[TABLE]

Using the last formula in (22) and changing the sign of $n$ in the second term gives the representation (21). ∎

Remark 3**.**

Ycart and Drouilhet formula can be derived from (21) by $T\rightarrow\infty$:

[TABLE]

Here $c$ is defined as $(a_{1}b_{1}-a_{2}b_{2})/2$ (see [6]).

4 Computational aspects

4.1 Convergence rates

Proposition 5**.**

Let’s denote by $K_{1,N}$ and $K_{2,N}$ the partial sums in series (16) and (22) respectively and consider the remainders:

[TABLE]

For $N>1$:

[TABLE]

Proof.

At first we note that $A_{n}$, $B_{n}$, $C_{n}$, $D_{n}$ are larger than $4(n-1)^{2}a_{+}b_{+}$ (see [6]) and apply second property of function $\Psi$ from (14):

[TABLE]

If quantity $a_{+}T/b_{+}$ is sufficiently large, then we can use Ycart and Drouilhet bounds [6] (this is second line in formula (25)) for the remainder $R_{1,N}$. If this term is sufficiently small, then we derive upper bound in the assumption $N>a_{+}T/b_{+}$.

Now we show that for any $n\geq N+1$ the arguments of any function $\Psi(\alpha,\beta)$ from the formula above are positive:

[TABLE]

for the second and third line:

[TABLE]

and for fourth line we apply the following inequality:

[TABLE]

Therefore,

[TABLE]

Application of the following inequality

[TABLE]

yields the following upper bound for the remainder $R_{1,N}$ :

[TABLE]

For the remainder $R_{2,N}$ we apply the following bounds:

[TABLE]

∎

In case of $T=\infty$ the bounds for remainders depend only on $a_{+}b_{+}$. In the finite horizon case we have one extra term $b_{+}^{2}/T$. We show that remainders $R_{1,N}$ and $R_{2,N}$ are characterized only by these terms. Denote by

[TABLE]

and employ these new variables in formulas (25) and (26):

[TABLE]

4.2 Convergence analysis and implementation

In this subsection we present numerical experiments illustrating convergence rates for the series (16) and (22). As we mentioned before, convergence rates depend only on values of $\alpha$ and $\beta$ which are defined in (27). Figure 1 illustrates the log-log plot of convergence rates. Color indicates minimal value of $N$ such that remainders $R_{1,N}$ and $R_{2,N}$ are less than $10^{-16}$. We can see that our formula convergences faster in case of $\alpha<1$ and $\beta<1$. In other cases Anderson’s formula converges faster than alternative.

Left sub-figure in 2 illustrates dependencies of minimal value $N$ on values of $\alpha$ and $\beta$. Let us mention, that minimal value of $N$ is equal to 5 for Anderson’s formula and to 6 for the alternative. Formula (21) should be used in the region defined by the curve which separates blue and red zones (see left Figure 2). Implementation of this rule can be difficult, so we suggest the following rule: we use our formula if and only if $\log(\alpha)<0$ and $\log(\beta)<0$. Values of $N$ are presented in right Figure 2. In this case minimal value of $N$ is increases up to 8 terms.

We implement our algorithm in C++11. All terms in formula (16) have been coded by using of functions from the standard library (std::exp and std::erf). Computation of definite integral in formula (22) have been implemented by the help of Simpson’ formula. Numerical experiments have been made using vectors of simulated entries with uniform distribution on $[0,10]$ for $a_{1}$, $b_{1}$, $a_{2}$, $b_{2}$ and on $[0,100]$ for $T$. The running time on a standard laptop is approx. 7 seconds for $10^{7}$ values.

Remark 4**.**

If we rearrange our formula (22) in the manner of Ycart and Drouilhet formula [6, formula 4] (i.e. apply product-to-sum identity for sines), we shall cut $N$ in half. Hence, in these terms we have $N=3$ (i.e. is equal to the infinite horizon case) for our formula.