我需要拟合如下状态空间模型:
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"]
(我在例子中都看到过)有什么区别?
感谢您的任何建议!