statsmodels 状态 space 模型中的平滑状态干扰
smoothed state disturbances in statsmodels' state space model
我正在构建自定义状态-space 模型，其中包含如下统计模型
from statsmodels.tsa.statespace import mlemodel
mod = mlemodel.MLEModel(YY, k_states=n_st, k_posdef=n_sh)
mod['design'] = Z
mod['transition'] = T
mod['selection'] = R
mod['state_cov'] = np.eye(n_sh)
mod['obs_intercept'] = d
mod.initialize_stationary()
我对平滑状态和平滑状态扰动很感兴趣，我得到了
results = mod.smooth([])
results.smoothed_state 中的平滑状态是正确的（我有我正在比较的真实值），但 results.smoothed_state_disturbance 中的平滑状态扰动向前移动了一个周期 - 第一列包含第二个周期的（正确的）平滑扰动等，而最后一列包含零，这是样本结束后一个周期的正确平滑扰动。
我的理解是，这与状态方程的时间有关，根据statsmodels docs here是
alpha(t+1) = T alpha(t) + R eta(t) (1)
因此意味着第一个观测值 y_{1} 与状态 alpha_{1} 有关，而状态 alpha_{1} 又取决于扰动 eta_{0}，并且该扰动的平滑值不会由顺畅。另一方面，在this statmodels docs中，状态方程的时序为
alpha(t) = T alpha(t-1) + R eta(t) (2)
并暗示状态 alpha_{1} 取决于 eta_{1}，而不是 eta_{0}。由于（未来 (1) 和同期 (2)）时间约定都出现在 statsmodels 文档中，我认为可以选择使用哪一个。不幸的是，一直无法找出方法。我尝试使用 results = mod.smooth([], filter_timing=1) 更改更平滑的时间，根据文档使用 Kim 和 Nelson (1999)（同期）时间，而不是默认的 Durbin 和 Koopman (2012)（未来）时间。但是后来我得到了完全不同的结果（而且是错误的，因为我知道真正的价值是什么）不仅来自更平滑的结果，还来自对数似然值的结果。我还在 the unit tests for smoothing 中查找示例，但只有针对 MATLAB 和 R 库的测试也使用未来时序，并且没有针对 STATA 的测试（用于干扰平滑），后者使用替代同期时序。
我的问题是，有没有办法用同期时间（上面的（2））写出状态方程，或者恢复与第一周期观测数据相关的平滑状态扰动。
这是以下带有测量误差的 AR(1) 模型的一些代码，对状态方程使用 同期 时序，使用平稳分布进行初始化。
alpha(0) ~ N(0, 1/(1-.5**2))
alpha(t) = .5 alpha(t-1) + eta(t), eta(t) ~ N(0, 1)
y(t) = alpha(t) + e(t), e(t) ~ N(0, 1)
from statsmodels.tsa.statespace import mlemodel
import numpy as np
import sys
from scipy.stats import multivariate_normal
from numpy.random import default_rng
gen = default_rng(42)
T = np.array([.5])
Z = np.array([1.])
Q = np.array([1.])
H = np.array([1.])
R = np.array([1.])
P0 = 1/(1-T**2)
# Simulate data for 2 periods
alpha0 = gen.normal(0, np.sqrt(P0))
eta1 = gen.normal(0, 1)
e1 = gen.normal(0, 1)
eta2 = gen.normal(0, 1)
e2 = gen.normal(0, 1)
alpha1 = .5*alpha0 + eta1
y1 = alpha1 + e1
alpha2 = .5*alpha1 + eta2
y2 = alpha2 + e2
首先，使用 statsmodels.statespace 计算平滑状态、平滑状态扰动和仅给定第一个数据点的对数似然
mod1 = mlemodel.MLEModel(y1, k_states=1, k_posdef=1)
mod1['design'] = Z
mod1['transition'] = T
mod1['selection'] = R
mod1['state_cov'] = Q
mod1['obs_cov'] = H
mod1.initialize_stationary()
results1 = mod1.smooth([])
results1.smoothed_state, results1.smoothed_state_disturbance, results1.llf
给予
(array([[-0.06491681]]), array([[0.]]), -1.3453530272821392)
注意，观察 y(1) 我们可以计算 eta(1) 的条件期望，然而，这里返回的是什么只是 eta(2) 的条件期望。由于模型是平稳的和高斯的，alpha(1) 和 eta(1) 的条件期望给定 y(1 )可以从它们的联合分布计算出来（相关公式见here），如下代码
# Define a matrix L1 which maps [alpha(0), eta(1), e(1)] into [alpha0, eta1, e1, alpha1, y1]
L1 = np.vstack((np.eye(3), # alpha(0), eta(1), e(1)
np.r_[T, 1, 0], # alpha(1)
np.r_[T, 1, 1], # y(1)
))
# check
np.testing.assert_array_equal(np.r_[alpha0, eta1, e1, alpha1, y1],
L1 @ np.r_[alpha0, eta1, e1])
# Compute Sigma1 as the covariance matrix of [alpha0, eta1, e1, alpha1, y1]
D1 = np.eye(3)
D1[0, 0] = P0
Sigma1 = L1 @ D1 @ L1.T
# [alpha0, eta1, e1, alpha1, y1] has a multivariate Normal distribution, and we can apply well-known formulae to compute conditional expectations and the log-likelihood
ind_e1 = 1
ind_eta1 = 2
ind_alpha1 = 3
ind_y1 = 4
smooth_eta1 = (Sigma1[ind_eta1, ind_y1]/Sigma1[ind_y1, ind_y1])*y1
smooth_alpha1 = (Sigma1[ind_alpha1, ind_y1]/Sigma1[ind_y1, ind_y1])*y1
loglik1 = multivariate_normal.logpdf(y1, cov=Sigma1[ind_y1, ind_y1])
smooth_alpha1, smooth_eta1, loglik1
这给出了
(array([-0.06491681]), array([-0.04868761]), -1.3453530272821392)
扩展到前 2 个周期，statsmodels
y = np.array([y1, y2])
mod2 = mlemodel.MLEModel(y, k_states=1, k_posdef=1)
mod2.ssm.timing_init_filtered = True
mod2['design'] = Z
mod2['transition'] = T
mod2['selection'] = R
mod2['state_cov'] = Q
mod2['obs_cov'] = H
mod2.initialize_stationary()
results2 = mod2.smooth([])
results2.smoothed_state, results2.smoothed_state_disturbance, results2.llf
给予
(array([[-0.25292213, -0.78447967]]),
array([[-0.65801861, 0. ]]),
-3.1092778246103645)
并计算联合分布的条件期望
# L2 maps [alpha(0), eta(1), e(1), eta(2), e(2)] into [alpha0, eta1, e1, eta2, e2, alpha1, alpha2, y1, y2]
L2 = np.vstack((np.eye(5), # alpha(0), eta(1), e(1), eta(2), e(2)
np.r_[T, 1, 0, 0, 0], # alpha(1)
np.r_[T**2, T, 0, 1, 0], # alpha(2)
np.r_[T, 1, 1, 0, 0], # y(1)
np.r_[T**2, T, 0, 1, 1], # y(2)
))
np.testing.assert_array_equal(np.r_[alpha0, eta1, e1, eta2, e2, alpha1, alpha2, y1, y2],
L2 @ np.r_[alpha0, eta1, e1, eta2, e2,])
# Sigma2 is the covariance of [alpha0, eta1, e1, eta2, e2, alpha1, alpha2, y1, y2]
D2 = np.eye(5)
D2[0, 0] = P0
Sigma2 = L2 @ D2 @ L2.T
ind_e = [2, 4]
ind_eta = [1, 3]
ind_alpha = [5, 6]
ind_y = [7, 8]
# compute smoothed disturbances and states, and loglikelihood
smooth_eta = Sigma2[ind_eta, :][:, ind_y] @ np.linalg.solve(Sigma2[ind_y, :][:, ind_y], y)
smooth_alpha = Sigma2[ind_alpha, :][:, ind_y] @ np.linalg.solve(Sigma2[ind_y, :][:, ind_y], y)
loglik2 = multivariate_normal.logpdf(y.flatten(), cov=Sigma2[ind_y, :][:, ind_y])
smooth_alpha.flatten(), smooth_eta.flatten(), loglik2
给予
(array([-0.25292213, -0.78447967]),
array([-0.1896916 , -0.65801861]),
-3.1092778246103636)
平滑状态 alpha(t) 和对数似然值相同。 statsmodels.statespace.mlemodel 返回的平滑干扰用于 eta(2) 和 eta(3).
最短答案
对您的问题的简短回答是，您可以执行以下操作来恢复您正在寻找的平滑估计：
r_0 = results1.smoother_results.scaled_smoothed_estimator_presample
R = mod1['selection']
Q = mod1['state_cov']
eta_hat_0 = Q @ R.T @ r_0
print(eta_hat_0)
这会如您所愿地给出 [-0.04868761]。这来自周期 0 的通常扰动平滑方程（例如 Durbin 和 Koopman 2012 年的书，方程 4.69）。在 Statsmodels 中，我们存储样本周期 1 - nobs 的平滑扰动估计，但我们确实存储 r_0，如上所示，所以你可以直接计算出你想要的值。
替代方法，也许更简单
这是获得此输出的另一种方法，我们可以通过以不同的方式编写您的问题来得出该输出。正如您所指出的，在 Statsmodels 的时间假设下，第一状态 alpha(1) ~ N(a1, P1) 的分布被指定为先验。我们可以将您想要的模型视为指定第 0 个状态 alpha(0) ~ N(a0, P0) 作为先验，然后将第一个状态视为由转移方程产生：alpha(1) = T alpha(0 ) + eta(0).
这是现在使用 Statsmodels 的错误约定编写的，我们可以使用 Statsmodels 来计算平滑的结果。唯一的技巧是没有与第一个状态 alpha(0) 关联的观察 y(0)，因为我们只包括 alpha(0) -> alpha(1) 转换步骤，以便我们可以像您一样指定先验通缉。但这没问题，我们可以简单地包含一个缺失值。
因此，如果我们只是通过将 nan 预采样值放入输入数据来修改您的原始模型，我们将得到我们想要的结果：
mod1 = mlemodel.MLEModel(np.r_[np.nan, y1], k_states=1, k_posdef=1)
mod1['design'] = Z
mod1['transition'] = T
mod1['selection'] = R
mod1['state_cov'] = Q
mod1['obs_cov'] = H
mod1.initialize_stationary()
results1 = mod1.smooth([])
print(results1.smoothed_state, results1.smoothed_state_disturbance, results1.llf)
产量：
[[-0.03245841 -0.06491681]] [[-0.04868761 0. ]] -1.3453530272821392
关于 filter_timing 标志的注意事项:
filter_timing 标志还允许您更改先验的计时约定，但不能更改状态。如果在备选时间（filter_timing=1）中设置先验 alpha(0|0) ~ N(a_00, P_00)，则设置 alpha(1) ~ N(a_1, P_1) with a_1 = T a_00 and P_1 = T P_00 T' + R Q R' in the default timing (filter_timing=0 ) 会给你完全相同的结果。
我正在构建自定义状态-space 模型，其中包含如下统计模型
from statsmodels.tsa.statespace import mlemodel
mod = mlemodel.MLEModel(YY, k_states=n_st, k_posdef=n_sh)
mod['design'] = Z
mod['transition'] = T
mod['selection'] = R
mod['state_cov'] = np.eye(n_sh)
mod['obs_intercept'] = d
mod.initialize_stationary()
我对平滑状态和平滑状态扰动很感兴趣，我得到了
results = mod.smooth([])
results.smoothed_state 中的平滑状态是正确的（我有我正在比较的真实值），但 results.smoothed_state_disturbance 中的平滑状态扰动向前移动了一个周期 - 第一列包含第二个周期的（正确的）平滑扰动等，而最后一列包含零，这是样本结束后一个周期的正确平滑扰动。
我的理解是，这与状态方程的时间有关，根据statsmodels docs here是
alpha(t+1) = T alpha(t) + R eta(t) (1)
因此意味着第一个观测值 y_{1} 与状态 alpha_{1} 有关，而状态 alpha_{1} 又取决于扰动 eta_{0}，并且该扰动的平滑值不会由顺畅。另一方面，在this statmodels docs中，状态方程的时序为
alpha(t) = T alpha(t-1) + R eta(t) (2)
并暗示状态 alpha_{1} 取决于 eta_{1}，而不是 eta_{0}。由于（未来 (1) 和同期 (2)）时间约定都出现在 statsmodels 文档中，我认为可以选择使用哪一个。不幸的是，一直无法找出方法。我尝试使用 results = mod.smooth([], filter_timing=1) 更改更平滑的时间，根据文档使用 Kim 和 Nelson (1999)（同期）时间，而不是默认的 Durbin 和 Koopman (2012)（未来）时间。但是后来我得到了完全不同的结果（而且是错误的，因为我知道真正的价值是什么）不仅来自更平滑的结果，还来自对数似然值的结果。我还在 the unit tests for smoothing 中查找示例，但只有针对 MATLAB 和 R 库的测试也使用未来时序，并且没有针对 STATA 的测试（用于干扰平滑），后者使用替代同期时序。
我的问题是，有没有办法用同期时间（上面的（2））写出状态方程，或者恢复与第一周期观测数据相关的平滑状态扰动。
这是以下带有测量误差的 AR(1) 模型的一些代码，对状态方程使用 同期 时序，使用平稳分布进行初始化。
alpha(0) ~ N(0, 1/(1-.5**2))
alpha(t) = .5 alpha(t-1) + eta(t), eta(t) ~ N(0, 1)
y(t) = alpha(t) + e(t), e(t) ~ N(0, 1)
from statsmodels.tsa.statespace import mlemodel
import numpy as np
import sys
from scipy.stats import multivariate_normal
from numpy.random import default_rng
gen = default_rng(42)
T = np.array([.5])
Z = np.array([1.])
Q = np.array([1.])
H = np.array([1.])
R = np.array([1.])
P0 = 1/(1-T**2)
# Simulate data for 2 periods
alpha0 = gen.normal(0, np.sqrt(P0))
eta1 = gen.normal(0, 1)
e1 = gen.normal(0, 1)
eta2 = gen.normal(0, 1)
e2 = gen.normal(0, 1)
alpha1 = .5*alpha0 + eta1
y1 = alpha1 + e1
alpha2 = .5*alpha1 + eta2
y2 = alpha2 + e2
首先，使用 statsmodels.statespace 计算平滑状态、平滑状态扰动和仅给定第一个数据点的对数似然
mod1 = mlemodel.MLEModel(y1, k_states=1, k_posdef=1)
mod1['design'] = Z
mod1['transition'] = T
mod1['selection'] = R
mod1['state_cov'] = Q
mod1['obs_cov'] = H
mod1.initialize_stationary()
results1 = mod1.smooth([])
results1.smoothed_state, results1.smoothed_state_disturbance, results1.llf
给予
(array([[-0.06491681]]), array([[0.]]), -1.3453530272821392)
注意，观察 y(1) 我们可以计算 eta(1) 的条件期望，然而，这里返回的是什么只是 eta(2) 的条件期望。由于模型是平稳的和高斯的，alpha(1) 和 eta(1) 的条件期望给定 y(1 )可以从它们的联合分布计算出来（相关公式见here），如下代码
# Define a matrix L1 which maps [alpha(0), eta(1), e(1)] into [alpha0, eta1, e1, alpha1, y1]
L1 = np.vstack((np.eye(3), # alpha(0), eta(1), e(1)
np.r_[T, 1, 0], # alpha(1)
np.r_[T, 1, 1], # y(1)
))
# check
np.testing.assert_array_equal(np.r_[alpha0, eta1, e1, alpha1, y1],
L1 @ np.r_[alpha0, eta1, e1])
# Compute Sigma1 as the covariance matrix of [alpha0, eta1, e1, alpha1, y1]
D1 = np.eye(3)
D1[0, 0] = P0
Sigma1 = L1 @ D1 @ L1.T
# [alpha0, eta1, e1, alpha1, y1] has a multivariate Normal distribution, and we can apply well-known formulae to compute conditional expectations and the log-likelihood
ind_e1 = 1
ind_eta1 = 2
ind_alpha1 = 3
ind_y1 = 4
smooth_eta1 = (Sigma1[ind_eta1, ind_y1]/Sigma1[ind_y1, ind_y1])*y1
smooth_alpha1 = (Sigma1[ind_alpha1, ind_y1]/Sigma1[ind_y1, ind_y1])*y1
loglik1 = multivariate_normal.logpdf(y1, cov=Sigma1[ind_y1, ind_y1])
smooth_alpha1, smooth_eta1, loglik1
这给出了
(array([-0.06491681]), array([-0.04868761]), -1.3453530272821392)
扩展到前 2 个周期，statsmodels
y = np.array([y1, y2])
mod2 = mlemodel.MLEModel(y, k_states=1, k_posdef=1)
mod2.ssm.timing_init_filtered = True
mod2['design'] = Z
mod2['transition'] = T
mod2['selection'] = R
mod2['state_cov'] = Q
mod2['obs_cov'] = H
mod2.initialize_stationary()
results2 = mod2.smooth([])
results2.smoothed_state, results2.smoothed_state_disturbance, results2.llf
给予
(array([[-0.25292213, -0.78447967]]),
array([[-0.65801861, 0. ]]),
-3.1092778246103645)
并计算联合分布的条件期望
# L2 maps [alpha(0), eta(1), e(1), eta(2), e(2)] into [alpha0, eta1, e1, eta2, e2, alpha1, alpha2, y1, y2]
L2 = np.vstack((np.eye(5), # alpha(0), eta(1), e(1), eta(2), e(2)
np.r_[T, 1, 0, 0, 0], # alpha(1)
np.r_[T**2, T, 0, 1, 0], # alpha(2)
np.r_[T, 1, 1, 0, 0], # y(1)
np.r_[T**2, T, 0, 1, 1], # y(2)
))
np.testing.assert_array_equal(np.r_[alpha0, eta1, e1, eta2, e2, alpha1, alpha2, y1, y2],
L2 @ np.r_[alpha0, eta1, e1, eta2, e2,])
# Sigma2 is the covariance of [alpha0, eta1, e1, eta2, e2, alpha1, alpha2, y1, y2]
D2 = np.eye(5)
D2[0, 0] = P0
Sigma2 = L2 @ D2 @ L2.T
ind_e = [2, 4]
ind_eta = [1, 3]
ind_alpha = [5, 6]
ind_y = [7, 8]
# compute smoothed disturbances and states, and loglikelihood
smooth_eta = Sigma2[ind_eta, :][:, ind_y] @ np.linalg.solve(Sigma2[ind_y, :][:, ind_y], y)
smooth_alpha = Sigma2[ind_alpha, :][:, ind_y] @ np.linalg.solve(Sigma2[ind_y, :][:, ind_y], y)
loglik2 = multivariate_normal.logpdf(y.flatten(), cov=Sigma2[ind_y, :][:, ind_y])
smooth_alpha.flatten(), smooth_eta.flatten(), loglik2
给予
(array([-0.25292213, -0.78447967]),
array([-0.1896916 , -0.65801861]),
-3.1092778246103636)
平滑状态 alpha(t) 和对数似然值相同。 statsmodels.statespace.mlemodel 返回的平滑干扰用于 eta(2) 和 eta(3).
最短答案
对您的问题的简短回答是，您可以执行以下操作来恢复您正在寻找的平滑估计：
r_0 = results1.smoother_results.scaled_smoothed_estimator_presample
R = mod1['selection']
Q = mod1['state_cov']
eta_hat_0 = Q @ R.T @ r_0
print(eta_hat_0)
这会如您所愿地给出 [-0.04868761]。这来自周期 0 的通常扰动平滑方程（例如 Durbin 和 Koopman 2012 年的书，方程 4.69）。在 Statsmodels 中，我们存储样本周期 1 - nobs 的平滑扰动估计，但我们确实存储 r_0，如上所示，所以你可以直接计算出你想要的值。
替代方法，也许更简单
这是获得此输出的另一种方法，我们可以通过以不同的方式编写您的问题来得出该输出。正如您所指出的，在 Statsmodels 的时间假设下，第一状态 alpha(1) ~ N(a1, P1) 的分布被指定为先验。我们可以将您想要的模型视为指定第 0 个状态 alpha(0) ~ N(a0, P0) 作为先验，然后将第一个状态视为由转移方程产生：alpha(1) = T alpha(0 ) + eta(0).
这是现在使用 Statsmodels 的错误约定编写的，我们可以使用 Statsmodels 来计算平滑的结果。唯一的技巧是没有与第一个状态 alpha(0) 关联的观察 y(0)，因为我们只包括 alpha(0) -> alpha(1) 转换步骤，以便我们可以像您一样指定先验通缉。但这没问题，我们可以简单地包含一个缺失值。
因此，如果我们只是通过将 nan 预采样值放入输入数据来修改您的原始模型，我们将得到我们想要的结果：
mod1 = mlemodel.MLEModel(np.r_[np.nan, y1], k_states=1, k_posdef=1)
mod1['design'] = Z
mod1['transition'] = T
mod1['selection'] = R
mod1['state_cov'] = Q
mod1['obs_cov'] = H
mod1.initialize_stationary()
results1 = mod1.smooth([])
print(results1.smoothed_state, results1.smoothed_state_disturbance, results1.llf)
产量：
[[-0.03245841 -0.06491681]] [[-0.04868761 0. ]] -1.3453530272821392
关于 filter_timing 标志的注意事项:
filter_timing 标志还允许您更改先验的计时约定，但不能更改状态。如果在备选时间（filter_timing=1）中设置先验 alpha(0|0) ~ N(a_00, P_00)，则设置 alpha(1) ~ N(a_1, P_1) with a_1 = T a_00 and P_1 = T P_00 T' + R Q R' in the default timing (filter_timing=0 ) 会给你完全相同的结果。