估计或指定状态空间模型中的参数¶
在本笔记本中，我们展示了如何在估计其他参数的同时，在statsmodels的状态空间模型中固定某些参数的特定值。
一般来说，状态空间模型允许用户：
通过最大似然估计所有参数
固定一些参数并估计其余部分
固定所有参数（以便不估计任何参数）
[1]:
%matplotlib inline
from importlib import reload
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from pandas_datareader.data import DataReader
为了说明，我们将使用服装的消费者价格指数，该指数具有时间变化的水平和强烈的季节性成分。
[2]:
endog = DataReader('CPIAPPNS', 'fred', start='1980').asfreq('MS')
endog.plot(figsize=(15, 3));
众所周知（例如，Harvey和Jaeger [1993]），在给定某些参数限制的情况下，HP滤波器的输出可以通过一个未观测分量模型生成。
未观测成分模型是：
\[\begin{split}\begin{aligned}
y_t & = \mu_t + \varepsilon_t & \varepsilon_t \sim N(0, \sigma_\varepsilon^2) \\
\mu_t &= \mu_{t-1} + \beta_{t-1} + \eta_t & \eta_t \sim N(0, \sigma_\eta^2) \\
\beta_t &= \beta_{t-1} + \zeta_t & \zeta_t \sim N(0, \sigma_\zeta^2) \\
\end{aligned}\end{split}\]
为了使趋势与HP滤波器的输出匹配，参数必须设置如下：
\[\begin{split}\begin{aligned}
\frac{\sigma_\varepsilon^2}{\sigma_\zeta^2} & = \lambda \\
\sigma_\eta^2 & = 0
\end{aligned}\end{split}\]
其中 \(\lambda\) 是相关HP滤波器的参数。对于我们在这里使用的月度数据，通常建议 \(\lambda = 129600\)。
[3]:
# Run the HP filter with lambda = 129600
hp_cycle, hp_trend = sm.tsa.filters.hpfilter(endog, lamb=129600)
# The unobserved components model above is the local linear trend, or "lltrend", specification
mod = sm.tsa.UnobservedComponents(endog, 'lltrend')
print(mod.param_names)
['sigma2.irregular', 'sigma2.level', 'sigma2.trend']
未观测成分模型（UCM）的参数写作：
\(\sigma_\varepsilon^2 = \text{sigma2.irregular}\)
\(\sigma_\eta^2 = \text{sigma2.level}\)
\(\sigma_\zeta^2 = \text{sigma2.trend}\)
为了满足上述限制，我们将设置 \((\sigma_\varepsilon^2, \sigma_\eta^2, \sigma_\zeta^2) = (1, 0, 1 / 129600)\)。
由于我们在这里固定了所有参数，我们根本不需要使用fit方法，因为该方法是用于执行最大似然估计的。相反，我们可以直接使用我们选择的参数运行卡尔曼滤波器和平滑器，使用smooth方法。
[4]:
res = mod.smooth([1., 0, 1. / 129600])
print(res.summary())
Unobserved Components Results
==============================================================================
Dep. Variable: CPIAPPNS No. Observations: 537
Model: local linear trend Log Likelihood -3005.996
Date: Wed, 16 Oct 2024 AIC 6017.992
Time: 18:27:47 BIC 6030.839
Sample: 01-01-1980 HQIC 6023.019
- 09-01-2024
Covariance Type: opg
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
sigma2.irregular 1.0000 0.009 115.625 0.000 0.983 1.017
sigma2.level 0 0.000 0 1.000 -0.000 0.000
sigma2.trend 7.716e-06 1.96e-07 39.281 0.000 7.33e-06 8.1e-06
===================================================================================
Ljung-Box (L1) (Q): 253.43 Jarque-Bera (JB): 1.65
Prob(Q): 0.00 Prob(JB): 0.44
Heteroskedasticity (H): 2.21 Skew: 0.04
Prob(H) (two-sided): 0.00 Kurtosis: 2.74
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
与HP滤波器的趋势估计相对应的估计值由水平的平滑估计给出（在上面的符号中为\(\mu_t\)）：
[5]:
ucm_trend = pd.Series(res.level.smoothed, index=endog.index)
可以看出，UCM对平滑水平的估计等于HP滤波器的输出：
[6]:
fig, ax = plt.subplots(figsize=(15, 3))
ax.plot(hp_trend, label='HP estimate')
ax.plot(ucm_trend, label='UCM estimate')
ax.legend();
添加季节性成分¶
然而，未观测成分模型比HP滤波器更加灵活。例如，上面显示的数据显然具有季节性，但季节性效应随时间变化（季节性在开始时比结束时弱得多）。未观测成分框架的优点之一是我们可以添加一个随机季节性成分。在这种情况下，我们将通过最大似然估计来估计季节性成分的方差，同时仍然包括上述对参数的限制，以使趋势对应于HP滤波器的概念。
添加随机季节性成分会增加一个新的参数，sigma2.seasonal。
[7]:
# Construct a local linear trend model with a stochastic seasonal component of period 1 year
mod = sm.tsa.UnobservedComponents(endog, 'lltrend', seasonal=12, stochastic_seasonal=True)
print(mod.param_names)
['sigma2.irregular', 'sigma2.level', 'sigma2.trend', 'sigma2.seasonal']
在这种情况下，我们将继续如上所述限制前三个参数，但我们希望通过最大似然估计来估计sigma2.seasonal的值。因此，我们将使用fit方法以及fix_params上下文管理器。
The fix_params 方法接受一个包含参数名称和关联值的字典。在生成的上下文中，这些参数将在所有情况下使用。在 fit 方法的情况下，只会估计未固定的参数。
[8]:
# Here we restrict the first three parameters to specific values
with mod.fix_params({'sigma2.irregular': 1, 'sigma2.level': 0, 'sigma2.trend': 1. / 129600}):
# Now we fit any remaining parameters, which in this case
# is just `sigma2.seasonal`
res_restricted = mod.fit()
RUNNING THE L-BFGS-B CODE
* * *
Machine precision = 2.220D-16
N = 1 M = 10
At X0 0 variables are exactly at the bounds
At iterate 0 f= 3.87461D+00 |proj g|= 2.27180D-01
At iterate 5 f= 3.25687D+00 |proj g|= 7.12208D-06
* * *
Tit = total number of iterations
Tnf = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip = number of BFGS updates skipped
Nact = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F = final function value
* * *
N Tit Tnf Tnint Skip Nact Projg F
1 5 12 1 0 0 7.122D-06 3.257D+00
F = 3.2568731346256001
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
This problem is unconstrained.
或者，我们可以简单地使用 fit_constrained 方法，该方法也接受一个约束字典：
[9]:
res_restricted = mod.fit_constrained({'sigma2.irregular': 1, 'sigma2.level': 0, 'sigma2.trend': 1. / 129600})
This problem is unconstrained.
RUNNING THE L-BFGS-B CODE
* * *
Machine precision = 2.220D-16
N = 1 M = 10
At X0 0 variables are exactly at the bounds
At iterate 0 f= 3.87461D+00 |proj g|= 2.27180D-01
At iterate 5 f= 3.25687D+00 |proj g|= 7.12208D-06
* * *
Tit = total number of iterations
Tnf = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip = number of BFGS updates skipped
Nact = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F = final function value
* * *
N Tit Tnf Tnint Skip Nact Projg F
1 5 12 1 0 0 7.122D-06 3.257D+00
F = 3.2568731346256001
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
摘要输出包括所有参数，但指出前三个参数是固定的（因此未被估计）。
[10]:
print(res_restricted.summary())
Unobserved Components Results
=====================================================================================
Dep. Variable: CPIAPPNS No. Observations: 537
Model: local linear trend Log Likelihood -1748.941
+ stochastic seasonal(12) AIC 3499.882
Date: Wed, 16 Oct 2024 BIC 3504.143
Time: 18:27:48 HQIC 3501.551
Sample: 01-01-1980
- 09-01-2024
Covariance Type: opg
============================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------------------
sigma2.irregular (fixed) 1.0000 nan nan nan nan nan
sigma2.level (fixed) 0 nan nan nan nan nan
sigma2.trend (fixed) 7.716e-06 nan nan nan nan nan
sigma2.seasonal 0.0924 0.007 12.672 0.000 0.078 0.107
===================================================================================
Ljung-Box (L1) (Q): 460.41 Jarque-Bera (JB): 38.25
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 2.39 Skew: 0.30
Prob(H) (two-sided): 0.00 Kurtosis: 4.18
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
作为比较，我们构建了无限制的最大似然估计（MLE）。在这种情况下，水平的估计将不再对应于HP滤波器的概念。
[11]:
res_unrestricted = mod.fit()
RUNNING THE L-BFGS-B CODE
* * *
Machine precision = 2.220D-16
N = 4 M = 10
At X0 0 variables are exactly at the bounds
At iterate 0 f= 3.63691D+00 |proj g|= 1.07263D-01
This problem is unconstrained.
At iterate 5 f= 1.99160D+00 |proj g|= 9.90158D-01
At iterate 10 f= 1.56578D+00 |proj g|= 5.36079D-01
At iterate 15 f= 1.50504D+00 |proj g|= 2.15906D-01
At iterate 20 f= 1.39950D+00 |proj g|= 2.61157D-01
At iterate 25 f= 1.38702D+00 |proj g|= 1.97478D-02
* * *
Tit = total number of iterations
Tnf = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip = number of BFGS updates skipped
Nact = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F = final function value
* * *
N Tit Tnf Tnint Skip Nact Projg F
4 29 55 1 0 0 4.866D-06 1.387D+00
F = 1.3868870900811454
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
最后，我们可以获取趋势和季节性成分的平滑估计值。
[12]:
# Construct the smoothed level estimates
unrestricted_trend = pd.Series(res_unrestricted.level.smoothed, index=endog.index)
restricted_trend = pd.Series(res_restricted.level.smoothed, index=endog.index)
# Construct the smoothed estimates of the seasonal pattern
unrestricted_seasonal = pd.Series(res_unrestricted.seasonal.smoothed, index=endog.index)
restricted_seasonal = pd.Series(res_restricted.seasonal.smoothed, index=endog.index)
比较估计的水平，可以清楚地看到，具有固定参数的季节性UCM仍然产生了一个与HP滤波器输出非常接近（尽管不再完全相同）的趋势。
同时，无参数限制模型（最大似然估计模型）的估计水平比这些要粗糙得多。
[13]:
fig, ax = plt.subplots(figsize=(15, 3))
ax.plot(unrestricted_trend, label='MLE, with seasonal')
ax.plot(restricted_trend, label='Fixed parameters, with seasonal')
ax.plot(hp_trend, label='HP filter, no seasonal')
ax.legend();
最后，具有参数限制的UCM仍然能够很好地捕捉到时间变化的季节性成分。
[14]:
fig, ax = plt.subplots(figsize=(15, 3))
ax.plot(unrestricted_seasonal, label='MLE')
ax.plot(restricted_seasonal, label='Fixed parameters')
ax.legend();
Last update:
Oct 16, 2024