State space models - concentrating the scale out of the likelihood function¶
[1]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
dta = sm.datasets.macrodata.load_pandas().data
dta.index = pd.date_range(start='1959Q1', end='2009Q4', freq='QE')
Introduction¶
(much of this is based on Harvey (1989); see especially section 3.4)
State space models can generically be written as follows (here we focus on time-invariant state space models, but similar results apply also to time-varying models):
Often, some or all of the values in the matrices \(Z, H, T, R, Q\) are unknown and must be estimated; in statsmodels, estimation is often done by finding the parameters that maximize the likelihood function. In particular, if we collect the parameters in a vector \(\psi\), then each of these matrices can be thought of as functions of those parameters, for example \(Z = Z(\psi)\), etc.
Usually, the likelihood function is maximized numerically, for example by applying quasi-Newton “hill-climbing” algorithms, and this becomes more and more difficult the more parameters there are. It turns out that in many cases we can reparameterize the model as \([\psi_*', \sigma_*^2]'\), where \(\sigma_*^2\) is the “scale” of the model (usually, it replaces one of the error variance terms) and it is possible to find the maximum likelihood estimate of \(\sigma_*^2\) analytically, by differentiating the likelihood function. This implies that numerical methods are only required to estimate the parameters \(\psi_*\), which has dimension one less than that of \(\psi\).
Example: local level model¶
(see, for example, section 4.2 of Harvey (1989))
As a specific example, consider the local level model, which can be written as:
In this model, \(Z, T,\) and \(R\) are all fixed to be equal to \(1\), and there are two unknown parameters, so that \(\psi = [\sigma_\varepsilon^2, \sigma_\eta^2]\).
Typical approach¶
First, we show how to define this model without concentrating out the scale, using statsmodels’ state space library:
[2]:
class LocalLevel(sm.tsa.statespace.MLEModel):
_start_params = [1., 1.]
_param_names = ['var.level', 'var.irregular']
def __init__(self, endog):
super(LocalLevel, self).__init__(endog, k_states=1, initialization='diffuse')
self['design', 0, 0] = 1
self['transition', 0, 0] = 1
self['selection', 0, 0] = 1
def transform_params(self, unconstrained):
return unconstrained**2
def untransform_params(self, unconstrained):
return unconstrained**0.5
def update(self, params, **kwargs):
params = super(LocalLevel, self).update(params, **kwargs)
self['state_cov', 0, 0] = params[0]
self['obs_cov', 0, 0] = params[1]
There are two parameters in this model that must be chosen: var.level \((\sigma_\eta^2)\) and var.irregular \((\sigma_\varepsilon^2)\). We can use the built-in fit method to choose them by numerically maximizing the likelihood function.
In our example, we are applying the local level model to consumer price index inflation.
[3]:
mod = LocalLevel(dta.infl)
res = mod.fit(disp=False)
print(res.summary())
Statespace Model Results
==============================================================================
Dep. Variable: infl No. Observations: 203
Model: LocalLevel Log Likelihood -457.632
Date: Thu, 14 Nov 2024 AIC 921.263
Time: 17:10:10 BIC 931.203
Sample: 03-31-1959 HQIC 925.285
- 09-30-2009
Covariance Type: opg
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
var.level 0.7447 0.156 4.766 0.000 0.438 1.051
var.irregular 3.3733 0.315 10.715 0.000 2.756 3.990
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 182.26
Prob(Q): 0.99 Prob(JB): 0.00
Heteroskedasticity (H): 1.75 Skew: -1.02
Prob(H) (two-sided): 0.02 Kurtosis: 7.18
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
We can look at the results from the numerical optimizer in the results attribute mle_retvals:
[4]:
print(res.mle_retvals)
{'fopt': np.float64(2.254343511382184), 'gopt': array([-7.10302928e-06, -9.72857350e-06]), 'fcalls': 27, 'warnflag': 0, 'converged': True, 'iterations': 7}
Concentrating out the scale¶
Now, there are two ways to reparameterize this model as above:
The first way is to set \(\sigma_*^2 \equiv \sigma_\varepsilon^2\) so that \(\psi_* = \psi / \sigma_\varepsilon^2 = [1, q_\eta]\) where \(q_\eta = \sigma_\eta^2 / \sigma_\varepsilon^2\).
The second way is to set \(\sigma_*^2 \equiv \sigma_\eta^2\) so that \(\psi_* = \psi / \sigma_\eta^2 = [h, 1]\) where \(h = \sigma_\varepsilon^2 / \sigma_\eta^2\).
In the first case, we only need to numerically maximize the likelihood with respect to \(q_\eta\), and in the second case we only need to numerically maximize the likelihood with respect to \(h\).
Either approach would work well in most cases, and in the example below we will use the second method.
To reformulate the model to take advantage of the concentrated likelihood function, we need to write the model in terms of the parameter vector \(\psi_* = [g, 1]\). Because this parameter vector defines \(\sigma_\eta^2 \equiv 1\), we now include a new line self['state_cov', 0, 0] = 1 and the only unknown parameter is \(h\). Because our parameter \(h\) is no longer a variance, we renamed it here to be ratio.irregular.
The key piece that is required to formulate the model so that the scale can be computed from the Kalman filter recursions (rather than selected numerically) is setting the flag self.ssm.filter_concentrated = True.
[5]:
class LocalLevelConcentrated(sm.tsa.statespace.MLEModel):
_start_params = [1.]
_param_names = ['ratio.irregular']
def __init__(self, endog):
super(LocalLevelConcentrated, self).__init__(endog, k_states=1, initialization='diffuse')
self['design', 0, 0] = 1
self['transition', 0, 0] = 1
self['selection', 0, 0] = 1
self['state_cov', 0, 0] = 1
self.ssm.filter_concentrated = True
def transform_params(self, unconstrained):
return unconstrained**2
def untransform_params(self, unconstrained):
return unconstrained**0.5
def update(self, params, **kwargs):
params = super(LocalLevelConcentrated, self).update(params, **kwargs)
self['obs_cov', 0, 0] = params[0]
Again, we can use the built-in fit method to find the maximum likelihood estimate of \(h\).
[6]:
mod_conc = LocalLevelConcentrated(dta.infl)
res_conc = mod_conc.fit(disp=False)
print(res_conc.summary())
Statespace Model Results
==================================================================================
Dep. Variable: infl No. Observations: 203
Model: LocalLevelConcentrated Log Likelihood -457.632
Date: Thu, 14 Nov 2024 AIC 921.263
Time: 17:10:10 BIC 931.203
Sample: 03-31-1959 HQIC 925.285
- 09-30-2009 Scale 0.745
Covariance Type: opg
===================================================================================
coef std err z P>|z| [0.025 0.975]
-----------------------------------------------------------------------------------
ratio.irregular 4.5297 1.226 3.694 0.000 2.126 6.933
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 182.26
Prob(Q): 0.99 Prob(JB): 0.00
Heteroskedasticity (H): 1.75 Skew: -1.02
Prob(H) (two-sided): 0.02 Kurtosis: 7.18
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
The estimate of \(h\) is provided in the middle table of parameters (ratio.irregular), while the estimate of the scale is provided in the upper table. Below, we will show that these estimates are consistent with those from the previous approach.
And we can again look at the results from the numerical optimizer in the results attribute mle_retvals. It turns out that two fewer iterations were required in this case, since there was one fewer parameter to select. Moreover, since the numerical maximization problem was easier, the optimizer was able to find a value that made the gradient for this parameter slightly closer to zero than it was above.
[7]:
print(res_conc.mle_retvals)
{'fopt': np.float64(2.2543435111703576), 'gopt': array([-6.71906974e-08]), 'fcalls': 12, 'warnflag': 0, 'converged': True, 'iterations': 5}
Comparing estimates¶
Recall that \(h = \sigma_\varepsilon^2 / \sigma_\eta^2\) and the scale is \(\sigma_*^2 = \sigma_\eta^2\). Using these definitions, we can see that both models produce nearly identical results:
[8]:
print('Original model')
print(f'var.level = {res.params.iloc[0]:.5f}')
print(f'var.irregular = {res.params.iloc[1]:.5f}')
print('\nConcentrated model')
print(f'scale = {res_conc.scale:.5f}' )
print(f'h * scale = {res_conc.params.iloc[0] * res_conc.scale:.5f}')
Original model
var.level = 0.74469
var.irregular = 3.37330
Concentrated model
scale = 0.74472
h * scale = 3.37338
Example: SARIMAX¶
By default in SARIMAX models, the variance term is chosen by numerically maximizing the likelihood function, but an option has been added to allow concentrating the scale out.
[9]:
# Typical approach
mod_ar = sm.tsa.SARIMAX(dta.cpi, order=(1, 0, 0), trend='ct')
res_ar = mod_ar.fit(disp=False)
# Estimating the model with the scale concentrated out
mod_ar_conc = sm.tsa.SARIMAX(dta.cpi, order=(1, 0, 0), trend='ct', concentrate_scale=True)
res_ar_conc = mod_ar_conc.fit(disp=False)
These two approaches produce about the same loglikelihood and parameters, although the model with the concentrated scale was able to improve the fit very slightly:
[10]:
print('Loglikelihood')
print(f'- Original model: {res_ar.llf:.4f}' )
print(f'- Concentrated model: {res_ar_conc.llf:.4f}')
print('\nParameters')
print(f'- Original model:')
for idx in res_ar.params.index:
print(f'{idx: <12} {res_ar.params[idx]:.4f}')
print('- Concentrated model:')
for idx in res_ar_conc.params.index:
print(f'{idx: <12} {res_ar_conc.params[idx]:.4f}')
print(f'{"sigma2": <12} {res_ar_conc.scale:.4f}')
Loglikelihood
- Original model: -245.8275
- Concentrated model: -245.8264
Parameters
- Original model:
intercept 0.4921
drift 0.0243
ar.L1 0.9808
sigma2 0.6490
- Concentrated model:
intercept 0.4864
drift 0.0242
ar.L1 0.9809
sigma2 0.6492
This time, about 1/3 fewer iterations of the optimizer are required under the concentrated approach:
[11]:
print('Optimizer iterations')
print(f'- Original model: {res_ar.mle_retvals["iterations"]}')
print(f'- Concentrated model: {res_ar_conc.mle_retvals["iterations"]}')
Optimizer iterations
- Original model: 36
- Concentrated model: 22