# initial conditions for statespace mlemodel in statsmodel

I am a bit puzzled by very sensitive dependence on the initial conditions in the statespace mlemodels in statsmodel. Let me take a concrete example here. I am trying to fit this Dynamic Linear model with 5 parameters $$(\phi_1, \phi_2, \sigma_{\epsilon}, \sigma_z, \sigma_u)$$. Here are the design and transition equations.

\begin{align*} Y_t &= \begin{bmatrix}1 & 0 & 1 & 0\end{bmatrix}\begin{bmatrix}Y_t^{(p)} \\ \delta_t \\ Y_t^{(g)} \\ \theta_{t,2}^{(g)} \end{bmatrix} + v_t, & v \sim \mathcal N (0, 0) \\ \begin{bmatrix}Y_t^{(p)} \\ \delta_t \\ Y_t^{(g)} \\ \theta_{t,2}^{(g)} \end{bmatrix} &= \begin{bmatrix}1 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \phi_1 & 1 \\ 0 & 0& \phi_2 & 0\end{bmatrix}\begin{bmatrix}Y_{t-1}^{(p)} \\ \delta_{t-1} \\ Y_{t-1}^{(g)} \\ \theta_{t-1,2}^{(g)} \end{bmatrix}+w_t, & w_t \sim \mathcal N\left(0, \begin{bmatrix}\sigma^2_{\epsilon}&0&0&0\\0&\sigma^2_z&0&0\\0&0&\sigma^2_u&0\\0&0&0&0\end{bmatrix}\right) \end{align*}

I am providing code to get the data which is US GDP. I try various 'start_params' values and different start.loglikelihood_burn values and get different results.

import pandas_datareader as pdr
dff = pdr.get_data_fred('GDPC1', '19500101', '20041001').apply(np.log)

class GDPGAP(sm.tsa.statespace.MLEModel):
start_params = [.1, .1, 100, 100, 100]
def __init__(self, endog):
super(GDPGAP, self).__init__(endog, k_states=4)
self['design'] = np.array([[1,0,1,0]])  # F
self['transition'] = np.array([[1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,1,0]])  # G
self['selection'] = np.eye(4)  # R=1
self['obs_cov', 0, 0] = np.array([[1e-8]])  # V
self['state_cov'] = np.diag([1, 1, 1, 0]) # W
self.initialize_approximate_diffuse()
self.loglikelihood_burn = 1

def update(self, params, **kwargs):
self['transition', 2,2] = params[0]
self['transition', 3,2] = params[1]
self['state_cov'] = np.diag([params[2]**2, params[3]**2, params[4]**2, 0])  # W

ggmodel = GDPGAP(dff)
res = ggmodel.fit(maxiter=1000)
res.summary()



I show various results of the parameters for different values tried

start_params, loglikelihood_burn
[.1, .1, 100, 100, 100],1 : [ 1.5525, -0.5525,  0.0059, -0.    ,  0.0061]
[.1, -.1, 100, 100, 100],1: [  1.9859,  -0.9858,  99.8563,  99.907 , 100.009 ]
[.1, .1, 100, 100, 100],2: [  1.9983,  -0.9982, 100.1304, 100.0526,  99.8045]


Why is it so? Is there a way to understand the stability of these and advice on how to set it up. The reference I am using shows that the right answer should be $$[1.481, -0.547, 0.00578, 0.00008, 0.00615]$$.

EDIT: following up on constraining the off diagonal terms. I tried the following and it does not seem to converge.

from statsmodels.tsa.statespace import initialization

class CAPM(sm.tsa.statespace.MLEModel):
def __init__(self, endog, exog):
super(CAPM, self).__init__(endog, exog=exog, k_states=8, k_posdef=4)
N, self.m = endog.shape
N, self.n = exog.shape
F = np.zeros((1, *exog.T.shape))
F[0] = exog.T
self['design']     = np.kron(F.transpose(2,0,1),np.eye(self.m)).transpose(1,2,0)                    # F
self['transition'] = np.kron(np.eye(self.n), np.eye(self.m))   # G
# selection matrix should make sure state_cov is not singular for numerical stability
self['selection', 4:,:]  = np.eye(4)   # R
# now initialize all the k_states, const with 'diffuse' and dynamics with 'stationary'
init = initialization.Initialization(self.k_states)
init.set((0,4), 'diffuse')
init.set((4, 8), 'stationary')
self.ssm.initialize(init)

@property
def param_names(self):
names = []
for i in range(self.m):
for j in range(self.m):
if i>=j:
names.append('var.V(%s,%s)'%(i,j))
for i in range(self.m):
for j in range(self.m):
if i>= j:
names.append('var.W(%s,%s)'%(i,j))
return names

