# Generalized Owen's T function

As Wikipedia teaches us https://en.wikipedia.org/wiki/Owen%27s_T_function the Owen's T function $$T(h,a)$$ defines a probability of a bivariate event $$X>h$$ and $$0 where $$X,Y$$ are standard, independent Gaussian random variables.

Now in the context of question Multivariate gaussian integral over positive reals a necessity appeared to deal with a slightly more general quantity. $$$$T(h,a,b):= {\bf P}\left(X>h \quad \wedge \quad a X+b > Y > 0 \left. \right| X = N(0,1) , Y=N(0,1) \right)$$$$ We have shown that : $$\begin{eqnarray} &&T(h,a,b)= \int\limits_h^\infty \frac{\exp(-1/2 \xi^2)}{\sqrt{2\pi}} \frac{1}{2} Erf(\frac{a \xi+b}{\sqrt{2}}) d\xi \quad (i1)\\ &&= \int\limits_0^a \frac{e^{-\frac{b^2}{2}-b h \xi -\frac{1}{2} h^2 \left(\xi ^2+1\right)}}{2 \pi \left(\xi ^2+1\right)} d\xi - \frac{b}{2\sqrt{2}\sqrt{\pi}} \int\limits_0^a \frac{\xi e^{-\frac{b^2}{2 \xi ^2+2}} \text{erfc}\left(\frac{\xi (b+h \xi )+h}{\sqrt{2} \sqrt{\xi ^2+1}}\right)}{\left(\xi ^2+1\right)^{3/2}} d\xi + \frac{1}{4} \text{erf}\left(\frac{b}{\sqrt{2}}\right) \text{erfc}\left(\frac{h}{\sqrt{2}}\right) \quad (i2) \end{eqnarray}$$

{a, b, h} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50]; b = 0;
NIntegrate[
Exp[-x^2/2]/Sqrt[2 Pi] 1/2 Erf[(a x + b)/Sqrt[2]], {x, h, Infinity},
WorkingPrecision -> 20]
NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
2 (1 + xi^2) \[Pi]) -
b  /(2 Sqrt[2] Sqrt[ \[Pi]]) (
xi  Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 +
xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a},
WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4


Update: Let $$A_j \in {\mathbb R}$$ for $$j=1,\cdots,3$$ and let $$x\in {\mathbb R}$$. Then we have: $$\begin{eqnarray} T(A_1 x, A_2, A_3 x) = \frac{1}{2\pi} \left(\arctan(A_2)-\arctan(A_2+\frac{A_3}{A_1})-\arctan(\frac{A_1+A_2 A_3+A_2^2 A_1}{A_3})\right) + \frac{1}{4} erf[\frac{A_3 x}{\sqrt{2} \sqrt{1+A_2^2}}] + T(A_1 x, \frac{A_2 A_1+A_3}{A_1})+ T(\frac{A_3 x}{\sqrt{1+A_2^2}},\frac{A_1+A_2 A_3 + A_2^2 A_1}{A_3}) \quad (ii) \end{eqnarray}$$

This identity comes from differentiating both sides with respect to $$x$$ then using the definition of the generalized Owen's T function to evaluate the derivative on the right hand side and having done this integrating both sides with respect to $$x$$ again.

