The distribution of the difference between two correlated non-central t distribution Suppose a binormal population $\{X, Y\}$ with means $\mathbf{\mu} = \{\mu_1,\mu_2\} \ne \{0,0\}$ and covariance $\Sigma= \sigma^2\begin{bmatrix}1 & \rho\\ \rho &1 \end{bmatrix}$. Let $S^2$ be an estimator of $\sigma^2$ with $f$ degrees of freedom.
It is known that the variates $\{X/S, Y/S\}$ follows a bivariate non-central $t$ distribution (e.g., Kshirsagar, 1961).
Walgren (1980) derived the distribution of the product $X/S \times Y/S$. Is there a derivation for their difference, $X/S - Y/S$?
Edit: In the present context, I use the pooled standard deviation, that is the mean of the separate standard deviations,
$$S_X = \sum_{i=1}^n (X_i-\bar{X})^2 / (n-1)$$
$$S_Y = \sum_{i=1}^n (Y_i-\bar{Y})^2 / (n-1)$$
with $ S = \sqrt{ (S_X^2 + S_Y^2)/2 }$. This estimate of $\sigma$ is independent of both $X$ and $Y$, and as shown here  is a chi-square distribution with degrees of freedom $2(n-1)/(1+\rho^2)$.

Here is an approach that shows how to obtain a numerical approximation to the probability density function of $Z=X/S-Y/S$.  (I haven't been successful in finding an analytic solution.)
Using the joint density of $X/S$ and $Y/S$ found in Kshirsagar 1961 (as given in the question):
r = {{1, \[Rho]}, {\[Rho], 1}}; (* Correlation matrix *)
\[Mu] = {\[Mu]x, \[Mu]y};
t = {x, y};
f = 2 (n - 1)/(1 + \[Rho]^2); (* Degrees of freedom for estimate of S^2 *)
jointPDF = (Exp[-\[Mu] . Inverse[r] . \[Mu]/(2 \[Sigma]^2)]/(\[Pi] f Sqrt[ Det[r]] Gamma[f/2]))*
  Sum[(2^(\[Alpha]/2) (t . Inverse[r] . \[Mu])^\[Alpha]  Gamma[(f + 2 + \[Alpha])/2])/
  (\[Sigma]^\[Alpha] f^(\[Alpha]/2) \[Alpha]! 
  (1 + t . Inverse[r] . t/f)^((f + 2 + \[Alpha])/2)), {\[Alpha], 0, \[Infinity]}];
jointPDF = FullSimplify[jointPDF, Assumptions -> {-1 < \[Rho] < 1, \[Sigma] > 0, n > 1, 
    n \[Element] Integers, \[Mu]x \[Element] Reals, \[Mu]y \[Element] 
     Reals, x \[Element] Reals, y \[Element] Reals}]

A more readable version of the code is below:

The result is

The pdf of the difference $Z = X/S - Y/S$ can be found numerically by replacing $y$ with $x-z$ and then integrating over $x$:
(* Numerical estimate of pdf of X/S - Y/S for a few values of n \
(sample size for estimating \[Sigma])*)
pdfz100 = 
  Table[{z, 
    NIntegrate[
     jointPDF /. {y -> x - z, 
       n -> 100, \[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y -> 
        3}, {x, -\[Infinity], \[Infinity]}]}, {z, -6, 3, 1/10}];
pdfz4 = Table[{z, 
    NIntegrate[
     jointPDF /. {y -> x - z, 
       n -> 4, \[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y -> 
        3}, {x, -\[Infinity], \[Infinity]}]}, {z, -6, 3, 1/10}];
pdfz2 = Table[{z, 
    NIntegrate[
     jointPDF /. {y -> x - z, 
       n -> 2, \[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y -> 
        3}, {x, -\[Infinity], \[Infinity]}]}, {z, -6, 3, 1/10}];
ListPlot[{pdfz100, pdfz4, pdfz2}, Joined -> True, ImageSize -> Large, 
 PlotLegends -> {"n = 100", "n = 4", "n = 2"},
 PlotLabel -> 
  Style["\[Sigma] = 2, \[Rho] = 1/2, \!\(\*SubscriptBox[\(\[Mu]\), \
\(x\)]\) = 1, \!\(\*SubscriptBox[\(\[Mu]\), \(y\)]\) = 3", Bold, 18]]

Again, a more readable version:

The results follow:

As a check one should perform some simulations.
parms = {\[Sigma] -> 2, \[Rho] -> 1/2, \[Mu]x -> 1, \[Mu]y -> 3, n -> 2};
nsim = 100000; (* Number of simulations *)
(* Data for x and y *)
data = RandomVariate[
  BinormalDistribution[{\[Mu]x, \[Mu]y}, {\[Sigma], \[Sigma]}, \[Rho]] /. parms, nsim];
(* Data to for estimating S *)
xy = RandomVariate[
   BinormalDistribution[{mu1, mu2}, {\[Sigma], \[Sigma]}, \[Rho]] /. 
    parms, {nsim, n /. parms}];
s = Sqrt[(Variance[#[[All, 1]]]/2 + Variance[#[[All, 2]]]/2) & /@ xy];
(* Z = X/S - Y/S *)
zz = data[[All, 1]]/s - data[[All, 2]]/s;

(* Numerically estimate the pdf of z *)
pdfz = Table[{z, NIntegrate[jointPDF /. y -> x - z /. parms, 
  {x, -\[Infinity], \[Infinity]}]},
   {z, Quantile[zz, 0.005], Quantile[zz, 0.995], (Quantile[zz, 0.995] - Quantile[zz, 0.005])/200}];

(* Plot the results *)
Show[Histogram[zz, "FreedmanDiaconis", "PDF"], 
 ListPlot[pdfz, Joined -> True, PlotRange -> All]]


There seems to be a match.
