Thanks @Charles_Driver !

Thinking about your suggestion the only way I see to do it is something like that:

y_t =\alpha z_t + \theta_t +\epsilon, \epsilon \sim N(0, \sigma^2)

\theta_t =\beta x_t + \delta \theta_{t-1} + \omega, \omega \sim N(0, \rho^2)

becomes

y_t - \alpha z_t =\begin{bmatrix} 1 & 0 \end{bmatrix}
\begin{bmatrix} \theta_t \\ 1 \end{bmatrix} +
\epsilon, \epsilon \sim N(0, \sigma^2)

\begin{bmatrix}
\theta_t \\ 1
\end{bmatrix} =
\begin{bmatrix}
\delta & \beta x_t \\
0 & 1
\end{bmatrix}
\begin{bmatrix}
\theta_{t-1} \\ 1
\end{bmatrix} +
\begin{bmatrix}
\omega \\ 0
\end{bmatrix},
\omega \sim N(0, \rho^2)

with m_0 = \begin{bmatrix} \theta_0 \\ 1 \end{bmatrix}, C_0 = \begin{bmatrix} C_0 & 0 \\ 0 & 0 \end{bmatrix}.

Somehow, setting variances to 0 makes me uncomfortable. Never sure if it creates computational problems.

Please, if you have other, most likely smarter, options fell free to share.

In any case it would be nice to have the ability to specify the whole model directly.

Cheers,

F