Let us present the proof of that in detail. Firstly we define $$f(x) := T[A_1 x, A_2, A_3 x]$$. Now we compute the derivative using the chain rule. We have: $$\begin{eqnarray} \frac{d }{d x} f(x) &=& \partial_1 T[A_1 x, A_2 , A_3 x] \cdot A_1 + \partial_3 T[A_1 x, A_2, A_3 x] \cdot A_3 \\ &=& - \left. \rho(h) \frac{1}{2} erf[\frac{a h+b}{\sqrt{2}}] \right|_{\begin{array}{r} h=A_1 x \\ a=A_2 \\ b=A_3 x \end{array}}\cdot A_1 + \left.\frac{1}{\sqrt{1+a^2}} \frac{1}{2} erf[\frac{h+a^2 h+a b}{\sqrt{2} \sqrt{1+a^2}}] \rho(\frac{b}{1+a^2})\right|_{\begin{array}{r} h=A_1 x \\ a=A_2 \\ b=A_3 x \end{array}} \cdot A_3 \\ &=& -\rho(A_1 x) \frac{1}{2} erf[\frac{A_1 A_2 + A_3}{\sqrt{2}} x] \cdot A_1 + \frac{1}{\sqrt{1+A_2^2}} \rho(\frac{A_3 x}{\sqrt{1+A_2^2}}) \frac{1}{2} erfc[\frac{A_1+A_2 A_3 +A_1 A_2^2}{\sqrt{2} \sqrt{1+A_2^2}} x] \cdot A_3 \end{eqnarray}$$

Now we integrate. We have: $$\begin{eqnarray} f(x)- f(0) &=& - \int\limits_0^x \rho(A_1 \xi) \frac{1}{2} erf[\frac{A_1 A_2 + A_3}{\sqrt{2}} \xi] d\xi \cdot A_1 + \\ &&\frac{1}{\sqrt{1+A_2^2}}\int\limits_0^x \rho(\frac{A_3 \xi}{\sqrt{1+A_2^2}}) \frac{1}{2} erfc[\frac{A_1+A_2 A_3 +A_1 A_2^2}{\sqrt{2} \sqrt{1+A_2^2}} \xi] d\xi \cdot A_3 \\ f(x) - \frac{1}{2\pi} \arctan(A_2) &=& - \frac{1}{2\pi} \arctan\left( \frac{A_1 A_2+A_3}{A_1}\right) + T(A_1 x, \frac{A_1 A_2 + A_3}{A_1}) + \\ && \frac{1}{4} erf\left( \frac{A_3}{\sqrt{2} \sqrt{1+A_2^2}} x\right) +\\ &&-\frac{1}{2\pi} \arctan\left( \frac{A_1+A_2 A_3 + A_1 A_2^2}{A_3}\right) + T\left( \frac{A_3}{\sqrt{1+A_2^2}} x, \frac{A_1+A_2 A_3 + A_1 A_2^2}{A_3}\right) \end{eqnarray}$$ where in the second line we used the results from An integral involving error functions and a Gaussian and the definition of the Owen's T function. This completes the proof.

(*A certain derivative. Used in Q869502.nb*)
T[h_, a_, b_] :=
NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
2 (1 + xi^2) \[Pi]) -
b  /(2 Sqrt[2] Sqrt[ \[Pi]]) (
xi  Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 +
xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a},
WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4;
{A1, A2, A3} = RandomReal[{-1, 1}, 3, WorkingPrecision -> 50];
u = Range[0, 1, 1/100];
mT = Interpolation[Transpose[{u, T[A1 u, A2, A3 u]}]];
u =.; u = RandomReal[{0, 1}, WorkingPrecision -> 50];
mT'[u]
-rho[A1 u] 1/2 Erf[(A1 A2 + A3)/Sqrt[2] u] A1 +
1/Sqrt[1 + A2^2]
rho[(A3 u)/Sqrt[1 + A2^2]] 1/
2 Erfc[(A1 + A2 A3 + A1 A2^2)/(Sqrt[2] Sqrt[1 + A2^2]) u] A3

T[A1 u, A2, A3 u]
1/(2 Pi) (ArcTan[A2] - ArcTan[(A2 A1 + A3)/A1] -
ArcTan[(A1 + A2 A3 + A2^2 A1)/A3]) +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
1/(2 Pi) (-ArcTan[A3/((A1 + A2 A3 + A2^2 A1))] -
ArcTan[(A1 + A2 A3 + A2^2 A1)/A3]) +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]
-1/(2 Pi) Pi/2 (Sign[A3/((A1 + A2 A3 + A2^2 A1))]) +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]

