Solution (source code)

= Solution

For fixed parameters $(a,b)$, every new error is independent of the previous observations. Put $\varphi_0(z)=(2\pi)^{-1/2}e^{-z^2/2}$. Since $X_0=X_1=0$, the first three conditional densities are
$$
\boxed{f_{X_2\mid X_1}(x_2\mid0)=\varphi_0(x_2),\qquad
f_{X_3\mid X_2,X_1}(x_3\mid x_2,0)=\varphi_0(x_3-ax_2),}
$$
and
$$
\boxed{f_{X_4\mid X_3,X_2,X_1}(x_4\mid x_3,x_2,0)=\varphi_0(x_4-ax_3-bx_2).}
$$
In general,
$$
\boxed{X_{t+2}\mid(X_{t+1},\ldots,X_1),a,b\sim N(aX_{t+1}+bX_t,1).}
$$
Its density is $\varphi_0(x_{t+2}-ax_{t+1}-bx_t)$. The recursion is interpreted from $t=0$, as required to define $X_2$ from the two given initial values. These are parameter-conditional sampling distributions; integrating over the prior would instead give predictive mixtures.