Exogenous variable in the state equation in statespace MLEmodel in statsmodels

Question

I'm trying to fit the following model:

$y_t = \left[\begin{matrix} (1-w) & 1 & w \end{matrix}\right] \left[\begin{matrix} d_t \\ \mu_t \\ m_t \end{matrix}\right] + \mathcal{N}(0,\sigma_\eta^2) $

$ \left[\begin{matrix} d_{t+1} \\ \mu_{t+1} \\ m_{t+1} \end{matrix}\right] = \left[\begin{matrix} (1-\alpha) & \alpha & 0 \\ 0 & 1 & 0 \\ 0& 0 & 1 \end{matrix}\right] \left[\begin{matrix} d_t \\ \mu_t\\ m_t \end{matrix}\right] + \mathcal{N}\left(0,\left[\begin{matrix} \sigma_\varepsilon^2 & 0 & 0 \\ 0 & \sigma_\xi^2 & 0 \\ 0 & 0 & 0\end{matrix}\right]\right) $

Where $y_t$ is a weighted average of a state variable $d_t$ and an exogenous variable $m_t$, but I have trouble putting this into the proposed custom statespace models as proposed in https://www.statsmodels.org/dev/examples/notebooks/generated/statespace_custom_models.html.

I don't understand how to implement the exogenous variable $m_t$ into the state equation. Below is the code I have got so far, where $m_t$ is not yet integrated.

class TVRegressionExtended(sm.tsa.statespace.MLEModel):
def __init__(self, y_t, m_t):
    exog = np.c_[m_t]
    
    super(TVRegressionExtended, self).__init__(
        endog=y_t, exog=exog, k_states=3,
        initialization='diffuse')
    
    #These have ok shape. Placeholders since I'm changing them
    #in the update() function
    self.ssm['design'] = np.array([1, 1, 1])
    self.ssm['selection',:2,:2] = np.eye(2)
    self.ssm['transition'] = np.array([[1, 1, 0],
                                       [0, 1, 0],
                                       [0, 0, 1]]) #{[0,0], [0,1]} = {1-alpha, alpha} 
                                                   #will be changed through update

    #Positive params
    self.positive_parameters = slice(0, 3)

@property
def param_names(self):
    return ['std.measure', 'std.return', 'std.mean',
            'alpha','weight']

@property
def start_params(self):
    params = np.r_[0.001, 0.001, 0.001, 0.8, 0.2]
    return params

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(TVRegressionExtended, self).update(params, **kwargs)
    self['obs_cov', 0, 0] = params[0]
    self['state_cov',:2,:2] = np.diag(params[1:3])
    self['transition', 0, 0] = params[3]
    self['transition', 0, 1] = 1-params[3]
    self['design', 0,0] = 1-params[4]
    self['design', 0,2] = params[4]

I think the solution is to time-vary the transition matrix in some way, but I can't seem to figure out how.

cfulton · Accepted Answer · 2021-06-16T14:56:40.477

There's no need to make the exogenous variable part of the state vector. Instead, you can write the model as:

$y_t = w m_t + \left[\begin{matrix} (1-w) & 1 \end{matrix}\right] \left[\begin{matrix} d_t \\ \mu_t \end{matrix}\right] + \mathcal{N}(0,\sigma_\eta^2) $

$ \left[\begin{matrix} d_{t+1} \\ \mu_{t+1} \end{matrix}\right] = \left[\begin{matrix} (1-\alpha) & \alpha \\ 0 & 1 \end{matrix}\right] \left[\begin{matrix} d_t \\ \mu_t \end{matrix}\right] + \mathcal{N}\left(0,\left[\begin{matrix} \sigma_\varepsilon^2 & 0 \\ 0 & \sigma_\xi^2\end{matrix}\right]\right) $

and put the $w m_t$ component into the observation intercept (self['obs_intercept']). For example:

Edit: this has now been updated to set k_exog and add a clone method, so that forecasting works.

class TVRegressionExtended(sm.tsa.statespace.MLEModel):
    def __init__(self, y_t, exog):
        exog = np.c_[exog]
        self.k_exog = 1
        
        super(TVRegressionExtended, self).__init__(
            endog=y_t, exog=exog, k_states=2,
            initialization='diffuse')
        
        #These have ok shape. Placeholders since I'm changing them
        #in the update() function
        self.ssm['design'] = np.array([1, 1])
        self.ssm['selection',:2,:2] = np.eye(2)
        self.ssm['transition'] = np.array([[1, 1],
                                           [0, 1]])

        #Positive params
        self.positive_parameters = slice(0, 3)
        
    def clone(self, endog, exog=None, **kwargs):
        return self._clone_from_init_kwds(endog, exog=exog, **kwargs)

    @property
    def param_names(self):
        return ['std.measure', 'std.return', 'std.mean',
                'alpha','weight']

    @property
    def start_params(self):
        params = np.r_[0.001, 0.001, 0.001, 0.8, 0.2]
        return params

    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(TVRegressionExtended, self).update(params, **kwargs)

        self['obs_cov', 0, 0] = params[0]
        self['obs_intercept'] = np.reshape(params[4] * self.exog, (1, self.nobs))

        self['state_cov',:2,:2] = np.diag(params[1:3])
        self['transition', 0, 0] = 1 - params[3]
        self['transition', 0, 1] = params[3]
        self['design', 0,0] = 1 - params[4]

Thanks for the quick response, very helpful. I was also wondering how I would go about forecasting? I read in another question that I would have to add the clone function, which I did. I split my exogenous variables into a training and test set. Then after fitting the model on the training set, I presumed using `res.get_forecast(m_t[n_train:])` but results in AttributeError: 'TVRegressionExtended' object has no attribute 'k_exog'. Any chance you could help me out here? — Thijs de Jongh, Jun 15 '21 at 10:20
You're on the right track, you basically just need to set the `k_exog` attribute in the model constructor. I've edited the answer to show the revised class that allows forecasting. — cfulton, Jun 16 '21 at 14:57
Again, many thanks for helping out! I have one more question. How would I go about incorporating stochastic volatility into the model? e.g. if I extend your model, by replacing $\mathcal{N}(0, \sigma_\eta^2)$ with $\sigma_\eta exp(\frac12 \theta_t) \eta_t$, where $\eta_t$ is now std. normal and $\theta_t$ the unobserved log-volatility. An extra unobserved equation needs to be added: $\theta_{t+1} = \theta_{t} + \psi_t$. But how do I go forward from here - how would I add $\theta_t$ into the update() function? — Thijs de Jongh, Jun 29 '21 at 21:00
The stochastic volatility cannot be fit as part of a linear Gaussian state space model. Often these models are fit using Bayesian Gibbs sampling techniques, which essentially decouple the problem into two separate state space models: (1) the original model, conditioning on the stochastic volatility terms, and (2) the random walk model for the stochastic volatility terms, conditioning on the states from the original model. — cfulton, Jul 01 '21 at 00:58

Exogenous variable in the state equation in statespace MLEmodel in statsmodels

1 Answers1