-(1/4) Sign[A3/((A1 + A2 A3 + A2^2 A1))] +
1/4 Erf[(A3 u)/(Sqrt[2] Sqrt[1 + A2^2])] +
OwenT[A1 u, (A2 A1 + A3)/A1] +
OwenT[A3/Sqrt[1 + A2^2] u, (A1 + A2 A3 + A2^2 A1)/A3]


Now by both taking $$x=1$$ and replacing $$A_1$$,$$A_2$$ and $$A_3$$ by $$h$$, $$a$$ and $$b$$ in $$(ii)$$ we express the generalized Owen's T function through Owen's T function itself. We have: $$\begin{eqnarray} T(h,a,b) = \frac{1}{2\pi} \left(\arctan(a)-\arctan(a+\frac{b}{h})-\arctan(\frac{h+a b+a^2 h}{b})\right) + \frac{1}{4} erf[\frac{b}{\sqrt{2(1+a^2)}}] + T\left( h,\frac{a h+b}{h}\right) + T\left( \frac{b}{\sqrt{1+a^2}},\frac{h+a b+a^2 h}{b}\right) \end{eqnarray}$$

As a sanity check we look at the limit $$b$$ going to zero. We have: $$\begin{eqnarray} \lim_{b \rightarrow 0_+} T(h,a,b) &=& \frac{1}{2\pi} \left(\arctan(a)-\arctan(a)-\frac{\pi}{2} sign(h))\right) + 0 + T(h,a) + \frac{1}{4} sign(h) \\ &=& T(h,a) \end{eqnarray}$$ as it should be.

As another sanity check we look at the case $$a=\imath$$. Going back to the calculations of the derivative above we have: $$\begin{eqnarray} \frac{d}{d x} f(x)= -\phi(A_1 x) \frac{1}{2} erf(\frac{A_1 A_2+A_3}{\sqrt{2}} x) A_1 + \frac{1}{2\pi \imath x} \exp(-\frac{1}{2} x^2 (2 A_1 \imath A_3+A_3^2)) \end{eqnarray}$$ where we used the asymptotic expansion for the complementary error function given in https://en.wikipedia.org/wiki/Error_function#Complementary_error_function . Now we take a number $$M$$ such that $$1< M$$ and we integrate the above from unity to $$M$$ and we get: $$\begin{eqnarray} f(1)-f(M)= \left.\left( T(A_1 \cdot \xi,A_2+\frac{A_3}{A_1}) + \frac{1}{4\pi \imath} Ei(-\frac{1}{2}(1+2\imath \frac{A_1}{A_3})(\xi A_3)^2\right)\right|_{\xi=M}^{\xi=1} \end{eqnarray}$$ where $$Ei()$$ is the exponential integral. Now it turns out that as $$M\rightarrow \infty$$ both $$f(M)$$ and $$T(\dots M,\dots)$$ tend to zero and $$$$\lim\limits_{M\rightarrow \infty} \frac{1}{4 \pi \imath} Ei((a+\imath b) M) = sign(b) \cdot \frac{1}{4} \cdot 1_{a<0} + \infty \cdot 1_{a>0}$$$$ Defining $$b:=b_1+\imath b_2$$ and taking $$h>0$$ this gives the final result: $$\begin{eqnarray} &&T(h,\imath , b) = \\ &&\left\{ \begin{array}{rr} T(h,\imath + \frac{b}{h}) + \frac{1}{4\pi \imath} Ei(\frac{1}{2}(-b_1^2+b_2^2+2 b_2 h-2\imath b_1(b_2+h))) + sign(b_1(b_2+h)) \cdot \frac{1}{4} & \mbox{if b_2<0 and -b_1^2 + b_2^2+2 b_2 h <0} \\ \infty & \mbox{otherwise} \end{array} \right. \end{eqnarray}$$

My question is the following. Has this quantity been ever analyzed in the literature before?

