MLE of parameters for a difference of two Exponential IID Suppose I have $X_1 \sim Exp(\theta_1)$ and $X_2\sim Exp(\theta_2)$. Then it is not difficult to show that $Y = X_1 - X_2$ will have density:
$f_Y(y) = \frac{1}{\theta_1 + \theta_2}e^{-y/\theta_1}\mathbb{1}_{y > 0} + \frac{1}{\theta_1+\theta_2}e^{y/\theta_2}\mathbb{1}_{y\leq 0}$.
Given $n$ iid observations from $Y$, we can find the sufficient statistics and the MLE of $\theta_1, \theta_2$. I am not sure if there is a shorter way to do this... I seem to arrive at quadratic equations here when I do this.
The likelihood is:
$\frac{1}{(\theta_1 + \theta_2)^n}e^{-\frac{1}{\theta_1}\sum y^+ +\frac{1}{\theta_2}\sum y^-}$; where the $y^+ = \max(0,y)$ and $y^- = \min(0,y)$. By factorization, we have $\left(\sum y^+, \sum y^-\right)$ as the sufficient statistics.
Now for the MLE, we have first the log of the likelihood which is:
$\ell = -n\log(\theta_1 + \theta_2)-\frac{1}{\theta_1}\sum y^+ + \frac{1}{\theta_2}\sum y^-$
Then the first derivatives which are set of 0, are:
$\frac{\partial \ell}{\partial \theta_1} = -\frac{n}{\theta_1 + \theta_2} + \frac{\sum y^+}{\theta_1^2}\equiv 0$ and $\frac{\partial \ell}{\partial \theta_2} = -\frac{n}{\theta_1 + \theta_2}-\frac{\sum y^-}{\theta_2^2} \equiv 0$.
So it follows that $\frac{1}{n} \sum y^+ = \left(\frac{\theta_1^2}{\theta_1+\theta_2}\right)_{MLE}$ and $-\frac{1}{n} \sum y^- = \left(\frac{\theta_2^2}{\theta_1+\theta_2}\right)_{MLE}$. This is where I am stuck.
How can we find the $\hat{\theta}_1^{MLE}$ and $\hat{\theta}_2^{MLE}$ individually? I can solve the system of equations on Mathematica, but I am unable to do so by hand.
A direction:
I could take $\frac{1}{n} \sum y^+ +\frac{1}{n} \sum y^-$ and by invariance, I get $\left(\frac{\theta_1^2-\theta_2^2}{\theta_1+\theta_2}\right)_{MLE} = (\theta_1 - \theta_2)_{MLE}$. However, I would still need to find a way to isolate the individual parameters.
Since we now have $(\theta_1 - \theta_2)$ we can also find $(\theta_1 + \theta_2)$. The denominator of $\frac{\theta_1^2}{\theta_1 + \theta_2}$ is telling that we should form $\frac{(\theta_1 + \theta_2)^2}{\theta_1 + \theta_2} \to (\theta_1 + \theta_2)$. Then we can solve an easier systems of equations that is linear on $\theta_1$ and $\theta_2$.
From here on, I use simpler notation, $S_p = \frac{1}{n}\sum y^+$ and $S_n = \frac{1}{n}\sum y^-$. By invariance, $\sqrt{-S_n}$ is the MLE estimate of $\frac{\theta_2}{\sqrt{\theta_1 + \theta_2}}$ and $\sqrt{S_p} = \frac{\theta_1}{\sqrt{\theta_1+\theta_2}}$, thus $2\sqrt{-S_n S_p}$ is the MLE of $\frac{2\theta_1 \theta_2}{\theta_1 + \theta_2}$. So
\begin{align*}
\frac{\theta_1^2 + 2\theta_1 \theta_2 + \theta_2^2}{\theta_1 + \theta_2} &= \frac{(\theta_1 + \theta_2)^2}{\theta_1 + \theta_2} \\
&= \theta_1 + \theta_2
\end{align*}
Now we can estimate the above by $S_p + S_n + 2\sqrt{-S_n S_p}$. Thus we have:
$\theta_1 - \theta_2 = S_p + S_n$ and $\theta_1 + \theta_2 = S_p + S_n +2\sqrt{-S_n S_p}$. This is now easily solvable for $(\theta_1, \theta_2)$.

Following from the direction at the end of the post above, the solution is now
\begin{align*}
\theta_1 - \theta_2 &= S_p + S_n \\
\theta_1 + \theta_2 &= S_p - S_n + 2\sqrt{-S_n S_p} \\
&\implies 2\theta_1 = 2S_p + 2\sqrt{-S_n S_p} \implies \theta_1^{MLE} = S_p + \sqrt{-S_n S_p} \\
&\implies -2\theta_2 = 2S_n - 2\sqrt{-S_n S_p} \implies \theta_2^{MLE} = -S_n + \sqrt{-S_n S_p}
\end{align*}
