For the multivariate model \(A_{0,t} y_t = u_t\) with \(u_t \sim N(0, \Omega_t)\)
the function produces a draw of the lower triangular part of \(A_{0,t}\) similar as in
Primiceri (2005), i.e., using $$y_t = Z_t \psi_t + u_t,$$
where
$$Z_{t} = \begin{bmatrix} 0 & \dotsm & \dotsm & 0 \\ -y_{1, t} & 0 & \dotsm & 0 \\ 0 & -y_{[1,2], t} & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dotsm & 0 & -y_{[1,...,K-1], t} \end{bmatrix}$$
and \(y_{[1,...,K-1], t}\) denotes the first to \((K-1)\)th elements of the vector \(y_t\).
The algorithm of Chan and Jeliazkov (2009) is used to obtain time varying coefficients.