This is not an answer to the question above but instead it is a generalization of the results above. Define $$\vec{a}:=(a_j)_{j=1}^d \in {\mathbb R}_+^d$$ and let us define as $$T^{(d)}(h,\vec{a})$$ the probability of a following multivariate event $$X>h$$ and $$0< Y_j < a_j X$$ for $$j=1,\cdots,d$$ where $$X$$ and $$\left( Y_j \right)_{j=1}^d$$ are standard, independent Gaussian random variables.

Now take another vector $$\vec{b}:=(b_j)_{j=1}^d \in {\mathbb R}_+^d$$ and define a slightly more general quantity: $$\begin{eqnarray} T^{(d)}(h,\vec{a},\vec{b}) &:=& P\left( X > h \quad \wedge \quad \begin{array}{rrr} a_1 X + b_1 >& Y_1& >0 \\ \vdots &\vdots& \vdots \\ a_d X + b_d >& Y_d& >0 \end{array}\left. \right| \begin{array}{rrr} X&=&N(0,1)\\Y_1 &=&N(0,1)\\ &\vdots& \\ Y_d&=&N(0,1) \end{array} \right)\\ &=& \int\limits_h^\infty \rho(\xi) \left[\prod\limits_{i=1}^d \frac{1}{2} erf(\frac{a_i \xi+b_i}{\sqrt{2}})\right] d\xi \end{eqnarray}$$ In what follows we will prove that if $$d\le 2$$ then the quantity $$T^{(d)}(h,\vec{a},\vec{b})$$ reduces to elementary functions and to $$T^{(d)}(h,\vec{a})$$ only.

As in the question above we consider a following quantity $$T^{(d)}(h \cdot x, \vec{a}, \vec{b} \cdot x)$$ which we differentiate with respect to $$x$$. We have: $$\begin{eqnarray} &&\frac{d}{d x} T^{(d)}(h \cdot x, \vec{a}, \vec{b} \cdot x)=\\ &&-h \cdot \rho(h \cdot x) \prod\limits_{i=1}^d \frac{1}{2} erf(\frac{a_i h x+b_i x}{\sqrt{2}})+\\ &&\sum\limits_{i=1}^d \frac{b_i}{\sqrt{2 \pi}} \int\limits_{h \cdot x}^\infty e^{-\frac{1}{2} (a_i \xi + b_i x)^2}\rho(\xi) \left[\prod\limits_{j=1,j\neq i}^d \frac{1}{2} erf(\frac{a_j h x+b_j x}{\sqrt{2}})\right] d\xi \end{eqnarray}$$ What we do now is to simplify the second term on the right hand side, i.e. absorb the exponential into the Gaussian density and extract a constant prefactor. After that we integrate the identity above over $$x$$ from zero to unity. The result reads: $$\begin{eqnarray} &&T^{(d)}(h,\vec{a},\vec{b}) - T^{(d)}(0,\vec{a},\vec{0})=\\ &&-T^{(d)}(0,\vec{a}+\frac{1}{h} \vec{b},\vec{0}) + T^{(d)}(h,\vec{a}+\frac{1}{h} \vec{b},\vec{0}) + \\ &&\sum\limits_{i=1}^d \int\limits_0^{\frac{b_i}{\sqrt{1+a_i^2}}} \rho(x) \cdot T^{(d-1)}\left([a_i + \frac{h}{b_i}(1+a_i^2)]x,\frac{(a_j)_{j=1,j\neq i}^d}{\sqrt{1+a_i^2}},\frac{(b_j(1+a_i^2)-b_i a_i a_j)_{j=1,j\neq i}^d}{b_i \sqrt{1+a_i^2}} \right) dx \end{eqnarray}$$ this clearly gave us a recurrence relation for the quantity in question subject to $$T^{(0)}(h,\vec{a},\vec{b})=T^{(0)}(h)= 1/2 erfc(h/\sqrt{2})$$.

