There are two parts.
(1) As defined, $f$ and $g$ are not only continuous, but also of bounded variation and absolutely continuous -- hence, differentiable almost everywhere. We can show that $f$ is Riemann-Stieltjes integrable with respect to $g$ with
$$\tag{1}\int_a^b f(x) \psi(x) \, dx = \int_a^b f \, dg$$
and similarly
$$\tag{2}\int_a^b g(x) \phi(x) \, dx = \int_a^b g \, df$$
If $\phi$ and $\psi$ were continuous then this would be easy, since we could conclude that $f$ and $g$ were differentiable and by an elementary proof
$$\int_a^b f(x) \psi(x) \, dx = \int_a^b f(x) g'(x) \, dx= \int_a^b f \, dg, \\ \int_a^b g(x) \phi(x) \, dx = \int_a^b g(x) f'(x) \, dx= \int_a^b g \, df.$$
However, we are only given that $\phi$ and $\psi$ are integrable. Nevertheless, (1) and (2) are still true. A proof under more general conditions using Riemann-Stieltjes sums is given here.
(2) Given the results in (1) we can invoke integration by parts for Riemann-Stieltjes integrals
$$\int_a^b f \, dg + \int_a^b g \, df = f(b)g(b) - f(a)g(a).$$
This, again, can be proved in a straightforward way using Riemann-Stieltjes sums.
As a telescoping sum, we have
$$f(b)g(b) - f(a)g(a) = \sum_{k=1}^n[f(x_k)g(x_k) - f(x_{k-1})g(x_{k-1})] \\ = \sum_{k=1}^ng(x_k)f(x_k) - \sum_{k=1}^ng(x_{k-1})f(x_{k-1}). $$
Hence, for any Riemann-Stieltjes sum $S(P;g,f) =\sum_{k=1}^n g(e_k)[f(x_k) - f(x_{k-1})$ we have
$$\tag{3}S(P;g,f) - \left[f(b)g(b) - f(a) g(a) \right] \\= \sum_{k=1}^n g(e_k)[f(x_k) - f(x_{k-1}) - \sum_{k=1}^ng(x_k)f(x_k) + \sum_{k=1}^ng(x_{k-1})f(x_{k-1}) \\ = -\left(\sum_{k=1}^nf(x_{k-1})[g(e_k) - g(x_{k-1})] + \sum_{k=1}^nf(x_{k})[g(x_k) - g(e_{k})] \right) $$
Subsituting with $x_k = y_{2k}$, $x_{k-1} = y_{2k-2}$ and $e_k = y_{2k-1}$ in the last line of (3) we obtain
$$\tag{4}S(P;g,f) - \left[f(b)g(b) - f(a) g(a) \right] \\=-\left(\sum_{k=1}^nf(y_{2k-2})[g(y_{2k-1}) - g(y_{2k-2})] + \sum_{k=1}^nf(y_{2k})[g(y_{2k}) - g(y_{2k-1})] \right)$$
Define special "intermediate points" $\xi_{2k-1} = y_{2k-2}$ and $\xi_{2k} = y_{2k}$. Substituting into (4) we obtain
$$\tag{5}S(P;g,f) - \left[f(b)g(b) - f(a) g(a) \right] \\=-\left(\sum_{k=1}^nf(\xi_{2k-1})[g(y_{2k-1}) - g(y_{2k-2})] + \sum_{k=1}^nf(\xi_{2k})[g(y_{2k}) - g(y_{2k-1})] \right)$$
Taking $j = 2k$, equation (5) is equivalent to
$$S(P;g,f) = f(b)g(b) - f(a) g(a) -\sum_{j=1}^{2n} f(\xi_j)[g(y_j) - g(y_{j-1})]$$
The sum on the RHS is a Riemann-Stieltjes sum of $f$ with respect to $g$ over $[a,b]$, say $S(P,f,g)$. We can rewrite this as
$$S(P;g,f) - \left[f(b)g(b) - f(a) g(a) \right] = S(P,f,g) $$.
As the partition is refined, the RHS converges to $\int_a^b f \, dg$. Whence, the LHS converges and
$$\int_a^b g \, df = f(b)g(b) - f(a) g(a) - \int_a^b f\, dg$$