State space models in python statsmodels: including lag state innovation in observation equation I am trying to fit the following state space model via the (excellent and highly useful) statsmodels state-space module.  The model is a standard local level model but where last period's innovation in the state equation loads into this period's observation equation via a coefficient to be estimated (a simpler version of Menkveld, Koopman and Lucas, JBES (2012) Eq (7) page 215):
$$
\alpha_{t+1} = \alpha_t + \eta_t \\
y_t = \alpha_t + \varepsilon_t + \theta \eta_{t-1}
$$
In state space form, the model can be written as:
$$
\mathbf{\delta_{t+1}} = \begin{pmatrix}  \alpha_{t+1} \\\eta_t \end{pmatrix} = 
\begin{pmatrix} 1,0 \\ 0,0 \end{pmatrix} \begin{pmatrix}  \alpha_{t} \\\eta_{t-1} \end{pmatrix} + 
\begin{pmatrix} 1 \\ 1 \end{pmatrix} \eta_t 
\\
y_t =  \begin{pmatrix} 1,  \theta\end{pmatrix} \delta_t + \varepsilon_t
$$
where $\delta_{t}$ is the state variable $(\alpha_t, \eta_{t-1})'$.
I've simulated the DGP as follows:
import numpy as np

T, k_obs = 5000, 1 # Sample length (rows of y) and number of variables (columns of y)
var_eta, var_eps = 2.8, 0.5  # Variance of the state disturbance and the obs. disturbance

rs = np.random.RandomState(seed=2017)
eta = rs.normal(scale=var_eta ** 0.5, size=T) # Sequence of state-varaible shocks  
epsil = rs.normal(scale=var_eps ** 0.5, size=T) # Sequence of obs-varaible shocks  

theta = -0.3 # Coefficient on lag state error in observation equation 

alpha, y  = np.zeros(T), np.zeros((T, k_obs)) 

# Update the first two observations of y and alpha:
y[0] = alpha[0] + epsil[0]
alpha[1] = alpha[0] + eta[0]
y[1] = alpha[1] + epsil[1] + theta*eta[0]

# Update remaining obervations as per the DGP
for t in range(2,T):
    alpha[t] = alpha[t-1] + eta[t-1] 
    y[t] = alpha[t] + epsil[t] + theta*eta[t-1]

Then I've specified and estimated a state space model as follows:
import statsmodels.api as sm

class llmodel_with_lag(sm.tsa.statespace.MLEModel):
    def __init__(self, endog):
        super(llmodel_with_lag, self).__init__(endog, 
                                              k_states=2, 
                                              k_posdef=1,
                                             initialization = 'diffuse')

        self['design', 0, 0] = 1.0 # Design matrix is a 1x2: [1,theta] 
        self['design', 0, 1] = 1.0 # Theta coefficient on lag state innovation in obs. equation
        
        self['transition', 0, 0] = 1.0 # Transition matrix is [[1,0],[0,0]]
        self['transition', 0, 1] = 0.0 
        self['transition', 1, 0] = 0.0 
        self['transition', 1, 1] = 0.0 
        
        self['selection', 0] = 1.0 # Selection matrix is [1,1]          
        self['selection', 1] = 1.0      

        self['state_cov'] = [0]
        self['obs_cov'] = [0]

        self.positive_parameters = slice(0,2)         
            
    @property
    def param_names(self):
        param_names =  [
            'state_var', 
            'obs_var',
            'design_ma'
        ]
        return param_names
    
    @property
    def start_params(self):
        params = [
            0.1,
            0.1,
            0
        ]
        return params
    
    # Constrain the variance parameters to be positive,
    def transform_params(self, unconstrained):
        constrained = unconstrained.copy()
        constrained[self.positive_parameters] = (
            constrained[self.positive_parameters] ** 2
        )
        return constrained

    def untransform_params(self, constrained):
        unconstrained = constrained.copy()
        unconstrained[self.positive_parameters] = (
            unconstrained[self.positive_parameters] ** 0.5
        )
        return unconstrained
    
    def update(self, params, **kwargs):
        params = super(llmodel_with_lag, self).update(params, **kwargs)
        
        self['state_cov', 0, 0] = params[0] # Variance of the state innov.        
        self['obs_cov', 0, 0] = params[1] # Variance of the obs. innov      
        self['design', 0, 1] = params[2] # Coefficient on lag state innov in obs. eqn

