Maximum Likelihood in a time series multi-population model Currently I am trying to implement a mortality rate model (Lee-Li), but I am having some trouble with the maximum likelihood procedure of a time series model given by:
\begin{align}
K_{t+1} &= K_t + \theta + \epsilon_{t+1}\\ 
\kappa_{t+1} &= a \kappa_{t} + \delta_{t+1}
\end{align}
Where $E_t := (\delta_t,\epsilon_t)$ are assumed to be independent and come from a bi variate normal distribution with mean $(0,0)$ and co-variance matrix $\mathbf C$. 
I have data available for both $K_t$ and $\kappa_t$, and the idea is to maximize for $a$, $\theta$ and $\mathbf C$ using maximum likelihood. However I have no clue as how to approach this problem. I tried obtaining the MLE for $\mathbf C$ using the MLE for the covariance matrix of a bivariate normal distribution, but this seems to be wrong. 
Furthermore I have no idea on how to approach $\theta$ and a for this matter either, since I don't know how to set up a (log)-likelihood for this expression.
Any help would be appreciated.
Edit: 
The co-variance matrix $\mathbf C = 
\begin{bmatrix} 
\sigma^{2}_{\epsilon} & \rho \sigma_{\epsilon} \sigma_{\delta}& \\
\rho \sigma_{\epsilon} \sigma_{\delta} & \sigma^{2}_{\delta}
\end{bmatrix}$
Where the correlation coefficient $\rho$ is non-zero.
I have tried to do the following steps to derive the maximum likelihood for $\theta$ :
1: Condition the distribution of $K_{t+1}$ on past data $t$ and the parameters
$\mathbf{\Theta}=[a,\theta,\sigma^{2}_{\epsilon},\sigma^{2}_{\delta},\rho \sigma_{\epsilon} \sigma_{\delta}]$, which has a bi variate normal distribution 
2: compute the conditional log-likelihood of this equation:
\begin{align}
\sum_{t}\log(L \left(\mathbf{\Theta} \right)) &= \log\left(\frac{1}{2 \pi \sigma_{\epsilon} \sigma_{\delta}\sqrt{ 1 - \rho^{2}}}\right) -  \frac{z}{2(1 - \rho^{2}}
\end{align}
where 
\begin{align}
z = \frac{(K_{t+1} - (K_{t} + \theta))^{2}}{\sigma_\epsilon^2} + \frac{(\kappa_{t+1} - a\kappa_{t}) ^{2}}{\sigma_\delta^2} - \frac{2\rho(K_{t+1} - (K_{t} + \theta))(\kappa_{t+1} - a\kappa_{t})}{\sigma_\epsilon \sigma_\delta}
\end{align}
3: take the derivative with respect to $\theta$. Since $z$ is the only expression involving $\theta$ we can ignore the rest:
 \begin{align}
\frac{\partial \log \left( L \right)}{\partial \theta} &= 0 \Leftrightarrow \\
\sum_{t}\frac{-2(K_{t+1} - (K_t+\theta))}{\sigma_{\epsilon}^2} + \sum_{t}\frac{2\rho(\kappa_{t+1} - a\kappa_t)}{\sigma_{\epsilon}\sigma_{\delta}}&=0 \Leftrightarrow \\
\sum_{t}\frac{K_{t+1} - (K_t+\theta)}{\sigma_{\epsilon}^2} &= \sum_{t}\frac{\rho(\kappa_{t+1} - a\kappa_t)}{\sigma_{\epsilon}\sigma_{\delta}} 
\Leftrightarrow \\
\sum_{t}K_{t+1} - K_t - n\theta &= \sigma_{\epsilon}^2\sum_{t}\frac{\rho(\kappa_{t+1} - a\kappa_t)}{\sigma_{\epsilon}\sigma_{\delta}} \Leftrightarrow \\
 \hat{\theta_{MLE}} =\sum_{t}\frac{(K_{t+1} - K_t)}{n} &-\sigma_{\epsilon}^2\sum_{t}\frac{\rho(\kappa_{t+1} - a\kappa_t)}{n\sigma_{\epsilon}\sigma_{\delta}} 
\end{align}
Similar calculations for $a,\sigma_{\epsilon},\sigma_{\delta}$ yield 
\begin{align}
\hat{a_{MLE}} &= \frac{\sum_{t} \kappa_{t+1} \kappa_{t} - \frac{\sigma_{\delta}}{\sigma_{\epsilon}} \rho \sum_{t} \kappa_{t}(K_{t+1} - K_{t} - \theta)}{ \sum_{t} \kappa_{t}^{2}} \\
\hat{\sigma_{\epsilon}^{MLE}} &= \frac{-\frac{\rho*(K_{t+1} - K_{t} - \theta)(\kappa_{t+1} - a*\kappa_{t})}{ \sigma_{\delta}} + \sqrt D_{\epsilon} )}{ 2n(1-\rho^2)} \\
\hat{\sigma_{\delta}^{MLE}} &= \frac{-\frac{\rho*(K_{t+1} - K_{t} - \theta)(\kappa_{t+1} - a*\kappa_{t})}{ \sigma_{\epsilon}} + \sqrt D_{\delta} )}{ 2n(1-\rho^2)}
\end{align}
Where 
\begin{align}
D_{\epsilon} &=  \left(\frac{\rho \sum_{t} (K_{t+1} - K_{t} - \theta)(\kappa_{t+1} - a*\kappa_{t})}{\sigma_{\epsilon}}\right)^{2} + 4 n(1-\rho^2) \sum_{t} (\kappa_{t+1} - a\kappa_{t})^2 \\
D_{\delta} &=\left(\frac{\rho \sum_{t} (K_{t+1} - K_{t} - \theta)(\kappa_{t+1} - a*\kappa_{t})}{\sigma_{\epsilon}}\right)^{2}  + 4  n(1-\rho^2) \sum_{t} (K_{t+1} - K_{t} - \theta)^2 \\
\end{align}
Finally, since the expression is too complicated analytically, one can use numerical techniques to solve the ML equations. I have used the quasi Newton-Raphson algorithm proposed by Goodman to determine the parameters.
 A: If you define $\mathbf{K}\equiv [K,\kappa]$ and $\boldsymbol{\mu}\equiv [K_t+\theta,a\kappa_t]$, then your model is
