问题描述
我需要拟合如下状态空间模型:
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,sparse=True)
self._obs_intercept_idx = ("obs_intercept",),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"]
(我在示例中都见过)有什么区别?
感谢您的任何建议!
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)