如何在statsmodels中定义含观测与未观测变量的观测矩阵？

我使用Python 3.9 + statsmodels 0.13包，通过卡尔曼滤波建模时间序列。状态空间转换矩阵如下：

观测矩阵依赖未观测变量v及观测变量p、n的历史值，如下：

我尝试自定义AR1类（继承自sm.tsa.statespace.MLEModel）实现模型，但无法添加观测方程中依赖psi/phi及过往观测值的首项，求可行解决思路。现有实现代码如下：

class AR1(sm.tsa.statespace.MLEModel):
    start_params = [-0.5, 0.5, 0.0, 0.05, 0.05]  # best guess at initial params
    param_names = ['-psi', '1-phi', 'r_t', 'e_t', 'w_t']
    
    def __init__(self, endog):
        # Initialize the state space model
        super(AR1, self).__init__(endog, k_states=2, k_posdef=1,
                                     initialization='stationary')

        # Setup the fixed components of the state space representation
        self['design'] = [[1., 1.],
                          [1., 0]]
        self['transition'] = [[1., 0],
                              [1., 0]]
        self['selection', 0, 0] = 1.

    # Describe how parameters enter the model
    def update(self, params, transformed=True, **kwargs):
        params = super(AR1, self).update(params, transformed, **kwargs)

        self['design', 0, 1] = params[0]        # param.0 is -psi
        self['design', 1, 0] = params[1]        # param.1 is (1-phi)
        self['selection', 0, 0] = params[2]     # param.2 is r_t  
        self['obs_cov', 0, 0] = params[3]       # param.3 is e_t
        self['obs_cov', 1, 1] = params[4]       # param.4 is w_t

    # Specify start parameters and parameter names
    @property
    def start_params(self):
        return self.start_params

# Create and fit the model
mod = AR1(DepVars)

# Display results
res = mod.fit()
print(res.summary())
res.plot_diagnostics(figsize=(13,7))

观测方程首项依赖过往观测值，属于动态观测方程（时变设计矩阵），statsmodels状态空间框架完全支持这类场景，下面提供两种实用方案：

方案1：扩展状态向量，纳入滞后观测值 #

1. 把滞后1期的p、n观测值加入状态向量，将原状态维度从2扩展到4（原v_t、x_t + 滞后p_{t-1}、n_{t-1}）
2. 调整转换矩阵：
   • 原v_t、x_t的转换规则保持不变
   • 滞后观测值的转换规则：当期滞后状态等于上一期的观测值，需要在update中用当期endog值动态更新这部分的转换逻辑
3. 调整设计矩阵：首项ψp_{t-1} + φn_{t-1}可通过设计矩阵对应滞后状态的系数（ψ、φ）直接实现

方案2：动态更新设计矩阵（更直接） #

1. 在__init__中将设计矩阵初始化为时变结构，比如提前创建一个形状为(T, k_endog, k_states)的空数组（T为时间序列长度）
2. 预处理得到滞后1期的观测值序列p_lag、n_lag，注意处理初始时间步的缺失值（可设为0或第一个观测值）
3. 在update方法中，遍历每个时间步，动态填充设计矩阵中依赖滞后观测值的部分：

   # 示例代码片段
   def __init__(self, endog):
       super().__init__(endog, k_states=2, k_posdef=1, initialization='stationary')
       # 初始化时变设计矩阵
       T = len(endog)
       self['design'] = np.zeros((T, 2, 2))
       # 设置固定项
       self['design', :, 0, 1] = 1.
       self['design', :, 1, 0] = 1.
       # 预计算滞后观测值
       self.p_lag = np.r_[endog[0,0], endog[:-1,0]]  # 用第一个值填充初始滞后项
       self.n_lag = np.r_[endog[0,1], endog[:-1,1]]
   
   def update(self, params, transformed=True, **kwargs):
       params = super().update(params, transformed, **kwargs)
       # 假设params新增psi、phi参数，索引为5、6
       self['design', :, 0, 0] = params[5] * self.p_lag + params[6] * self.n_lag
       # 其余参数更新逻辑保持不变
       self['design', 0, 1] = params[0]
       self['design', 1, 0] = params[1]
       self['selection', 0, 0] = params[2]
       self['obs_cov', 0, 0] = params[3]
       self['obs_cov', 1, 1] = params[4]
   

关键调整要点 #

• 新增psi、phi参数到start_params和param_names中，对应观测方程的系数
• 时变矩阵会增加少量计算量，但常规长度的时间序列完全可以承受
• 滞后观测值的初始值要合理设置，避免影响滤波初始步的稳定性

内容的提问来源于stack exchange，提问作者tom_servo