Gibbs Sampler, Bivariate Normal, Subchain In Example 7.1 of "Introducing Monte Carlo Methods with R", the authors write
$(X,Y)\sim N\Bigg((0,0),\begin{pmatrix}1 &\rho \\ \rho & 1\end{pmatrix}\Bigg)$
Then,


*

*Given $x_t$,

*$Y_{t+1}\mid x_t \sim N(\rho x_t, 1-\rho^2)$

*$X_{t+1}\mid y_{t+1} \sim N(\rho y_{t+1}, 1-\rho^2)$


Finally, they say: The subchain $(X_t)_t$ satisfies 
$$ X_{t+1} \mid X_t = x_t\sim N(\rho ^2 x_t,1-\rho^4) $$.
I'm at a loss for how to show this.  I did
$$
\begin{align}
f(x_{t+1}\mid x_{t}) &= \int_\mathbb{R} f(x_{t+1}, y_{t+1} \mid x_t)dy_{t+1}\\
&=\frac{1}{2\pi(1-\rho^2)}\int_{\mathbb{R}}exp\Big[ \frac{-1}{2(1-\rho^2)}\big((y_{t+1} - \rho x_t)^2+ (x_{t+1}-\rho y_{t+1})^2\big) \Big] dy_{t+1}
\end{align}
$$
where I have assumed  $X_{t+1}$ is conditionally independent of $X_t$ given $Y_{t+1}$.
But I just don't know where to go from there.  Any help would be appreciated.
 A: Sorry if this is unclear in our book, but when
\begin{align*}
Y_{t+1}\mid X_t=x_t &\sim \mathrm{N}(\rho x_t, 1-\rho^2)\\
X_{t+1}\mid Y_{t-1}=y_{t+1},X_t=x_t &\sim \mathrm{N}(\rho y_{t+1}, 1-\rho^2)
\end{align*}
one gets
$$X_{t+1}|X_t=x_t \sim \mathrm{N}(\rho \times \rho x_t,1-\rho^2+\rho^2(1-\rho^2))$$
as
$$\mathbb{E}[X_{t+1}|X_t]=\underbrace{\mathbb{E}[\mathbb{E}[X_{t+1}|Y_{t+1}]|X_t]}_{\mathbb{E}[\rho Y_{t+1}|X_t]=\rho^2 X_t}$$
and
$$\text{var}(X_{t+1}|X_t)=\underbrace{\mathbb{E}[\text{var}(X_{t+1}|Y_{t+1})|X_t)}_{\mathbb{E}[1-\rho^2|X_t]=1-\rho^2}+\underbrace{\text{var}(\mathbb{E}[X_{t+1}|Y_{t+1}]|X_t)}_{\rho^2\text{var}(Y_{t+1}|X_t)=\rho^2(1-\rho^2)}$$
If one wants to solve by the integral in the question,
\begin{align*}
&\exp\left\{\frac{-1}{2(1-\rho^2)}\big((y - \rho x_t)^2+ (x_{t+1}-\rho y)^2\big)\right\}\\
&\quad = \exp\left\{\frac{-1}{2(1-\rho^2)}\big(y^2(1+\rho^2)-2y(\rho x_t+\rho x_{t+1})+\rho^2 x_t^2+x_{t+1}^2\big)\right\}\\
&\quad = \exp\left\{\frac{-(1+\rho^2)}{2(1-\rho^2)}\big(y^2-2\rho (x_t+x_{t+1})y/(1+\rho^2)+\rho^2(x_t+x_{t+1})^2/(1+\rho^2)^2\big)\right\}\\
&\quad\times\exp\left\{\frac{-1}{2(1-\rho^2)}\big(\rho^2 x_t^2+x_{t+1}^2-\rho^2(x_t+x_{t+1})^2/(1+\rho^2)\big)\right\}\\
\end{align*}
which makes the first exponential term a perfect normal density in $y$ and the second term can be simplified into
\begin{align*}
&\exp\left\{\frac{-1}{2(1-\rho^2)}\big(\rho^2 x_t^2+x_{t+1}^2-\rho^2(x_t+x_{t+1})^2/(1+\rho^2)\big)\right\}\\
&\quad=\exp\left\{\frac{-1}{2(1-\rho^2)(1+\rho^2)}\big(
(1+\rho^2)[\rho^2 x_t^2+x_{t+1}^2]-\rho^2(x_t+x_{t+1})^2\big)\right\}\\
&\quad=\exp\left\{\frac{-1}{2(1-\rho^4}\big(
x_{t+1}^2[1+\rho^2-\rho^2]-2x_tx_{t+1} \rho^2+x_t^2\rho^2[1+\rho^2-\rho^2]\big)\right\}\\&\quad=\exp\left\{\frac{-1}{2(1-\rho^4)}\big(
x_{t+1}-\rho^2 x_t\big)^2\right\}\\\end{align*}
which leads to the same result (!)
A: Perhaps this view helps:
\begin{align}
Y_{t + 1} &= \rho x_t + \sqrt{1 - \rho^2} N_t \\
X_{t + 1} 
&= \rho Y_{t + 1} + \sqrt{1 - \rho^2} N_{t + 1} \\
&= \rho^2 x_t + \rho \sqrt{1 - \rho^2} N_t + \sqrt{1 - \rho^2} N_{t + 1},
\end{align}
where $N_t$ and $N_{t + 1}$ are standard normal distributed. It is then easy to see that 
$$\mathbb{E}[X_{t + 1} | x_t] = \rho^2 x_t + 0 + 0$$
and
$$\mathbb{V}[X_{t + 1} \mid x_t] = \rho^2(1 - \rho^2) + (1 - \rho^2) = (1 + \rho^2)(1 - \rho^2) = 1 - \rho^4,$$
where we have used that the variance of the sum is the sum of the variances for uncorrelated variables and the third binomial formula.