Now we state the result for $$d=2$$. Firstly we define auxiliary quantities: $$\begin{eqnarray} \delta&:=& h^2+(a_1 h+b_1)^2+(a_2 h+b_2)^2\\ \delta_1&:=&h(1+a_1^2+a_2^2) + a_1 b_1 + a_2 b_2\\ \delta_2&:=&1+a_1^2+a_2^2\\ \hline\\ (m_1,m_2)&:=& (b_1(1+a_2^2)-a_1 a_2 b_2,b_2(1+a_1^2)-a_1 a_2 b_1)\\ (n_1,n_2)&:=&(h+h a_1^2+a_1 b_1,h+h a_2^2+a_2 b_2)\\ (o_1,o_2)&:=&(h a_1+b_1,h a_2 +b_2)\\ (p_1,p_2)&:=&(1+a_1^2,1+a_2^2)\\ (k_1,k_2)&:=&(\frac{\sqrt{p_1} \delta_1}{m_2},\frac{\sqrt{p_2} \delta_1}{m_1})\\ (l_1,l_2)&:=&(\frac{m_1}{\sqrt{p_2 \delta_2}},\frac{m_2}{\sqrt{p_1 \delta_2}}) \end{eqnarray}$$ Then the result reads: $$\begin{eqnarray} &&4 \pi T^{(2)}(h,\vec{a},\vec{b})=\\ &&\arctan(\frac{a_1 a_2}{\sqrt{\delta_2}}) - \arctan(\frac{o_1 o_2}{h \sqrt{\delta}})+\\ &&\arctan(\frac{m_2}{\sqrt{\delta_2} b_1}) + \arctan(\frac{m_1}{\sqrt{\delta_2} b_2})+\\ &&\arctan(\frac{b_1 o_2}{n_1 \sqrt{\delta}}) + \arctan(\frac{b_2 o_1}{n_2 \sqrt{\delta}})+\\ &&\arctan(\frac{b_1 \delta_1}{\sqrt{m_2^2 \delta}}) + \arctan(\frac{b_2 \delta_1}{\sqrt{m_1^2 \delta}})+\\ &&\left( \arctan(\frac{a_2}{\sqrt{p_1}}) - \arctan(\frac{\sqrt{p_1} o_2}{n_1})-\arctan(k_1)\right) \cdot erf(\frac{b_1}{\sqrt{2 p_1}})+\\ &&\left( \arctan(\frac{a_1}{\sqrt{p_2}}) - \arctan(\frac{\sqrt{p_2} o_1}{n_2})-\arctan(k_2)\right) \cdot erf(\frac{b_2}{\sqrt{2 p_2}})+\\ &&2\pi \left( T(\frac{n_1}{\sqrt{p_1}},\frac{\sqrt{p_1} o_2}{n_1}) + T(l_2,k_1)\right)\cdot erf(\frac{b_1}{\sqrt{2 p_1}})+\\ &&2\pi \left( T(\frac{n_2}{\sqrt{p_2}},\frac{\sqrt{p_2} o_1}{n_2}) + T(l_1,k_2)\right)\cdot erf(\frac{b_2}{\sqrt{2 p_2}})+\\ &&-2\pi\left( T(\frac{b_1}{\sqrt{p_1}},\frac{m_2}{\sqrt{\delta_2} b_1}) + T(\frac{b_2}{\sqrt{p_2}},\frac{m_1}{\sqrt{\delta_2} b_2})\right)+\\ &&4\pi \left( T^{(2)}(h,(a_j+\frac{b_j}{h})_{j=1}^2) -T^{(2)}(\frac{n_1}{\sqrt{p_1}},(\frac{b_1}{n_1},\frac{\sqrt{p_1} o_2}{n_1})) -T^{(2)}(\frac{n_2}{\sqrt{p_2}},(\frac{b_2}{n_2},\frac{\sqrt{p_2} o_1}{n_2})) -T^{(2)}(l_2,(\frac{\sqrt{\delta_2}b_1}{m_2},k_1)) -T^{(2)}(l_1,(\frac{\sqrt{\delta_2}b_2}{m_1},k_2)) \right) \end{eqnarray}$$ As usual I include a piece of code that verifies this expression:

