DIGRESSION
This is a good example to showcase that densities should better be defined for the whole real line using indicator functions and not branches. Because, if one looks at the likelihood, one could, at least for a moment, say "hey, this likelihood will be maximized for the value from the sample that is positive and closest to zero -why not take this as the MLE"?
The density for one typical uniform in this case is
$$f \left( x, \theta_x \right) =\frac{1}{2\theta_x}\cdot \mathbf 1\{x_i \in [-\theta_x,\theta_x] \},\qquad \theta_x >0$$
Note that the interval is (and should be) closed, and that we define the parameter as positive because, defining it as belonging to the real line would a) include the value zero which would make the setup meaningless and b) add nothing to the case except heavy dead-burden notation.
The likelihood of an (ordered) sample of size $i=1,...,n_1$ from i.i.d such r.v.s is
$$L(\theta_x \mid \{x_1,...,x_{n_1}\}) = \frac{1}{2^{n_1}\theta_x^{n_1}}\cdot \prod_{i=1}^{n_1}\mathbf 1\{x_i \in [-\theta_x,\theta_x] \}$$
$$=\frac{1}{2^{n_1}\theta_x^{n_1}}\cdot \min_i\left\{\mathbf 1\{x_i \in [-\theta_x,\theta_x]\}\right\}$$
The existence of the indicator function tells us that if we select a $\hat \theta$ such that even one realized sample value $x_i$ will be outside $[-\hat \theta_x,\hat \theta_x]$, the likelihood will equal zero. Now, respecting this constraint, this likelihood is always higher for positive values of the parameter, and it has a singularity at zero so it is "maximized" (tends to plus infinity) as $\theta_x \rightarrow 0^+$. It is then the constraint to choose a $\theta_x$ such that all realized values of the sample are inside $[-\hat \theta_x,\hat \theta_x]$ that guides us to move away from zero the minimum possible (reducing the value of the likelihood as little as possibly permitted by the constraint), and this is the actual reason why we arrive at the estimator $\hat{\theta_x} = \max \{-X_1,X_{n_1} \}$ and the estimate $\hat{\theta_x} = \max \{-x_1,x_{n_1} \}$.
MAIN ISSUE
To arrive at the joint distribution of $V$ and $U$ as defined in the question, we need first to derive the distribution of the ML estimator. The cdf of $\hat \theta_x$, respecting also the relation $X_1 \le X_{n_1}$, is
$$F_{\theta_x}(\hat{\theta_x}) = P(-X_1 \le \hat{\theta_x}, X_{n_1}\le \hat{\theta_x}\mid X_1 \le X_{n_1}) = P(-\hat{\theta_x} \le X_1 \le X_{n_1}\le \hat{\theta_x})$$
Denoting the joint density of $(X_1, X_{n_1})$ by $f_{X_1X_{n_1}}(x_1,x_{n_1})$, to be derived shortly, the density of the MLE therefore is
$$f_{\theta_x}(\hat{\theta_x}) = \frac {d}{d\hat{\theta_x}}F_{\theta_x}(\hat{\theta_x}) = \frac {d}{d\hat{\theta_x}}\int_{-\hat \theta_x}^{\hat \theta_x}\int_{-\hat \theta_x}^{x_{n_1}}f_{X_1X_{n_1}}(x_1,x_{n_1})dx_1dx_{n_1}$$
Applying (carefully) Leibniz's rule we have
$$f_{\theta_x}(\hat{\theta_x}) = \int_{-\theta_x}^{\hat \theta_x}f_{X_1X_{n_1}}(x_1, \hat \theta_x)dx_1 -(-1)\cdot \int_{-\hat \theta_x}^{-\hat \theta_x}f_{X_1X_{n_1}}(x_1,-\hat \theta_x)dx_1 + \\ +\int_{-\hat \theta_x}^{\hat \theta_x}\left(\frac {d}{d\hat{\theta_x}}\int_{-\theta_x}^{x_{n_1}}f_{X_1X_{n_1}}(x_1,x_{n_1})dx_1\right)dx_{n_1}$$
$$= \int_{-\theta_x}^{\hat \theta_x}f_{X_1X_{n_1}}(x_1, \hat \theta_x)dx_1+0-(-1)\cdot \int_{-\hat \theta_x}^{\hat \theta_x}f_{X_1X_{n_1}}(-\hat \theta_x,x_{n_1})dx_{n_1}$$
$$\Rightarrow f_{\theta_x}(\hat{\theta_x}) =\int_{-\theta_x}^{\hat \theta_x}f_{X_1X_{n_1}}(x_1, \hat \theta_x)dx_1+\int_{-\hat \theta_x}^{\hat \theta_x}f_{X_1X_{n_1}}(-\hat \theta_x,x_{n_1})dx_{n_1}$$
The general expression for the joint distribution of two order statistics is
In our case this becomes
$$f_{X_1,X_{n_1}}(x_1,x_{n_1}) = \frac {n_1!}{(n_1-2)!}f_X(x_1)f_X(x_{n_1})\cdot\left[F_X(x_{n_1})-F_X(x_1)\right]^{n_1-2}$$
$$\Rightarrow f_{X_1,X_{n_1}}(x_1,x_{n_1})=n_1(n_1-1)\left(\frac 1{2\theta_x}\right)^2 \left[\frac {x_{n_1}+\theta_x}{2\theta_x} - \frac {x_1+\theta_x}{2\theta_x}\right]^{n_1-2}$$
$$\Rightarrow f_{X_1,X_{n_1}}(x_1,x_{n_1}) = n_1(n_1-1)\left(\frac 1{2\theta_x}\right)^{n_1}(x_{n_1}-x_1)^{n_1-2}$$
Plugging this into the expression for the density of the MLE and performing the simple integration we get
$$f_{\theta_x}(\hat{\theta_x})=n_1\left(\frac 1{2\theta_x}\right)^{n_1}\cdot \left[-(\hat \theta_x-x_1)^{n_1-1} + (x_{n_1}+\hat \theta_x)^{n_1-1}\right]_{-\hat \theta_x}^{\hat \theta_x}$$
$$=n_1\left(\frac 1{2\theta_x}\right)^{n_1}\cdot 2\cdot (2\hat \theta_x)^{n_1-1}$$
$$\Rightarrow f_{\theta_x}(\hat{\theta_x}) = \frac {n_1}{\theta_x^{n_1}}\hat \theta_x^{n_1-1}$$
By the way , this is a "extended-support" Beta distribution, Beta$(\alpha = n_1, \beta =1, min=0, max = \theta_x)$. For the $Y$ r.v.s the expression would be the same, using $\hat \theta_y,\, \theta_y,\, n_2$.
We turn now to the joint density $g(u,v)$ as defined in the question. $U$ and $V$ are defined as the extreme order statistics of just ...two random variables, $\hat \theta_x$ and $\hat \theta_y$. If we want to apply again the theorem above in order to derive the joint density of $U$ and $V$ under the null $H_0$ (which only makes $\theta_x = \theta_y=\theta$), we need the densities of the MLEs to be identical, and for this we need in addition that $n_1=n_2=n$ (and the mystery is solved). Under these assumptions, application of the theorem (the $n$ of the theorem is now set equal to $2$) we get
$$g(u,v) = 2f_{\theta}(u)f_{\theta}(v)$$
$$=2\frac {n}{\theta^{n}}u^{n-1}\frac {n}{\theta^{n}}v^{n-1} = 2n^2u^{n-1}v^{n-1}/\theta^{2n}$$
Working with extrema requires care, but it doesn't have to be difficult. The crucial question, found near the middle of the post, is
... we need to show that $\frac{n+1}{n}{\hat \theta}^{n}$ is unbiased.
Earlier you obtained
$$\hat\theta = \max(-y_1, y_n/2) = \max\{-\min\{y_i\}, \max\{y_i\}/2\}.$$
Although that looks messy, the calculations become elementary when you consider the cumulative distribution function $F$. To get started with this, note that $0\le \hat\theta \le \theta$. Let $t$ be a number in this range. By definition,
$$\eqalign{
F(t) &= \Pr(\hat\theta\le t) \\&= \Pr(-y_1 \lt t\text{ and }y_n/2 \le t) \\
&= \Pr(-t \le y_1 \le y_2 \cdots \le y_n \le 2t).
}$$
This is the chance that all $n$ values lie between $-t$ and $2t$. Those values bound an interval of length $3t$. Because the distribution is uniform, the probability that any specific $y_i$ lies in this interval is proportional to its length:
$$\Pr(y_i \in [-t, 2t]) = \frac{3t}{3\theta} = \frac{t}{\theta}.$$
Because the $y_i$ are independent, these probabilities multiply, giving
$$F(t) = \left(\frac{t}{\theta}\right)^n.$$
The expectation can immediately be found by integrating the survival function $1-F$ over the interval of possible values for $\hat\theta$, $[0, \theta]$, using $y=t/\theta$ for the variable:
$$\mathbb{E}(\hat\theta) = \int_0^\theta \left(1 - \left(\frac{t}{\theta}\right)^n\right)dt = \int_0^1 (1-y^n)\theta dy = \frac{n}{n+1}\theta.$$
(This formula for the expectation is derived from the usual integral via integration by parts. Details are given at the end of https://stats.stackexchange.com/a/105464.)
Rescaling by $(n+1)/n$ gives
$$\mathbb{E}\left(\frac{n+1}{n}\,{\hat \theta}\right) = \theta,$$
QED.
Best Answer
Since$$(y_1,\ldots,y_r)\sim\frac{n!\theta^{-r}}{(n-r)!}e^{-\frac{1}{\theta}[\sum_{i=1}^{r}y_i+(n-r)y_r]}\mathbb{I}_{y_\le y_2\le \ldots \le y_r}$$you have the joint pdf of $(y_1,\ldots,y_r)$. From there, you can deduce the pdf of $$s_r=\sum_{i=1}^{r}y_i+(n-r)y_r\,.$$Indeed, because the Jacobian of the transform is constant,\begin{align*}f_s(y_1,\ldots,y_{r-1},s_r) &\propto f_Y\left(y_1,\ldots,\left\{s_r-\sum_{i=1}^{r-1}y_i\right\}\Big/(n-r+1)\right) \\&\propto \theta^{-r} \exp\{-s_r/\theta\}\mathbb{I}_{y_\le y_2\le \ldots \le\left\{s_r-\sum_{i=1}^{r-1}y_i\right\}/(n-r+1)}\end{align*}implies by integration in $y_1,\ldots,y_{r-1}$ that$$f_s(s_r)\propto\theta^{-r} \exp\{-s_r/\theta\}s_r^{r-1}$$ Indeed, \begin{align*} f_s(s_r)&=\int\cdots\int f_s(y_1,\ldots,y_{r-1},s_r)\text{d}y_1\cdots\text{d}y_{r-1}\\ &= \theta^{-r} \exp\{-s_r/\theta\}\int\cdots\int \mathbb{I}_{y_\le y_2\le \ldots \le\left\{s_r-\sum_{i=1}^{r-1}y_i\right\}/(n-r+1)}\text{d}y_1\cdots\text{d}y_{r-1} \end{align*} leads to constraint $y_{r-1}$ by $y_{r-2}\le y_{r-1}$ and by $$y_{r-1}\le \left\{s_r-\sum_{i=1}^{r-1}y_i\right\}/(n-r+1)=\left\{s_r-\sum_{i=1}^{r-2}y_i\right\}/(n-r+1)-\frac{y_{r-1}}{n-r+1}$$ which simplifies into $$y_{r-1}\le \left\{s_r-\sum_{i=1}^{r-2}y_i\right\}/(n-r+2)$$ If one starts integrating in $y_{r-1}$, the most inner integral is \begin{align*}\int_{y_{r-2}}^{\{s_r-\sum_{i=1}^{r-2}y_i\}/(n-r+2)}\text{d}y_{r-1}&=\left\{s_r-\sum_{i=1}^{r-2}y_i\right\}/(n-r+2)-y_{r-2}\\ &=\left\{s_r-\sum_{i=1}^{r-3}y_i\right\}/(n-r+2)-\frac{(n-r+1)y_{r-2}}{n-r+2} \end{align*} and from there one can proceed by recursion.
Hence$$s_r\sim\mathcal{G}a(r,1/\theta)$$
Here is an R simulation to show the fit:
obtained as follows