Kalman/HMM for (short) multivariate time series from a sample with missing values

The problem in short:

I want to estimate (?) a lag-1 Markovian hidden process for offline multi-variate discrete-time time series with continuous distributions via smoothing, with no dimensionality reduction. The emission model is time invariant.

The data that I have:

I have a dataset of multivariate time series of $$N = 100$$ participants, each of $$8$$ variables (i.e. $$Y_{t_{8 \times 1}} = [y^1_{t}, y^2_{t}, ..., y^8_{t}]'$$) for $$T=50$$ time points, where I have around $$5\%$$ missing data. Missingness is most probably MCAR, but position and number of missig data differ from person to person.

The data is (rather highly) correlated ($$Cov(Y_t)$$ is absolutely non-diagonal) and the time series are autocorrelated but it drops significantly/sufficiently after lag-1.

What I need:

I have to find, for each participant, time series of $$X_{t_{8 \times 1}} = [x^1_{t}, x^2_{t}, ..., x^8_{t}]'$$ such that each $$x^i_{t}$$ is the lag-1 hidden Markov process of $$y^i_{t}$$, and it does not explain (or, is not contaminated by) any information of other $$y^j_{t}$$s.

I am skeptical that I can use HMM for this problem, at least since estimating the high number of parameters seems over-kill for this amount of data. So I assume Kalman filter would be a candidate here.

I have the following questions

1. For this problem, do you suggest Kalman filter or something else?

2. Does fitting a multivariate Kalman filter (or another model, based on your answer to the first question) prefered over multiple separate filters/mdels for each $$y^j_{t}$$? I think that it may help in estimating the distribution of $$\eta_t$$, $$\epsilon_t$$, and/or $$\alpha_t$$ (based on the notation in this answer.) However, as I said, I have to avoid the contamination mentioned above.

3. Can/should I fit them separately to each participant, or I shoul/have to fit to all the data of all people? (and how?)

4. Based on your answers so far, is this filtering/smoothing/model fitting feasible for the kind of data I have? I mean with this number of variables, data points, and missingness. I am in fact worried about overfitting or poor model fit.

5. What R package/function can I use for this problem?