d = 2; M = 3;
mj =.; mj[i_] := If[i == 1, 2, 1];
Clear[phi]; phi[x_] := Exp[-1/2 x^2]/Sqrt[2 Pi];
Clear[T];
T[h_, a_, b_] :=
1/(2 Pi) (ArcTan[a] - ArcTan[a + b/h] - ArcTan[a + h/b + a^2 h/b]) +
1/4 Erf[b/Sqrt[2 (1 + a^2)]] + OwenT[h, a + b/h] +
OwenT[b/Sqrt[1 + a^2], a + h/b + a^2 h/b];
gT[h_, a_] :=
NIntegrate[
phi[xi] Product[
1/2 Erf[a[[i]]/Sqrt[2] xi], {i, 1, Length[a]}], {xi, h,
Infinity}];
For[count = 1, count <= 100, count++,
h = RandomReal[{0, M}, WorkingPrecision -> 50];
Clear[a];
For[i = 1, i <= d, i++,
a[i] = RandomReal[{0, M}, WorkingPrecision -> 50]];
Clear[b];
For[i = 1, i <= d, i++,
b[i] = RandomReal[{0, M}, WorkingPrecision -> 50]];

I1 = NIntegrate[
phi[xi] Product[
1/2 Erf[(a[i] xi + b[i])/Sqrt[2]], {i, 1, d}], {xi, h,
Infinity}];

NIntegrate[
phi[xi] Product[1/2 Erf[a[i]/Sqrt[2] xi], {i, 1, d}], {xi, 0,
Infinity}] -
NIntegrate[
phi[xi] Product[
1/2 Erf[(a[i] + 1/h b[i] )/Sqrt[2] xi], {i, 1, d}], {xi, 0,
Infinity}] +
NIntegrate[
phi[xi] Product[
1/2 Erf[(a[i] + 1/h b[i] )/Sqrt[2] xi], {i, 1, d}], {xi, h,
Infinity}] +
Sum[NIntegrate[
phi[xi] T[(a[i] + h/b[i] (1 + a[i]^2)) xi, a[mj[i]]/Sqrt[
1 + a[i]^2], (b[mj[i]] (1 + a[i]^2) - b[i] a[i] a[mj[i]])/(
b[i] Sqrt[1 + a[i]^2]) xi], {xi, 0, b[i]/Sqrt[
1 + a[i]^2]}], {i, 1, d}];
dd = h^2 + (a[1] h + b[1])^2 + (a[2] h + b[2])^2;
dd1 = h (1 + a[1]^2 + a[2]^2) + a[1] b[1] + a[2] b[2];
dd2 = 1 + a[1]^2 + a[2]^2;

{m1, m2} = {b[1] (1 + a[2]^2 ) - a[1] a[2] b[2],
b[2] (1 + a[1]^2 ) - a[1] a[2] b[1]};
{n1, n2} = {h + h a[1]^2 + a[1] b[1], h + h a[2]^2 + a[2] b[2]};
{o1, o2} = {h a[1] + b[1], h a[2] + b[2]};
{p1, p2} = {1 + a[1]^2, 1 + a[2]^2};
{k1, k2} = {(Sqrt[p1] (dd1))/m2, (Sqrt[p2] (dd1))/m1};
{l1, l2} = {m1/Sqrt[(p2) (dd2)], m2/Sqrt[(p1) (dd2)]};

I2 = 1/(
4 \[Pi]) (ArcTan[(a[1] a[2])/Sqrt[dd2]] -
ArcTan[((o1) (o2))/(h Sqrt[dd])] +
ArcTan[m2/(Sqrt[dd2] b[1])] + ArcTan[m1/(Sqrt[dd2] b[2])] +
ArcTan[(b[1] (o2))/((n1) Sqrt[dd])] +
ArcTan[(b[2] (o1) )/((n2) Sqrt[dd])] +
ArcTan[(b[1] (dd1))/Sqrt[(m2)^2 (dd)]] +
ArcTan[(b[2] (dd1))/Sqrt[(m1)^2 (dd)]] +
(ArcTan[a[2]/Sqrt[p1]] - ArcTan[(Sqrt[p1] (o2))/n1] -
ArcTan[k1]) Erf[b[1]/(
Sqrt[2] Sqrt[p1])] + (ArcTan[a[1]/Sqrt[p2]] -
ArcTan[(Sqrt[p2] (o1))/n2] - ArcTan[k2]) Erf[b[2]/(
Sqrt[2] Sqrt[p2])] +
2 \[Pi] (OwenT[n1/Sqrt[p1], (Sqrt[p1] (o2))/n1] +
OwenT[l2, k1]) Erf[b[1]/(Sqrt[2] Sqrt[p1])] +
2 \[Pi] (OwenT[n2/Sqrt[p2], (Sqrt[p2] (o1))/n2] +
OwenT[l1, k2]) Erf[b[2]/(Sqrt[2] Sqrt[p2])] -
2 \[Pi] (OwenT[b[1]/Sqrt[p1], m2/(Sqrt[dd2] b[1])] +
OwenT[b[2]/Sqrt[p2], m1/(Sqrt[dd2] b[2])]) +
4 \[Pi] (gT[h, {a[1] + b[1]/h, a[2] + b[2]/h}] +
-gT[n1/Sqrt[p1], {b[1]/n1, (Sqrt[p1] (o2))/n1}] -
gT[n2/Sqrt[p2], {b[2]/n2, (Sqrt[p2] (o1))/n2}] -
gT[l2, {(Sqrt[dd2] b[1])/m2, k1}] -
gT[l1, {(Sqrt[dd2] b[2])/m1, k2}]));
If[Abs[I2/I1 - 1] > 10^(-2),
Print["Results do not match..", {count, {a[1], a[2], b[1], b[2],
h}, {I1, I2}}]; Break[]];
PrintTemporary[{count, I1, I2}];
];


