Подбор линейной модели пространства состояний с наблюдаемыми и ненаблюдаемыми факторами в Python ⇐ Python

[Anonymous]

Я хочу подобрать линейную модель гауссовского пространства состояний, которая имеет три скрытых/ненаблюдаемых фактора и один наблюдаемый фактор. Оно задается следующим образом:
как здесь
и ненаблюдаемые факторы следуют авторегрессионному процессу:
как здесь
здесь G — некоторая матрица перехода параметров, предполагается, что эта имеет нулевое среднее и постоянную ковариационную матрицу. Размерность дельты x_t равна n, b^o — вектор параметров n-1, f^o_t+1 — скаляр и наблюдается, B^u — матрица (n-1)x3 факторных нагрузок, которая известна и указано, а f^u_t+1 — это вектор в R^3 для трех ненаблюдаемых факторов.
Я пытался использовать https://www.statsmodels.org/stable/generated/ statsmodels.tsa.statespace.dynamic_factor_mq.DynamicFactorMQ.html#statsmodels.tsa.statespace.dynamic_factor_mq.DynamicFactorMQ
но проблема в наблюдаемом факторе, похоже, это работает только тогда, когда все f незаметно. Я хочу подогнать это под алгоритм EM и застрял. Это должно быть стандартно и применимо к существующему пакету, возможно, я просто неправильно понимаю.
Я пробовал это с нуля, но это не работает.
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)


Подробнее здесь: https://stackoverflow.com/questions/787 ... ctors-in-p