# A maxdet problem with three variables - maximum likelihood estimation in Mathlab/Mata/R

## Problem

Given $n$ known vectors $y_j = (y_1, y_2, ..., y_{K_j})', \forall j=1,n$ and a constant $K$; determine three vectors ($\alpha, \beta, \theta^2)$ that maximise:

$$\text{T} (\alpha, \beta, \theta^2) = \sum\limits_{j=1}^{n} [\log(\det \Omega) + (y_j - \alpha)'\Omega^{-1}(y_j - \alpha)]$$

where $\alpha = (\alpha_1, \alpha_2, ..., \alpha_{K_j})'$, $\beta = (\beta_1, \beta_2, ..., \beta_{K_j})'$ and $\theta^2 = (\theta^2_1, \theta^2_2, ..., \theta^2_{K_j})'$; $\text{B}$, $\Theta$ and $\Omega$ are three $K_j \times K_j$ matrices as: $$\text{B} = \begin{bmatrix} \beta_1 & 0 & ... & 0 & 0 \\ 0 & \beta_2 & ... & 0 & 0 \\ ... & ... & ... & ... & ... \\ 0 & 0 & ... & 0 & \beta_{K_j} \end{bmatrix}$$ $$\Theta = \begin{bmatrix} \theta^2_1 & 0 & ... & 0 & 0 \\ 0 & \theta^2_2 & ... & 0 & 0 \\ ... & ... & ... & ... & ... \\ 0 & 0 & ... & 0 & \theta^2_{K_j} \end{bmatrix}$$ $$\Omega = \beta \beta' + B \Theta B'$$

## Application

This maximisation is used to solve the Unobserved Component Models. I happen to know that if $\beta_i \beta_j = 0$ then Stata can run GLS to maximise $\text{T}$ given the $n$ vectors $y_j$.

In a general case, I am not aware of any mathematical/statistical programme that can yield/perform the algorithm to maximise $\text{T}$. If anyone happens to come across this type of problem, please let me know too.

I am from economics major and just want to understand the algorithms behind the maximum-likelihood estimators. I find this problem from a working paper of the World Bank and attempt to tackle.

Thank you!

• You haven't written this in a way which makes sense. Do you have a link to the actual paper? You've written B and Θ matrices as being diagonal, which also makes ′Ω diagonal. I presume you meant B * B' in the definition of Ω, not ββ′. If all the β 's and θ 's are positive, the second term in the sum will be convex. The first term in the sum is concave. If this is truly the function, it's a "nasty" non-convex problem.It may be solvable with general purpose local optimizer, and maybe with global optimizer, but not necessarily solvable to completion. Ability to solve can depend on starting value. – Mark L. Stone Jun 24 '15 at 8:55
• If those matrices really are diagonal, then multiplying them by their transpose, as written, just amounts to squaring each entry on the diagonal. This is part of why it looks suspicious that you didn't write this out correctly. – Mark L. Stone Jun 24 '15 at 8:58
• Are there constraints, such as nonnegativity, on the vectors which are being optimized? Or maybe constraints on matrices being positive definite or positive semi-definite (distinction only matters if those matrices really aren't diagonal). – Mark L. Stone Jun 24 '15 at 9:09
• Edit to my first comment: ββ' could make sense. That's a rank 1 outer product, and must be positive semi-definite, and if j > 1, will be singular. – Mark L. Stone Jun 24 '15 at 9:16
• Hi Mark, Thank you for your attention and my apology for the confusion. In the definition of $\Omega$, it is $\beta \beta'$, that is a symmetric matrices not diagonal (that is why in the application, I am able to understand the GLS command in Stata when $\beta_i \beta_j = 0$). I think we can reasonably assume $\theta > 0$ because the problem state to find $\beta^2$; but there is no constraint for $\beta$. The paper I mentioned can be downloaded/view from papers.ssrn.com/sol3/papers.cfm?abstract_id=1148386 Please refer to the Appendix, page 97-99. – Khan Jun 24 '15 at 9:25

I think you need nonnegativity constraints on all the parameters (variables). MAXDET by Boyd et al, and its later more flexible evolution CVX http://cvxr.com/cvx/ , can NOT be used because even though the log determinant term is concave, the other term is convex, hence, you have a "nasty" non-convex problem. There could be multiple local maxima, and there may even be local minima and saddle points. Different starting points (initial values (guesses) for the parameters could result in different solutions.

You can not maximize one term at a time. It will not work.

Forget about maximum likelihood estimation software, and use a general purpose local nonlinear optimizer. You need to be able to evaluate the objective function. That is easy to do in MATLAB. If you have the Optimization toolbox, you could use FMINCON. Because FMINCON only minimizes, you'll have to minimize the negative of the function you want to maximize. Be sure to specify the lower bound constraints. If you don't know how to evaluate the gradient, you can tell the optimizer to use finite differences for the gradient, but that will make it take longer if there are a lot of parameters to be estimated. If you can get access to KNITRO http://www.ziena.com/knitro.htm , use that instead of FMINCON, because even though the algorithms are nominally basically the same as FMINCON, it's a much higher quality implementation, which could make the difference between solving it in a minute in KNITRO vs. having FMINCON labor for hours or days and still not solve it (on the other hand, FMINCON might work as well on this problem, but likely not). You should try several different starting values.

You could also call an optimizer from C/C++, Python, Fortran, or maybe R.

Also, the question (thread) title is misleading. It is NOT a maxdet problem. If if were a maxdet of a positive definite matrix problem, with no other term in the objective function, then it would be a concave maximization, which is considered to be a convex optimization problem, which is easy and fast to formulate and solve (unless there are several thousand or more parameters being estimated) using a tool such as CVX. But because you are adding a convex function to maxdet, it is decidedly non-convex, and not so easy, and neither maxdet nor CVX can be used at all.

If you want to understand how MLE and nonlinear optimization algorithms work, and how they SHOULD work, take a course in nonlinear optimization.

• This answer has been unaccepted. Did you encounter difficulty in trying to follow my advice? – Mark L. Stone Jul 3 '15 at 17:31
• I am sorry, it was my fault (clicking around with my iPad). I've been reading your answer for a while. I am able to solve the maximisation by hand now (using matrix determinant lemma) and derivations. However you are right that I need to take further course in coding/econometrics to put those manual calculations into a statistical package. Thank you very much! – Khan Jul 4 '15 at 9:22
• Hi Mark, following up this post, some of my friends suggest package stats4 in R that can deal with the problem. I am learning the language now but just to let you know how I get on. Thank you! – Khan Jul 6 '15 at 12:43
• I think they may be talking about mle in stats4, which calls optim. You can specify bounds on the parameters being estimated, so be sure to do so. If your problem is well behaved and easy to solve, which I think it may not be, using mle with optim might be o.k., presuming you can get the likelihood correctly modeled in mle. But optim is not a very good optimizer. The best general purpose (i.e.,, not specific to mle) nonlinear optimizers are much better. You already know the objective function and constraints. You're better off avoiding mle software and using another nonlinear optimizer. – Mark L. Stone Jul 6 '15 at 12:58
• I will keep learning and let you know how I get on. Thank you for your advice!!! – Khan Jul 6 '15 at 13:03