Generalized Owen's T function

Question

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<Y<a X$ 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. \begin{equation} 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) \end{equation} 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 \begin{equation} \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} \end{equation} 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?

Answer 1

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 : \begin{equation} 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) \end{equation} as it should be (see An integral involving error functions and a Gaussian).