Showing that double integral of indicator random variable is the product of the random variables Let $X$ and $Y$ be non-negative random variables with an arbitrary joint probability distribution
function. Let $I(x, y) = 1$ if $X > x, Y > y$ and $0$ otherwise.
Show that $\int_0^\infty \int_0^\infty I(x, y) dxdy = XY$.
I'm stuck, but if I assumed what is to be proven, then can it be concluded that both $X$ and $Y$ are constant? Is part of the problem to first show that $X$ and $Y$ are constant?

Because notation is the crux of the matter, let's be rigorous and explicit about it.
To say that $(X,Y)$ has a joint distribution means $(X,Y)$ is a (measurable) function from a measure space $(\Omega, \mathfrak F)$ to $\mathbb R^2.$  The definition of the integrand $I$ involves four quantities: $x,$ $y,$ $X,$ and $Y,$ so let's make this explicit by writing the integral in the question as
$$\int_0^\infty\int_0^\infty I(x,y)\,\mathrm d x \mathrm d y = \int_0^\infty\int_0^\infty I(x,y,X,Y)\,\mathrm d x \mathrm d y.$$
What does this mean?  The integrand $I$ is a (measurable) function of four real variables $z_1,z_2,z_3,z_4$ defined by the (usual) indicator function $\mathscr I$
$$I(z_1,z_2,z_3,z_4) = \mathscr{I}(z_3\gt z_1,z_4\gt z_2) =  \left\{\begin{aligned} 1 & \text{ if } z_3 \gt z_1\text{ and } z_4 \gt z_2 \\ 0 & \text{ otherwise.}\end{aligned}\right.$$
Now, integration of a (measurable) function $g:\mathbb R^n\to\mathbb R$ over any (measurable) region $\mathscr R \subset \mathbb R^n$ is defined in terms of an integral over the entire space as
$$\iint_{\mathscr R} g(\mathbf x)\,\mathbf{\mathrm{d}}\mathbf x = \iint_{\mathbb R^n} \mathscr{I}(\mathbf x \in \mathscr R) g(\mathbf x)\,\mathbf{\mathrm{d}}\mathbf x.$$
Let's pause to note a basic property of the indicator function: multiplication corresponds to intersection in the sense that for any two sets $\mathscr R$ and $\mathscr S,$
$$\mathscr{I}(x \in \mathscr R)\mathscr{I}(x \in \mathscr S) = \mathscr{I}(x\in \mathscr R \cap \mathscr S).$$
Use this to define a function $f:\mathbb R^2\to \mathbb R$ as follows:
$$\begin{aligned}
f(u,v) &= \int_0^\infty\int_0^\infty I(x,y, u,v)\,\mathrm d x \mathrm d y\\
&=\iint_{\mathbf R^2}\mathscr{I}(0 \le x, 0 \le y)\mathscr{I}(x \lt u, y \lt v) \,\mathrm d x \mathrm d y \\
&=\iint_{\mathbf R^2}\mathscr{I}(0 \le x \lt u, 0 \le y \lt v) \,\mathrm d x \mathrm d y \\
&= \int_0^{u}\int_0^{v}\mathrm d x \mathrm d y = uv 
\end{aligned}$$
(assuming $u$ and $v$ are both non-negative; the result is $0$ otherwise).
That's the area of a rectangle of width $v$ and height $u,$ equal to $uv.$  (Don't read too much into these calculations.  This is all elementary set theory and geometry at this point, having nothing to do with random variables or even with any theory of integration.  The only conception of integration you need is that the integral of the constant $c$ over a rectangle is $c$ times the area of the rectangle: this is the starting point of all theories of integration.)
Composing $f$ with $(X,Y)$ when both $X$ and $Y$ always have non-negative values gives the (necessarily measurable) function
$$f\circ (X,Y):(\Omega,\mathfrak F) \to \mathbb R$$
where, by definition, the value of the composition at any $\omega\in\Omega$ is
$$(f\circ (X,Y))(\omega) = f(X(\omega), Y(\omega)) = X(\omega)Y(\omega).$$
That's what "$XY$" means, QED.
