Expectation of product of compositions Suppose $(\Omega,\mathcal{A},P)$ be a probability space, and let $X,Y$ be two random variables on $\Omega$. Also, consider the (measurable) indicator function $\chi_{(-\infty,x]}$ given any $x\in \mathbb{R}$. I know that
$$\int_\Omega (\chi_{(-\infty,x]}\circ X)(w) P(dw)=\int_{\mathbb{R}}\chi_{(-\infty,x]}(w')P_X(dw')=P_X((-\infty,x])=P(\{X\leq x\})$$
where $P_X$ is the probability distribution of $X$.

Now I'm struggling to show that

$$E[(\chi_{(-\infty,x]}\circ X)(\chi_{(-\infty,y]}\circ Y)]=\int_\Omega [(\chi_{(-\infty,x]}\circ X)(\chi_{(-\infty,y]}\circ Y)](w)P(dw)=P(\{X\leq x,Y\leq y\})$$

It looks simple, but I cannot progress. Can you help me?

*I noticed the identity $(\chi_{(-\infty,x]}\circ X)(\chi_{(-\infty,y]}\circ Y)=\chi_{(-\infty,x]\times (-\infty,y]}(X,Y)$. So I could use

$$\int_\Omega \chi_{(-\infty,x]\times (-\infty,y]}(X(w),Y(w))P(dw)=\int_{\mathbb{R}^2} \chi_{(-\infty,x]\times (-\infty,y]}(w_1,w_2)P_{XY}(dw_1\otimes d2_2)=P(\{X\leq x,Y\leq y\})$$

What do you think?
I'm not sure if the transformation theorem works for the bivariate case as I did...
Let $(\Omega, \mathcal F, \mathbb P)$ be probability space, let $X: \Omega \to \mathbb R^n$ be random variable. (We consider $\mathbb R^n$ with borel $\mathcal B(\mathbb R^n)$ sigma field.) Moreover, let $f:\mathbb R^n \to \mathbb R$ be arbitrary function and $\mu_X$ the distribution of $X$. Then we have:

$$\mathbb E[f\circ X] = \int_\Omega (f \circ X)(\omega)d\mathbb P(\omega) = \int_{\mathbb R^n} f(s)d\mu_X(s)$$

Proof:

1) case: Let $f = \chi_A$, where $A \in \mathcal B(\mathbb R^n)$. 

Then $$\mathbb E[ f \circ X] = 1 \cdot \mathbb P(X \in A) + 0 \cdot \mathbb P(X \notin A) = \mu_X(A) = \int_{\mathbb R^n}f(s)d\mu_X(s)$$

We only need that one above in your problem, but let's just at least point out, how to prove the rest:

2) case: $f$ is a simple function, that is there are $A_1,...,A_m \in \mathcal B(\mathbb R^n)$ and $a_1,...,a_m \in \mathbb R$ such that $f = \sum_{i=1}^m a_i\chi_{A_i}$. Then we can use 1) case + linearity of integrals and we're done.

3) case: $f$ is any non-negative borel function. Then there exists an inscreasing sequence of non-negative simple functions $f_n$ converging to $f$. Use Lebesgue-Levy monotone convergence theorem and you're done.

4) case: any borel $f$. Then $f = f^+ - f^-$ and use your result from case 3) to $f^+$ and $f^-$ (note that at least one $\mathbb E[f^+ \circ X], \mathbb E[f^- \circ X]$ is finite).

Now in your problem, as our random variable $X$ take vector $(X,Y): \Omega \to \mathbb R^2$ and as $f$ take $\chi_A$, where $A = (-\infty, x] \times (-\infty,y]$.

Use case 1) here: $\mathbb E[ f \circ (X,Y) ] = \mu_{(X,Y)}(A) = \mathbb P( (X,Y) \in A) = \mathbb P(X \le x, Y \le y)$
The integrand is $1$ on $\{X\leq x, Y\leq y\}$ (which is measurable since $X,Y$ are random variables) and zero elsewhere, so the result immediately follows by the definition of the integral of a simple function.