Assumes a normal-inverse gamma model,
$$x_t = \phi_0 + \phi_1 x_{t-1} + \dots + \phi_p x_{t-p} + \sigma z_t ,$$
where \(z_t\) is standard Gaussian noise.
With \(\Phi\) being the (p+1)-dimensional vector of the \(\phi\)s,
the priors are
\(\Phi \mid \sigma \sim N(0, \sigma^2 V_0)\) and
\(\sigma^2 \sim IG(a,b)\), where \(V_0 = \gamma^2 I\).
Defaults are given for the hyperparameters, but the user
may choose \((a,b)\) as (prior_sig_a, prior_sig_b)
and \(\gamma^2\) as prior_var_phi.
The algorithm is efficient and converges quickly.
Further details can be found in Chapter 6 of the 5th edition of the
Springer text.