Generalized Owen's T function

Solution 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).