Choleskly constraint in mlemodel in statsmodel

Question

I want to constraint the off diagonal terms in the covariance matrix in a dynamic linear model. I tried using Cholesky method but it does not seem to converge.

I am trying to fit a multivariate CAPM model on 4 assets (XOM, IBM, WY, C) with static alpha and dynamic beta. The observation and state covariance matrix are the variables we want to estimate and alpha and beta are the states. To make the covariance matrix positive definite, i use Cholesky construction $\Sigma = LL^T$. However the loglikelihood does not seem to converge when i do that.

For each asset we can write the DLM \begin{align*} y_{it}&=\begin{bmatrix}1 & x_t\end{bmatrix}\begin{bmatrix}\alpha_{it}\\\beta_{it}\end{bmatrix} + v_{it}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ v_{it}\sim\mathcal N(0, V_i) \\ \begin{bmatrix}\alpha_{it} \\ \beta_{it}\end{bmatrix} &= \begin{bmatrix}1&0\\0&1\end{bmatrix} \begin{bmatrix}\alpha_{i,t-1} \\ \beta_{i,t-1}\end{bmatrix} + \begin{bmatrix}w_{1,it}&0\\0&w_{2,it}\end{bmatrix} \ \ \ \ \ \ \ w_{1,it}\sim\mathcal N(0, W_{1i}), w_{2,it}\sim\mathcal N(0, W_{2i}) \end{align*} We now assume that for the four assets, the intercepts and slopes are correlated. We can write the combined DLM for the 4 stocks as \begin{align*} \underset{(4 \times 1)}{y_t} &= \underset{(4 \times 8)}{(F \otimes I_4)}\underset{(8 \times 1)}{\theta_t} + \underset{(4 \times 1)}{v_t}, \ \ \ & v_t \sim \mathcal N(0, \underset{(4\times 4)}{V}) \\ \underset{(8\times 1)}{\theta_t} &= \underset{(8\times 8)}{(G \otimes I_4)}\underset{(8\times 1)}{\theta_{t-1}} + \underset{(8\times 1)}{w_t}, \ \ \ & w_t \sim \mathcal N(0, \underset{(8\times 8)}{W}) \end{align*} with $$y_t = \begin{bmatrix}y_{1t}\\ \vdots \\ y_{4t}\end{bmatrix}, \theta_t = \begin{bmatrix}\alpha_{1t}\\ \vdots \\ \alpha_{4t} \\ \beta_{1t} \\ \vdots \\ \beta_{4t}\end{bmatrix}, v_t = \begin{bmatrix}v_{1t} \\ \vdots \\ v_{4t}\end{bmatrix}, w_t =\begin{bmatrix}w_{\alpha 1t}\\ \vdots \\ w_{\alpha 4t} \\ w_{\beta 1t} \\ \vdots \\ w_{\beta 4t}\end{bmatrix}, \\ F=\begin{bmatrix}1 & x_t \end{bmatrix}, G=\begin{bmatrix}1&0\\0&1\end{bmatrix}, W=\begin{bmatrix}\underset{(4\times 4)}{W_{\alpha}}&0\\0&\underset{(4\times 4)}{W_{\beta}}\end{bmatrix}$$ We assume that $\alpha_{it}$ are time-invariant, i.e. $W_{\alpha}=0$. The correlation between the different excess returns is explained in terms of the non-diagonal variance matrices $V$ and $W_{\beta}$.

Here is the code that encodes the model using the tsa.statespace module.

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))
            mask = np.tri(n, dtype=bool, k=0)
            out = np.zeros((n, n), dtype=int)
            out[mask] = a
            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 I show the dimensions of the various matrices as a check

mod.param_names
Out[140]: 
['var.V(0,0)',
 'var.V(1,0)',
 'var.V(1,1)',
 'var.V(2,0)',
 'var.V(2,1)',
 'var.V(2,2)',
 'var.V(3,0)',
 'var.V(3,1)',
 'var.V(3,2)',
 'var.V(3,3)',
 'var.W(0,0)',
 'var.W(1,0)',
 'var.W(1,1)',
 'var.W(2,0)',
 'var.W(2,1)',
 'var.W(2,2)',
 'var.W(3,0)',
 'var.W(3,1)',
 'var.W(3,2)',
 'var.W(3,3)']
mod['design'].shape
Out[141]: (4, 8, 120)
mod['transition'].shape
Out[142]: (8, 8)
mod['selection'].shape
Out[143]: (8, 4)
mod['obs_cov'].shape
Out[144]: (4, 4)
mod['state_cov'].shape
Out[145]: (4, 4)

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')
        df[item] = dff['Adj Close']
    _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']
    exog = sm.add_constant(x)
    # OLS model ==============
    mod = sm.OLS(y, exog, intercept=False)
    res = mod.fit()
    res.params
```

Answer 1

I think there are two main issues here, and then a third minor thing:

The first is that transition equation implies a random walk for every element of the state vector, so you can't use 'stationary' initialization for the last four elements. You should use:

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

The second is that in fill_lower_diag, you're currently forcing the data into the integer datatype, and when you do that you lose all of the floating point resolution, so it should be:

out = np.zeros((n, n), dtype=a.dtype)

When I make these two changes, the optimizer converges to a solution that's pretty similar to the OLS case.

The third thing you might consider is improving the start_params a bit. The ones you're using imply very large variances, which doesn't really fit this setup. In particular, good starting variances for the time-varying regression parameters tend to be pretty small (i.e. close to the baseline scenario where the variance is zero, which corresponds to the usual regression estimates). So I would recommend something like:

@property
def start_params(self):
    V = np.diag(self.endog.var(axis=0))
    W = np.eye(self.ssm.k_posdef) * 1e-1
    
    ix_V = np.tril_indices_from(V)
    ix_W = np.tril_indices_from(W)

    params = np.r_[
        np.linalg.cholesky(V)[ix_V],
        np.linalg.cholesky(W)[ix_W]
    ]

    return params