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
```