Solution (source code)

= Solution

First integrate out $\eta_i$: conditional on $\xi_i$, $y_i\sim N(\alpha+\beta\xi_i,s^2)$ with $s^2=\sigma^2+\sigma_y^2$. Since $(x_i,y_i)$ is an affine transformation of independent normal variables, it is <multivariate normal distribution>[bivariate normal] with mean
$$
m=\binom{\mu}{\alpha+\beta\mu}
$$
and covariance
$$
\Sigma=
\begin{pmatrix}
A&B\\B&C
\end{pmatrix},
\quad
A=\tau^2+\sigma_x^2,
\quad B=\beta\tau^2,
\quad C=\beta^2\tau^2+s^2.
$$
Its determinant simplifies to
$$
D=AC-B^2=(\tau^2+\sigma_x^2)s^2+\beta^2\tau^2\sigma_x^2.
$$
For $r_i=(x_i-\mu,y_i-\alpha-\beta\mu)^T$, the observed-data likelihood is therefore
$$
L=(2\pi)^{-N}D^{-N/2}
\exp\!\left[-\frac12\sum_{i=1}^N
r_i^T\Sigma^{-1}r_i\right],
\qquad
\Sigma^{-1}=\frac1D\begin{pmatrix}C&-B\\-B&A\end{pmatrix}.
$$