Help to understand martingale example from Billingsley I am studying Billingsley's section 35 about martingales and I have some difficulties understanding one of the examples.
This is the example 35.4.
Suppose we have a measurable space $(\Omega,\mathcal F)$ and let $Q$ and $P$ be probability measures on the latter. Let $\mathcal F_n=\sigma (Y_1,\cdots,Y_n)$ be the sigma algebra generated by the random variables $Y_1,\cdots,Y_n$.
Suppose that under  $P$ the distribution of the random vector $(Y_1,\cdots,Y_n)$ has density $p_n(y_1\cdots,y_n)$ with respect to the n-dimensional Lebesgue measure, and under $Q$ it has density $q_n(y_1\cdots,y_n)$. (To avoid technicalities assume that $p_n,q_n$ are everywhere positive.
Then the Radon-Nikodym derivative of $Q$ with respect to $P$ on $\mathcal F_n$ is 
$$X_n=\frac{q_n(y_1\cdots,y_n)}{p_n(y_1\cdots,y_n)}$$
So far so good, but the author then concludes that $X_n$ is a martigale under $P$ w.r.t. the filtration $\mathcal F_n$.
How can I see this last point, 
why is that the following holds?
$$\int_F X_n dP=\int_F X_{n+1} dP, \forall F\in\mathcal F_n$$
Thanks in advance.
EDIT:
I may have  found a possible solution, but I am not sure about the validity of it.
I can assume that there is a random variable $\Phi\in\mathcal M(\mathcal F)$ such that
$$Q(F)=\int_F \Phi dP,\forall F\in\mathcal F_n$$
(One problem might be that then $Q$ is not a probability measure on $\mathcal F$).
Then I can interpret the Radon-Nikodym derivative $X_n$ as the condition expectation $E[\Phi|\mathcal F_n]$.
Finally I have that 
 $$E[X_{n+1}|\mathcal F_n]=E[E[\Phi|\mathcal F_{n+1}]|\mathcal F_n]=E[\Phi|\mathcal F_n]=X_n$$
$$E[X_{n+1}]=E[E[X_{n+1}|\mathcal F_n]]=E[X_n]$$
 A: It's a little strange to say that $X_n=\frac{q_n(y_1\cdots,y_n)}{p_n(y_1\cdots,y_n)}$ is the "the Radon-Nikodym derivative of $Q$ with respect to $P$ on $\mathcal{F}_n$."
In general, a Radon-Nikodym derivative $\frac{dQ}{dP}$ is a random variable defined on $(\Omega,\mathcal{F})$, whereas $X_n=\frac{q_n(y_1\cdots,y_n)}{p_n(y_1\cdots,y_n)}$ is defined on $\mathbb{R}^n$.
I will assume a specific $\Omega$ so that $X_n$ is actually $\frac{dQ}{dP}$. Take $\Omega = \mathbb{R}^{\infty}$, with the Borel $\sigma$-algebra given by the product topology, and $Y_n$ is the $n$-th coordinate map. So $\mathcal F_n=\sigma (Y_1,\cdots,Y_n)$ consists of cylinder sets given by the first $n$ coordinates. In this case, $X_n = \frac{q_n(y_1\cdots,y_n)}{p_n(y_1\cdots,y_n)} = \frac{dQ}{dP}$.
The answer to your question now simply follows from the factorization
$$
q_{n+1}(y_1\cdots,y_{n+1}) = q_{n}(y_1\cdots,y_{n}) q_{n+1|n}(y_{n+1}|y_1\cdots,y_{n})
$$
where $q_{n+1|n}(y_{n+1}|y_1\cdots,y_{n})$ is the conditional density of $y_{n+1}$ given $y_1\cdots,y_{n}$.
To be more explicit, a general element of $\mathcal{F}_n$ takes the form
$F_n = B_n \times \mathbb{R} \times \mathbb{R} \cdots$, where $B_n$ is a Borel set in  $\mathbb{R}^n$. So
\begin{align*}
\int_{F_n} X_{n+1} dP &= \int_{B_n \times \mathbb{R}} \frac{q_{n+1}}{p_{n+1}} p_{n+1} dy_1 \cdots dy_{n+1}\\
&= \int_{B_n \times \mathbb{R}} q_{n+1} dy_1 \cdots dy_{n+1}\\
&= \int_{B_n \times \mathbb{R}} q_{n} q_{n+1 | n} dy_1 \cdots dy_{n+1}\\
&= \int_{B_n} (\int_{\mathbb{R}} q_{n+1 | n} ) q_{n}  dy_1 \cdots dy_{n}\\
&= \int_{B_n} q_{n}  dy_1 \cdots dy_{n}\\
&= \int_{B_n} \frac{q_{n}}{p_n} p_n dy_1 \cdots dy_{n}\\
&= \int_{F_n} X_{n} dP.\\
\end{align*}
(This really just a martingale of the form $E[X|\mathcal{F}_n]$ where $X$ is an $L^1$ random variable and $\{ \mathcal{F}_n \}$ is a filtration.)
