Gaussian conjugacy for an initialized AR(2) regression (source code)

= Gaussian conjugacy for an initialized AR(2) regression
{c}
{title2=$\Sigma=\Lambda^{-1},\quad\mu=\Lambda^{-1}h$}

Condition on initial observations $x_0,x_1$ and model
$$
x_{t+2}=a x_{t+1}+b x_t+\varepsilon_t,\qquad
\varepsilon_t\overset{\rm iid}\sim N(0,\sigma^2),
$$
with known $\sigma^2>0$ and <independent> standard-<normal> <prior distributions> for $a,b$. Put $v_t=(x_{t+1},x_t)^{\mathsf T}$ and $\beta=(a,b)^{\mathsf T}$. <Gaussian conjugacy for a normal linear model> gives posterior <precision matrix> and <posterior mean>
$$
\Lambda=I+\sigma^{-2}\sum_t v_tv_t^{\mathsf T},\qquad
h=\sigma^{-2}\sum_t v_t x_{t+2},\qquad
\beta\mid x\sim N(\Lambda^{-1}h,\Lambda^{-1}).
$$
The <prior distribution> makes $\Lambda$ <positive-definite> even with a deficient <design matrix>. The unscaled sums occur specifically when $\sigma^2=1$.

If $\Lambda=\left(\begin{smallmatrix}A&C\\C&D\end{smallmatrix}\right)$ and $h=(r,s)^{\mathsf T}$, the two <full conditional distributions> have means $(r-Cb)/A$ and $(s-Ca)/D$ and <variances> $1/A$ and $1/D$. These formulas use the conditional <autoregressive model> of order two <likelihood function>. Giving the initial pair its stationary law, or imposing a stationary-parameter restriction, changes that <posterior distribution>.