Update: It might be interesting to know if it is possible to express the quantities $$T^{(2)}(h,(a_1,a_2))$$ in some alternative way. As a matter of fact by starting from the integral definition of this quantity , then differentiating with respect to $$a_1$$ and then integrating by parts and finally integrating with respect to $$a_1$$ from zero to $$a_1$$ we stumbled on a following formula: $$\begin{eqnarray} T^{(2)}(h,(a_1,a_2)) = \frac{2 \pi \text{erf}\left(\frac{\text{a2} h}{\sqrt{2}}\right) T(h,\text{a1})+\arctan\left(\frac{\text{a1} \text{a2}}{\sqrt{\text{a1}^2+\text{a2}^2+1}}\right) \text{erfc}\left(\frac{h \sqrt{\text{a1}^2+\text{a2}^2+1}}{\sqrt{2}}\right)}{4 \pi } + \frac{h \sqrt{1+a_2^2}}{\pi^{3/2} 2^{3/2}} \int\limits_0^{arccosh(\sqrt{\frac{1+a_1^2+a_2^2}{1+a_2^2}})} \sinh(\theta) \cdot \arctan\left(a_2 \frac{\sinh(\theta)}{\cosh(\theta)} \right) \cdot e^{-\frac{h^2}{2} (1+a_2^2) \cosh(\theta)^2} d\theta \end{eqnarray}$$

In particular for $$h=0$$ we have : $$$$T^{(2)}(0,(a_1,a_2)) = \frac{1}{4 \pi} \arctan\left(\frac{\text{a1} \text{a2}}{\sqrt{\text{a1}^2+\text{a2}^2+1}}\right)$$$$ as it should be (see An integral involving error functions and a Gaussian).