# Choleskly constraint in mlemodel in statsmodel

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))
out = np.zeros((n, n), dtype=int)
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.

df = {}
for item in ['XOM', 'IBM', 'WY', 'C', '^TYX', '^NYA']:
dff = pdr.get_data_yahoo(item, '19780101', '19871201')
_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']
# OLS model ==============
mod = sm.OLS(y, exog, intercept=False)
res = mod.fit()
res.params
$$`$$
• 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? May 10 at 19:43
• @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. May 11 at 2:30

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
• I can't express my gratitude enough. thanks! May 11 at 17:10
• No problem, glad it was useful May 12 at 4:09