Parameter estimation in Dynamic Linear Models

I am currently developing a DLM of the following form

$$\underset{k \times 1} {y_t} = \underset{k \times n}A \underset{n \times 1}{\theta_t} + \epsilon_t$$ $$\theta_t = \mu + \underset{n \times n}B\theta_{t-1} + \eta_t$$

where $$\epsilon_t$$ and $$\eta_t$$ are uncorrelated normally distributed errors.

I know how to get filtered and smoothed estimates of the hidden parameters $$(\theta_t)_t$$ when matrices $$A$$ and $$B$$ are known. However, assume matrix $$B$$ has unknown parameters in it, how am I supposed to estimate them using MLE? Without knowing $$B$$, I cannot have any filtered or smoothed estimates of $$(\theta_t)_t$$ yet, so I cannot figure out how to use just the observed data $$(y_t)_t$$ to have an estimation of those parameters.

Given some set of parameter values $$\boldsymbol{\beta}$$ you can construct $$B$$ and run the Kalman filter. It turns out the Kalman filter can be used to recursively calculate the marginal likelihood (see this paper by Tusell). The general idea is to factorize the likelihood function into a product of conditional likelihoods, where we sequentially evaluate the likelihood of each observation given all previous observations.
In R you could write a function, taking parameter values $$\boldsymbol{\beta}$$ as inputs, which
1. builds matrix $$B$$;