0

我需要拟合如下状态空间模型:

y_t = Z\alpha_t + d + \epsilon_t; \quad \epsilon_t ~ N(0, diag(\sigma^2_\epsilon))

\alpha_{t+1} = T \alpha_t + \eta_t; \eta ~ N(0, I)

使用y_t一个d维向量(在我的情况下为d ~100)和 \alpha_t k维。

即,具有时不变矩阵和对角噪声观测协方差的“简单”线性状态空间模型。如果不是对角协方差,pykalman就可以完成这项工作。但是,据我了解,它不允许将观察协方差限制为对角线(您可以完全修复它或让它在没有约束的情况下进行估计)。

因此,我转向statsmodels并尝试自定义sm.tsa.statespace.MLEModel,但我不知道如何进行。我能找到的所有示例都导致要估计的标量参数很少,并且它们是通过该param_names方法显式声明的,但我不知道如何处理矩阵参数。

我正在做这样的事情:

class myStateSpace(MLEModel):
    """
    linear state space with diagonal observation covariance
    """

    def __init__(self, endog, k_states):

        # initialize state space model
        super().__init__(endog=endog, k_states=k_states, k_posdef=k_states,
        initialization='approximate_diffuse', loglikelihood_burn=2)

        self.ssm["state_intercept"] = np.zeros(self.k_states)
        self.ssm["state_cov"] = np.eye(self.k_states)
        self.ssm["selection"] = np.eye(self.k_states)
        # remain to be optimized:
        # design, obs_intercept, transition, and diag components of obs_cov 

        # tell the model that obs_cov will be diag
        self._obs_cov_idx = ("obs_cov", ) + np.diag_indices(self.k_endog)
        self._design_idx = ("design",) + np.indices((self.k_endog, self.k_states), sparse=True)
        self._transition_idx = ("transition",) + np.indices((self.k_states, self.k_states), sparse=True)
        self._obs_intercept_idx = ("obs_intercept",) + np.indices((self.k_endog,), sparse=True)


        # Which parameters need to be positive?
        # the last k_endog, namely the variances of obs_cov,
        # see param_names()
        self.num_params = self.k_states + self.k_states**2 + self.k_states*self.k_endog + self.k_endog 
        self.positive_parameters = np.array(range((self.num_params - self.k_endog), self.num_params))


    @property
    def param_names(self):
        design_names = ["design_%i%i" % (i, j) for i in range(self.K_endog) for j in range(self.k_states)]
        transition_names = ["transition_%i%i" % (i, j) for i in range(self.K_states) for j in range(self.k_states)]
        obs_intercept_names = ["state_intercept_%i" % i for i in range(self.k_endog)]
        obs_cov_names = ["obs_cov_%i" % i for i in range(self.K_endog)]

        return design_names + transition_names + obs_intercept_names + obs_cov_names 



    def transform_params(self, unconstrained):
            """
            constraint the variance parameters
            to be positive,
            """
            constrained = unconstrained.copy()
            constrained[self.positive_parameters] = constrained[self.positive_parameters]**2
            return constrained

    def untransform_params(self, constrained):
            """
            Need to unstransform all the parameters you transformed
            in the `transform_params` function
            """
            unconstrained = constrained.copy()
            unconstrained[self.positive_parameters] = unconstrained[self.positive_parameters]**0.5
            return unconstrained

    def update(self, params, **kwargs):
        params = super().update(params, **kwargs)
        design_n = (self.k_states * self.k_endog)
        transition_n = design_n + self.k_states**2
        obs_intercept_n = transition_n + self.K_endog


        self.ssm[self._design_idx] = params[: design_n].reshape((self.k_endog, self.k_states))
        self[self._transition_idx] = params[design_n : transition_n].reshape((self.k_states, self.k_states))
        self[self._obs_intercept_idx] = params[transition_n : obs_intercept_n]
        self[self._obs_cov_idx] = params[-self.k_endog:]

抛开我仍然怀念一个set_params方法,感觉这种做法真的很牵强:有没有更直接的方法来处理矩阵参数,而不必给矩阵的每个条目起个名字?

self.ssm["transition"]此外,使用,比如说,而不是self["transition"](我在例子中都看到过)有什么区别?

感谢您的任何建议!

4

0 回答 0