$$
\mathbf{K}_{t+1}\vert(\mathbf{K}_t,a,\theta,\mathbf{C})\sim \mathrm{N}_{\boldsymbol{\mu},\mathbf{C}}$$
i.e. the distribution of your data at $t+1$ conditioned on both the data at $t$, and on your parameters, is a bivariate normal with the specified non-zero conditional mean. The equation above is the likelihood.
If your covariance matrix were isotropic ($\mathbf{C}=\sigma^2\mathbf{I}$) this would be simple linear regression (i.e. with $\mathbf{X}=\mathbf{K}_t$ and $\mathbf{Y}=\mathbf{K}_{t+1}$). So the only change is your have a general covariance matrix, so three components (e.g. two variances and a correlation coefficient).
For MLE you want to determine the parameter values that make the gradient of the log-likelihood zero. Your parameters are $\mathbf{\Theta}=[a,\theta,C_{K,K},C_{\kappa,\kappa},C_{K,\kappa}]$, where I have expanded out the diagonal and off-diagonal components of the symmetric covariance matrix. For each parameter you get an equation by setting the gradient to zero, i.e. $\frac{\partial L}{\partial \Theta_k}=0$ for $k=1,\ldots,5$, where $L$ is the negative log likelihood.
Note that for the first two parameters $a,\theta$ you can get the gradient from the chain rule, i.e. 
$$\frac{\partial L}{\partial a}=\frac{\partial L}{\partial \mu_{\kappa}}\frac{\partial \mu_{\kappa}}{\partial a}$$ and $$\frac{\partial L}{\partial \theta}=\frac{\partial L}{\partial \mu_K}\frac{\partial \mu_K}{\partial \theta}$$
where $\frac{\partial \mu_{\kappa}}{\partial a}=\kappa_t$ and $\frac{\partial \mu_K}{\partial \theta}=1$.
