I want to fit a linear Gaussian state space model that has three latent/unobservable factors and one observable factor. It is given as follows:
enter image description here

and the unobserved factors follow an autoregressive process:
enter image description here

here, G is some transition matrix of parameters, eta is assumed to have zero mean and constant covariance matrix. Dimension of delta x_t is n, b^o is an n-1 vector of parameters, f^o_t+1 is a scalar and is observed, B^u is the matrix of factor loadings which is known and specified, and f^u_t+1 is a vector in R^3 for the three unobservable factors.

I tried to use https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.dynamic_factor_mq.DynamicFactorMQ.html#statsmodels.tsa.statespace.dynamic_factor_mq.DynamicFactorMQ

but the issue is the observable factor, this appears to only work for when all f are unobserved. I want to fit this using the EM algorithm and am stuck. This should be standard and able to do with an existing package, maybe I am just misunderstanding.

So far I tried this from scratch but it doesn’t work.

from numpy.linalg import inv

def em_algorithm_state_space(delta_I, z_test, B_u, max_iter=100, tol=1e-3):
    """
    EM algorithm for fitting a linear Gaussian state space model.

    Args:
        delta_I: (n, k) array of observed data.
        z_test: (n, 1) array of the observed factor.
        B_u: (k, m) matrix of unobserved factor loadings.
        max_iter: Maximum number of EM iterations.
        tol: Convergence tolerance.

    Returns:
        G: Estimated transition matrix.
        Sigma_eta: Estimated covariance matrix of state innovations.
        b: Estimated intercept vector.
        Sigma_epsilon: Estimated covariance matrix of observation noise.
    """

    n, k = delta_I.shape
    m = B_u.shape[1]  # Number of unobserved factors

    # Initialize parameters
    G = np.zeros((m, m))
    Sigma_eta = np.eye(m)
    b = np.zeros((k, 1))
    Sigma_epsilon = np.eye(k)
    f_t = np.zeros((m, 1))  # Initial state vector

    for _ in range(max_iter):
        # E-step: Calculate conditional expectations of states
        f_t_cond = np.zeros((m, n))
        for t in range(1, n):
            # Prediction step (Kalman filter)
            f_t_pred = G @ f_t
            P_t_pred = G @ Sigma_eta @ G.T + Sigma_eta

            # Update step (Kalman filter)
            y_t = delta_I[t, :].reshape(-1, 1) - b - B_u @ f_t_pred - z_test[t, :].reshape(-1, 1)
            S_t = B_u @ P_t_pred @ B_u.T + Sigma_epsilon
            K_t = P_t_pred @ B_u.T @ inv(S_t)
            f_t = f_t_pred + K_t @ y_t
            P_t = (np.eye(m) - K_t @ B_u) @ P_t_pred

            f_t_cond[:, t] = f_t.reshape(-1)

        # M-step: Update parameters
        G = np.dot(f_t_cond[:, 1:] @ f_t_cond[:, :-1].T, inv(np.dot(f_t_cond[:, :-1] @ f_t_cond[:, :-1].T, G) + Sigma_eta))
        Sigma_eta = np.dot((f_t_cond[:, 1:] - G @ f_t_cond[:, :-1]) @ (f_t_cond[:, 1:] - G @ f_t_cond[:, :-1]).T, 1 / (n - 1))
        b = (delta_I - (B_u @ f_t_cond).T - z_test).mean(axis=0).reshape(-1, 1)
        residual = delta_I - b - B_u @ f_t_cond.T - np.tile(z_test.T, (1, k))
        Sigma_epsilon = np.dot(residual, residual.T) / n

    return G, Sigma_eta, b, Sigma_epsilon

# Example usage (replace with your actual data)
# delta_I = ...
# z_test = ...
# B_u = ...

# Fit the model
G_est, Sigma_eta_est, b_est, Sigma_epsilon_est = em_algorithm_state_space(delta_I_test, z_test, B_u)

print("Estimated G:", G_est)
print("Estimated Sigma_eta:", Sigma_eta_est)
print("Estimated b:", b_est)
print("Estimated Sigma_epsilon:", Sigma_epsilon_est)