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.

  • $\begingroup$ I am glad if my answer was able to help you. However note that it is fine for you to answer your own question on StackExchange. So I might suggest that you could put your edit as an answer? (And you can accept that one, then which is certainly more complete!) $\endgroup$
    – GeoMatt22
    Commented Aug 29, 2016 at 15:29
  • $\begingroup$ I've solved it yesterday. I have used Goodman's algorithm to determine the parameters. Your suggestion to put the correlation to zero as a first guess converged pretty fast. $\endgroup$
    – Stephan
    Commented Aug 31, 2016 at 11:36

1 Answer 1


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$.

  • $\begingroup$ Sorry for the late reply, I was caught up at work. I have tried using your approach and wrote down the conditional likelihood with respect to a bi variate normal distribution. However, when I tried taking the derivative with respect to \bf{\mu} I still do not have the required parameters that I want to maximize, since I am maximizing $\bf{\mu}$ and not, for instance $\theta$ $\endgroup$
    – Stephan
    Commented Aug 28, 2016 at 15:54
  • $\begingroup$ I am writing a new comment since I had no more time to edit the old one. How do you suggest estimating these parameters. I have tried taking the derivative, for example, with respect to $\theta$ but when I plugged these in they did not seem to make any sense. I'll update my post with some more information regarding the covariance matrix and what I have tried in a short while. $\endgroup$
    – Stephan
    Commented Aug 28, 2016 at 16:03
  • $\begingroup$ I updated my answer on using the chain rule to get the derivatives w.r.t. your parameters from the derivatives w.r.t. $\mu$. $\endgroup$
    – GeoMatt22
    Commented Aug 28, 2016 at 16:17
  • $\begingroup$ @Stephan one check on your formulas is the following. If the covariance is diagonal then then $K$ and $\kappa$ are uncoupled. In this case I believe MLE should give $\theta$ as the sample mean of $K_{t+1}-K_t$, and $a$ as the projection of $\kappa_{t+1}$ onto $\kappa_t$. And the diagonal variances should be the RMS error of these predictors. So whatever formulas you get from MLE should give the same result in this case. $\endgroup$
    – GeoMatt22
    Commented Aug 28, 2016 at 17:54
  • $\begingroup$ I've been able to derive the right estimators, or at least I think. The only problem at this moment is that I can't seem to estimate anything since the equations are all dependent on each other.On top of that I can not seem to get the right correlation coefficient. Do you have any idea on how to solve this. It seems that the non negative correlation parameter is messing things up. $\endgroup$
    – Stephan
    Commented Aug 29, 2016 at 12:37

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.