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')

        def param_names(self):
            names = []
            for i in range(self.m):
                for j in range(self.m):
                    if i>=j:
            for i in range(self.m):
                for j in range(self.m):
                    if i>= j:
            return names

        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

Out[141]: (4, 8, 120)
Out[142]: (8, 8)
Out[143]: (8, 4)
Out[144]: (4, 4)
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()
  • $\begingroup$ I'm having a hard time seeing what the model is supposed to be, so it would be helpful to write out the state space form you're trying to achieve. One thing that jumps out is that you are estimating covariances between all of the state innovation terms, but it seems to me like there are two separate blocks of state innovations that should be uncorrelated? $\endgroup$
    – cfulton
    Commented May 10, 2021 at 19:43
  • $\begingroup$ @cfulton I have updated the question will the details. Regarding the state covariances, I am only estimating the lower block $W_{\beta}$ in code, which is one of the blocks, while assuming the upper block $W_{\alpha}$ is zero. $\endgroup$
    – manav
    Commented May 11, 2021 at 2:30

1 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')

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:

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_[

    return params
  • $\begingroup$ I can't express my gratitude enough. thanks! $\endgroup$
    – manav
    Commented May 11, 2021 at 17:10
  • $\begingroup$ No problem, glad it was useful $\endgroup$
    – cfulton
    Commented May 12, 2021 at 4:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.