Including parameters in state space model from statsmodels

Building up the model from a previous post, and the helpful answer, I've subclassed the MLEModel to encapsulate the model. I'd like to allow for two parameters q1 and q2 so that the state noise covariance matrix is generalized as in Sarkka (2013)'s example 4.3 (terms re-arranged for my convention):

I thought I would accomplish this with the update method below, but I'm running into problems with the fit method, as it returns a UFuncTypeError: Cannot cast ufunc 'multiply' output from dtype('complex128') to dtype('float64') with casting rule 'same_kind'. What am I missing here?

import numpy as np
import scipy.linalg as linalg
import statsmodels.api as sm

class Tracker2D(sm.tsa.statespace.MLEModel):
    """Position tracker in two dimensions with four states

    """
    start_params = [1.0, 1.0]
    param_names = ["q1", "q2"]

    def __init__(self, endog):
        super(Tracker2D, self).__init__(endog, k_states=4)
        self.endog = endog
        self._state_names = ["x1", "dx1/dt",
                             "x3", "dx3/dt"]
        # dt: sampling rate; s = standard deviation of the process noise
        # common to both dimensions
        dt, s = 0.1, 0.5
        # dynamic model matrices A and Q
        A2d = [[1, dt],
               [0, 1]]
        A = linalg.block_diag(A2d, A2d)
        Q2d = [[dt ** 3 / 3, dt ** 2 / 2],
               [dt ** 2 / 2, dt]]
        Q = linalg.block_diag(Q2d, Q2d)
        # measurement model matrices H and R
        H = np.array([[1, 0, 0, 0],
                      [0, 0, 1, 0]])
        R = s ** 2 * np.eye(2)
        self["design"] = H
        self["obs_cov"] = R
        self["transition"] = A
        self["selection"] = np.eye(4)
        self["state_cov"] = Q

    def update(self, params, **kwargs):
        self["state_cov", :2, :2] *= params[0]
        self["state_cov", 2:, 2:] *= params[1]


# Initialization
m0 = np.array([[0, 1, 0, -1]]).T  # state vector column vector
P0 = np.eye(4)                   # process covariance matrix

# With object Y below being the simulated measurements in downloadable
# data file from previous post
with open("measurements_2d.npy", "rb") as f:
    Y = np.load(f)

tracker2D = Tracker2D(pd.DataFrame(Y.T))
tracker2D.initialize_known((tracker2D["transition"] @ m0.flatten()),
                           (tracker2D["transition"] @ P0 @
                            tracker2D["transition"].T +
                            tracker2D["state_cov"]))
# Below throws the error
tracker2D.fit()

The error message you are receiving is about trying to set a complex value in a dtype=float matrix. You would get the same error from:

A = np.eye(2)
A *= 1.0j

The error is showing up in:

def update(self, params, **kwargs):
    self["state_cov", :2, :2] *= params[0]
    self["state_cov", 2:, 2:] *= params[1]

because you are modifying the "state_cov" in place. When params is a complex vector but the existing "state_cov" matrix has dtype float, then the error will occur. Statsmodels will set the parameter vector to be complex when computing the standard errors of the parameters, because it uses complex step differentiation.

You could use something like

def update(self, params, **kwargs):
    self["state_cov", :2, :2] = params[0] * self["state_cov", :2, :2]
    self["state_cov", 2:, 2:] = params[1] * self["state_cov", 2:, 2:]

Although I should point out that this will not give you what I think you want, because it will modify the "state_cov" based on whatever it previously was. I think instead, you want something like:

class Tracker2D(sm.tsa.statespace.MLEModel):
    """Position tracker in two dimensions with four states

    """
    start_params = [1.0, 1.0]
    param_names = ["q1", "q2"]

    def __init__(self, endog):
        super(Tracker2D, self).__init__(endog, k_states=4)
        self.endog = endog
        self._state_names = ["x1", "dx1/dt",
                             "x3", "dx3/dt"]
        # dt: sampling rate; s = standard deviation of the process noise
        # common to both dimensions
        dt, s = 0.1, 0.5
        # dynamic model matrices A and Q
        A2d = [[1, dt],
               [0, 1]]
        A = linalg.block_diag(A2d, A2d)
        Q2d = [[dt ** 3 / 3, dt ** 2 / 2],
               [dt ** 2 / 2, dt]]

        # First we save the base Q matrix
        self.Q = linalg.block_diag(Q2d, Q2d)

        # measurement model matrices H and R
        H = np.array([[1, 0, 0, 0],
                      [0, 0, 1, 0]])
        R = s ** 2 * np.eye(2)
        self["design"] = H
        self["obs_cov"] = R
        self["transition"] = A
        self["selection"] = np.eye(4)
        self["state_cov"] = self.Q.copy()

    def update(self, params, **kwargs):
        # Now update the state cov based on the original Q
        # matrix, and set entire blocks of the matrix, rather
        # than modifying it in-place.
        self["state_cov", :2, :2] = params[0] * self.Q[:2, :2]
        self["state_cov", 2:, 2:] = params[1] * self.Q[2:, 2:]