@property
def start_params(self):
return 10*np.ones(len(self.param_names))

def fill_lower_diag(self, a):
n = int(np.sqrt(len(a) * 2))
out = np.zeros((n, n), dtype=int)
return out

def update(self, params, **kwargs):
params = super(CAPM, self).update(params, **kwargs)
Lv = self.fill_lower_diag(params[:len(params)//2])
Lw = self.fill_lower_diag(params[len(params)//2:])
self['obs_cov']   = Lv @ Lv.T  # V
self['state_cov'] = Lw @ Lw.T  # W

n = len(x)
mod = CAPM(y.values, exog=exog.values)
preliminary = mod.fit(maxiter=1000)
res = mod.fit(preliminary.params, method='nm', disp=0, maxiter=1000)


Here is some data I am using and OLS solution to get an idea of what might be a reasonable answer.

    import pandas_datareader as pdr

df = {}
for item in ['XOM', 'IBM', 'WY', 'C', '^TYX', '^NYA']:
dff = pdr.get_data_yahoo(item, '19780101', '19871201')
_r = 1 + pd.DataFrame(df).ffill().pct_change().dropna()
r = _r.groupby([_r.index.year, _r.index.month]).prod() - 1
y = r[['IBM', 'XOM', 'WY', 'C']].sub(r['^TYX'], axis=0)
x = r['^NYA'] - r['^TYX']
# OLS model ==============
mod = sm.OLS(y, exog, intercept=True)
res = mod.fit()
res.params
$$$$


This is a nice example that shows a number of potential issues when using state space models, so there are a number of suggestions that I can make that might be helpful. First, I'll describe what I think is the solution here, and then I'll describe a few more things that may be helpful.

Parameter estimation instability

The main issue you brought up was instability in the estimated parameters. The main problem that is leading to this instability is that by default, parameters are unconstrained. As a result, the optimizer can suggest $$\phi_1$$ and $$\phi_2$$ parameters that lead to a non-stationary process for $$Y_t^{(g)}$$. This is problematic for identification because it leads to two unit root processes, and it makes it difficult for the optimizer to determine which direction to go. The solution here is to constrain the parameters so that the implied processes for $$Y_t^{(g)}$$ is stationary.

A second related issue is that you are estimating standard deviations but you are not constraining them to be positive. This is another identification issue that can cause problems for the optimizer, since, for example, the $$\sigma_\epsilon = 0.01$$ produces the same loglikelihood value as $$\sigma_\epsilon = -0.01$$. A solution here is to constrain the standard deviation parameters to be positive.

In Statsmodels, these constraints are imposed using parameter transformations. In practice, you define two new functions transform_params, which converts from unconstrained parameters used by the optimizer to parameters used in the model (i.e. the constrained parameters would imply the stationary process and would have positive standard deviations), and untransform_params, which does the inverse operation.

When I use the following code, all of your examples produce estimates that are very similar to the ones from your reference.

First, we need to import two new functions:

from statsmodels.tsa.statespace.tools import (
constrain_stationary_univariate,
unconstrain_stationary_univariate)


def transform_params(self, unconstrained):
constrained = np.zeros_like(unconstrained)
constrained[:2] = constrain_stationary_univariate(unconstrained[:2])
constrained[2:] = unconstrained[2:]**2
return constrained

def untransform_params(self, constrained):
unconstrained = np.zeros_like(constrained)
unconstrained[:2] = unconstrain_stationary_univariate(constrained[:2])
unconstrained[2:] = constrained[2:]**0.5
return unconstrained


Once you use parameter transformations, then you also have to modify the the update method as follows to make sure they get applied:

def update(self, params, **kwargs):
params = super().update(params, **kwargs)
# ...


Handling singular covariance matrices

In your model, both the observation and state innovations have singular covariance matrices. I have two suggestions:

First, there is no problem when the observation innovation covariance matrix is singular, i.e. $$Var(v_t) = 0$$, so you can simply remove the following line:

self['obs_cov', 0, 0] = np.array([[1e-8]])  # V


Second, it is better not to define a singular state innovation covariance matrix, but instead to use the selection matrix, because this is more numerically stable. For example, you can rewrite your __init__ function as follows:

def __init__(self, endog):
super(GDPGAP, self).__init__(endog, k_states=4, k_posdef=3)
self['design'] = np.array([[1,0,1,0]])  # F
self['transition'] = np.array([[1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,1,0]])  # G
self['selection', :3, :3] = np.eye(3)  # R=1
self['state_cov'] = np.diag([1, 1, 1]) # W


The k_posdef argument defines the dimension of the state innovation with a positive definite (hence "posdef") covariance matrix, and then you use the selection matrix to select the appropriate innovation for each element of the state vector. Then your last line of the update method would be:

self['state_cov'] = np.diag([params[2]**2, params[3]**2, params[4]**2])


In my tests, this also improved the estimated parameters.

Initialization

Currently, you are initializing the state space model to be approximately diffuse. For a few years now, Statsmodels has had exact diffuse initialization, so it would be better to use that. You could do this by removing the both the initialize_approximate_diffuse line and the loglikelihood_burn line.

However, we can do even better. Because the first two elements of the state vector are non-stationary, a diffuse initialization is appropriate. But because the last two elements are stationary (which we have now imposed via the parameter transformations), we have the option of initializing using the unconditional distribution.

This can be done as follows. First, import the initialization package:

from statsmodels.tsa.statespace import initialization


Then, at the end of your __init__ method (where you deleted the initialize_approximate_diffuse and loglikelihood_burn lines), use:

    init = initialization.Initialization(self.k_states)
init.set((0, 2), 'diffuse')
init.set((2, 4), 'stationary')
self.ssm.initialize(init)


Working harder to find the maximum likelihood estimates

After all this, the parameters estimated by Statsmodels still do not generate a likelihood as high as the parameters in from your reference. What this suggests is that we may need to work a little harder to find the maximum likelihood estimates. One way to do this would be to try alternative starting parameters, but another way is to try to chain optimizers to get a better fit.

We can do this as follows:

# As a preliminary fit, the default method uses BFGS
preliminary = ggmodel.fit(maxiter=1000)

# Next we'll use Nelder-Mead, starting at the parameter
# estimates from the preliminary estimation
res = ggmodel.fit(preliminary.params, method='nm', disp=0, maxiter=1000)


This finally consistently generates parameter estimates that are actually a little better (generate a higher likelihood) than those from your reference text.

Full code:

Here's the full code resulting from all of the above steps:

import pandas_datareader as pdr
dff = pdr.get_data_fred('GDPC1', '19500101', '20041001').apply(np.log)

from statsmodels.tsa.statespace.tools import (
constrain_stationary_univariate,
unconstrain_stationary_univariate)

from statsmodels.tsa.statespace import initialization

class GDPGAP(sm.tsa.statespace.MLEModel):
start_params = [.1, .1, 100, 100, 100]
def __init__(self, endog):
super(GDPGAP, self).__init__(endog, k_states=4, k_posdef=3)
self['design'] = np.array([[1,0,1,0]])  # F
self['transition'] = np.array([[1,1,0,0],[0,1,0,0],[0,0,1,1],[0,0,1,0]])  # G
self['selection', :3, :3] = np.eye(3)  # R=1
self['state_cov'] = np.diag([1, 1, 1]) # W
init = initialization.Initialization(self.k_states)
init.set((0, 2), 'diffuse')
init.set((2, 4), 'stationary')
self.ssm.initialize(init)

def transform_params(self, unconstrained):
constrained = np.zeros_like(unconstrained)
constrained[:2] = constrain_stationary_univariate(unconstrained[:2])
constrained[2:] = unconstrained[2:]**2
return constrained

def untransform_params(self, constrained):
unconstrained = np.zeros_like(constrained)
unconstrained[:2] = unconstrain_stationary_univariate(constrained[:2])
unconstrained[2:] = constrained[2:]**0.5
return unconstrained

def update(self, params, **kwargs):
params = super().update(params, **kwargs)
self['transition', 2,2] = params[0]
self['transition', 3,2] = params[1]
self['state_cov'] = np.diag([params[2]**2, params[3]**2, params[4]**2])  # W

ggmodel = GDPGAP(dff)
preliminary = ggmodel.fit(maxiter=1000)
res = ggmodel.fit(preliminary.params, method='nm', disp=0, maxiter=1000)

• This is massively appreciated @cfulton. This answer cleared up a lot of things for me and taught me a lot of intricacies of a very useful package. Thanks again! Mar 27 at 3:11
• How would one introduce constraints on off diagonal covariance parameters? Apr 11 at 1:09
• Typically you would parameterize the model in terms of the Cholesky factor. For example, if you had a 2x2 matrix, you would have 3 parameters that would define the lower triangle of a matrix L, and then you would set self['state_cov'] = L @ L.T. For this case, you don't need to do anything to those specific parameters in the transform_params and untranform_params` methods. Apr 12 at 1:13
• Following up on the Cholesky factor suggestion I tried the following and it does not seem to work numerically, i.e. does not converge to a solution. Putting it as an EDIT. Apr 28 at 11:59