# Estimation
mod = llmodel_with_lag(y) 
modinit = mod.fit(maxiter = 1000)
modfit = mod.fit(start_params = modinit.params,
                 method = 'nm', 
                 disp = 0,
                 maxiter = 1000) 


When doing this, I am repeatedly getting a very strange result where the coefficients on the lagged state innovation $\theta$ in the observation equation, and the variance of the observation shock, are not estimated correctly, but more curiously, the $\theta$ coefficient always has a z-statistic of exactly 1.000 or -1.000 (so the magnitude of the std. error is equal to the magnitude coefficient).  I've tried this for a number of different parameter values and estimation methods and continue to generate the same z-value of 1 or -1 on $\theta$. The model converges and the condition number on the parameter covariance matrix is finite. I also obtained a similar result using real data so I don't believe the issue is in the DGP itself.
I am not sure what is driving this, but since the coefficient and standard error have virtually identical magnitudes and the condition numbers are fine, this didn't appear to be a problem that the parameters are not identified or that the model being incorrectly specified (though I could be mistaken about this).
Any ideas or help regarding what I am doing wrong or what should be changed to get the estimation working would be much appreciated!

The problem here is that the $\theta$ and $\sigma_\varepsilon^2$ parameters are not separately identified. In particular, if you work through the Kalman filtering equations, it can be shown that the loglikelihood depends on these two parameters essentially only through the sum $(\theta + \theta^2) \sigma_\eta^2 + \sigma_\varepsilon^2$. (I say "essentially", because depending on what "initialization" method is used, the loglikelihood for the first few periods can depend on the parameters in other ways, but this effect disappears quickly - in the example below, it disappears after 2 observations).
In other words, suppose that the true parameters used to generate the data are $(\bar \sigma_\eta^2, \bar \sigma_\varepsilon^2, \bar \theta)$. Then any other $\sigma_\varepsilon^2, \theta$ pair such that
$(\theta + \theta^2) \bar \sigma_\eta^2 + \sigma_\varepsilon^2 = (\bar \theta + \bar \theta^2) \bar \sigma_\eta^2 + \bar \sigma_\varepsilon^2$
will yield the same loglikelihood as the true parameters.
So, for example, if I am given a particular $\sigma_\varepsilon^2$, then I could solve for the associated $\theta$ by solving the following quadratic equation
$$\theta^2 + \theta - \frac{1}{\bar \sigma_\eta^2} \left ( (\bar \theta + \bar \theta^2) \bar \sigma_\eta^2 + \bar \sigma_\varepsilon^2 - \sigma_\varepsilon^2 \right ) = 0$$
The code snippet below shows how this works in practice:
# Construct the model
mod = llmodel_with_lag(y)
# Ignore the loglikelihood associated with the
# first two periods
mod.loglikelihood_burn = 2

# Collect the true parameters
bar_sigma2_eta = 2.8
bar_sigma2_eps = 0.5
bar_theta = -0.3

# Compute the true value of the key sum that
# appears in the loglikelihood
key_sum = (bar_theta + bar_theta**2) * bar_sigma2_eta + bar_sigma2_eps

# Specify an alternative \sigma_\varepsilon^2,
sigma2_eps = 0.3

# Compute the corresponding alternative \theta
# by solving the quadratic equation
a = 1
b = 1
c = -(key_sum - sigma2_eps) / bar_sigma2_eta

theta_option1 = (-b + (b**2 - 4 * a * c)**0.5) / (2 * a)
theta_option2 = (-b - (b**2 - 4 * a * c)**0.5) / (2 * a)

# Check that these all yield the same loglikelihood
print('Loglikelihood at true parameters')
print(mod.loglike([bar_sigma2_eta, bar_sigma2_eps, bar_theta]))

print('Loglikelihood at sigma2_eps, theta_option1')
print(mod.loglike([bar_sigma2_eta, sigma2_eps, theta_option1]))

print('Loglikelihood at sigma2_eps, theta_option2')
print(mod.loglike([bar_sigma2_eta, sigma2_eps, theta_option2]))

Which prints out:
Loglikelihood at true parameters
-9394.989471741343
Loglikelihood at sigma2_eps, theta_option1
-9394.989471741343
Loglikelihood at sigma2_eps, theta_option2
-9394.989471741343

