 Solution

At $\sigma =1$, complete the square:For fixed $y$, all factors not involving $x$ enter its conditional normalizing constant. Interchanging the roles of the coordinates givesThese are the full conditional distributions of the standard bivariate normal distribution with correlation $\rho$.

A systematic Gibbs sampler starts at a chosen $(x_0,y_0)$ and, at each sweep, draws independent standard normals $Z_{1,n},Z_{2,n}$ and setsThe second update uses the newly drawn $x_{n+1}$. Both coordinate updates preserve the target by part (a), hence so does the sweep. The Gaussian two-coordinate Gibbs recursion is explicitlyan AR(1) recursion with innovation variance $1-{\rho }^4$. Its stationary variance is one, and its correlation at sweep lag $k$ is ${\rho }^{2k}$, displaying the dependence of the sample.

The printed general-$\sigma$ expression has a separate normalization defect: for $\sigma >0$ its integral is $\sigma$, not one, since the exponent has covariance ${\sigma }^2\left(\begin{array}{cc}1& \rho \\ \rho & 1\end{array}\right)$, whose normalizing denominator is $2\pi {\sigma }^2\sqrt{1-{\rho }^2}$. Retaining that exponent requires ${\sigma }^2$ in the prefactor. The requested $\sigma =1$ target is already normalized, so this defect does not alter the conditional distributions above.