Proof: $Y$ stochastically dominates $X$ implies $E[\phi(Y)]\geq E[\phi(X)]$ for increasing $\phi$ 
Suppose $X$ and $Y$ are real random variables with CDF $F$ and $G$ such that $F(x)\geq G(x)$ (i.e. $Y$ exhibits (first-order) stochastic dominance over $X$). Then, for all increasing function $\phi$, we have $E[\phi(Y)]\geq E[\phi(X)]$. I would like to find a rigorous proof of this. Can you please provide one or point me to a reference?


For the special case in which both $F$ and $G$ are continuous, one can consider the integral
$$
\int \phi(y)dG(y)
$$
and employ the substitution $y\to\varphi(x)$ where $\varphi(x)\equiv G^{-1}[F(x)]$ to infer that
$$
E[\phi(Y)]=\int\phi(y)dG(y){\color{red}=}\int\phi[\varphi(x)]dF(x)\geq\int\phi(x)dF(x)=E[\phi(X)].
$$
The inequality above is because
$$
F(x)\geq G(x)\implies\varphi(x)=G^{-1}[F(x)]\geq G^{-1}[G(x)]{\color{blue}=}x;\quad \phi\text{ is increasing.}
$$
I've made some attempt for when $F$ and $G$ are general by considering $\varphi(x)=G^*[F(x)]$ where $G^*$ now is the (left-continuous) quantile function of $Y$. I've run into at least 2 problems: (i) I don't think the equality ${\color{red}=}$ above continues to hold; (ii) $G^*[G(x)]\leq x$ so ${\color{blue}=}$ above is turned into $\leq$ which breaks the inequality chain.
A rigorous proof can also follow from a different direction. But this different direction relies on a construction I don't know how to do.
 A: The idea is to construct a measure $\mu$on an product space $\Bbb{R}\times \Bbb{R}$ such that $\mu((-\infty,a),\Bbb{R}) = F(a)$ and $\mu(\Bbb{R},(-\infty,b)) = G(b)$ such that $\mu((x,y):  x>y ) = 0$. 
Once you have this, then $Y = X + Z$ with $Z\geq 0$ and noting that $\phi$ is increasing 
$$\Bbb{E}[\phi(Y)] = \Bbb{E}[\phi(X+Z)] \geq \Bbb{E}[\phi(X)] $$
claim: such a $\mu$ exists. 
Indeed, consider 
$X_n = \sum_{i= 0}^{n-1} \delta_{a^n_i} \frac{1}{n}$ where $a^n_i = \inf\big\{a, F(a)\geq\frac{i}{N}\big\}$ and 
$Y_n = \sum_{i= 0}^{n-1} \delta_{b^n_i} \frac{1}{n}$ where $b^n_i = \inf\big\{b, G(b)\geq\frac{i}{n}\big\}$. 
Note that $b^n_i \geq a^n_i$. remark also that one might have $a^n_i = a^n_{i+1}$ 
Construct $\mu^n$ the following way $X^n = a_i^n \Rightarrow Y^n = b_i^n$ (if $X^n = a^n_i$ then $Y^n = b^n_i$). This way $Y^n \geq X^n$. 
Now note for every $x$ there is a $k^n(x): \frac{k^n(x)}{n} < F(x) \leq \frac{k^n(x)+1}{n}$ 
$$F^n(x) = \mu^n(Y^n \leq x) = \frac{\sup\{i, a^n_i \leq x\}}{n} = \frac{k^n(x)}{n}$$
So $X^n\Rightarrow X$ -$X^n$ converges in distribution to $X$ - once the cumulative distributions converge (for every (continuity) point of $F$).
Simmilarly for every $y$ there is a $\tilde{k}^n(y): \frac{\tilde{k}^n(y)}{n} \leq G(y) < \frac{\tilde{k}^n(y)+1}{n}$ 
$$G^n(y) =\mu^n(Y^n \leq y) = \frac{\sup\{i, b^n_i \leq y\}}{n} = \frac{\tilde{k}^n(y)}{n}$$
So $Y^n\Rightarrow Y$ -$Y^n$ converges in distribution to $Y$ - once the cumulative distributions converge (for every (continuity) point of $G$).
Note further that $\mu^n$ is tight. Indeed let $\alpha$ be such that $F(\alpha)< \epsilon$ and let $\beta$ be such that $G(\beta)> 1 - \epsilon$.
Let $N$ be such that for $n > N$
$$ |F(\alpha) - F^n(\alpha)|< \epsilon\\
|G(\beta) - G^(\beta)| \leq \epsilon $$
Then $$\mu^n((\alpha, \beta)^2) =  1 - \mu^n (X^n \leq \alpha, Y^n \geq \beta) \geq 1 - \mu^n (X^n \leq \alpha) - \mu^n (Y^n \geq \beta) \geq 1 - 4 \epsilon $$
Now, since the sequence of measures is tight (by Prohorov) it admits a limit $\mu$ for some subsequence.
Moreover (by Portemanteau) $\liminf \mu^n(A) \geq \mu(A)$ for every open $A$ 
Take $A = \{(x,y): x>y\}$ and note (to conclude) that $\mu^n(A) = 0$
A: Here is an easier proof.
Since increasing functions are differentiable almost surely, we may consider all $\phi$ that is both differentiable and increasing without loss of generality.
Define
$$
H(x)=G(x)-F(x)\leq 0,\quad H(0)=0,\quad \lim_{x\rightarrow \infty}H(x)=0.
$$ 
Apply the Integration by Parts for the generalized distribution function
$$
\int \phi(x)dH(x)=[\phi(x)H(x)]^\infty_{0}-\int \phi'(x)H(x-)dx
$$
where
$$
H(x-)=\lim_{t\uparrow x}H(t)\leq 0.
$$
The first term vanishes by construction. The second term is evidently positive. It follows that
$$
\int \phi(x)dG(x)-\int \phi(x)dF(x)=\int \phi(x)dH(x)\geq 0.
$$
This completes the proof.
