Calculating Log-likelihood using Raphson and Jacobian matrices?

Question

I am reading the following paper:

http://www.ntuzov.com/Nik_Site/Niks_files/Research/papers/stat_arb/Ahmed_2009.pdf

and in particular page 13. I want to try and calculate lambda_t(p) = exp^(Beta^T x_t) but I am having problems with the fact beta is derived from a log-likelihood estimation and also jacobian matrices, which require newton raphson:

Now I never did a degree in maths (I did engineering) but I am reasonably mathematical. I need therefore quite a bit of help trying to work out how to actually implement the above in code. I am trying to calculate the order rate, lambda_t(p) = exp^(Beta^T x_t).

I found this, which I think may be answering my question:

http://en.wikipedia.org/wiki/Jacobi_method#Example

If the above link is no good, would someone be able to give me a really basic example of how I could calculate beta, so that I can then calculate lambda_t(p) = exp^(Beta^T x_t)?

Does anyone know what T is representing?

I wasn't sure if I need to implement newton raphson, as well as the jacobian method, or whether the latter encompasses that.

Help most appreciated,

Thanks

Answer 1

The wiki link above is okay, but when using Newton-Raphson and Jacobi with functions and vectors like yours, I would actually read this one instead. A little more complex, but that's because it's more generally applicable then the link you referenced. Your trying to find a $\beta$ that solves the MLE equation:

For the NR iteration, your iterating from $\beta$-value to $\beta$-value along the function until you reach the solution, which is the desired $\beta$-value. The Jacobian represents the derivative of the function you need to solve, in this case the MLE; it allows you to iterate properly along the curve (derivative = a function's rate of change). On the wiki page, $\beta: MLE(\beta)=0$ analogous to $x: f(x)=0$, and $f'(x)$ is analogous to $J(\beta)$.

Does anyone know what T is representing?

T is simply the transpose operation (note that it says $\beta\in\mathbb{R}^p$)- the vector $\mathbf{x_t}$ looks to me like it should be $\in\mathbb{R}^{n\times p}$, otherwise the next equation (the summation) would not have enough $\mathbf{x_i}$ column vectors.

I wasn't sure if I need to implement newton raphson, as well as the jacobian method, or whether the latter encompasses that.

Other way around: Newton-Raphson (your last equation) encompasses the Jacobian (denoted $J$; it's inverse is used).

Note The Newton Raphson technique doesn't iterate until the error goes to zero: it takes an infinite amount of time to do so. Instead, set an upper-bound (a given tolerance $\epsilon$) on the error between $\beta$-values (using the 2-norm of the difference), and use that as a while-loop condition for your $\beta_n$ iterations.

The first thing you'll need to do is come up with an "educated guess" of an initial value, call it $\beta_0$. Whenever possible, try to think of typical $\beta$ values for the given scenario, and pick something around there.

From here, you need to calculate the (1) Jacobian matrix and (2) partial derivative for that $\beta$-value, where each element is given by Eq3 and Eq2 (respectively). You can use for-loops to iterate through $jac$ and $mle$ and create the matrix/vector element-by-element (since each element requires a summation). Define functions where the only parameter you have to pass to them is $\beta$, and let them output $J(\beta_n)$ and $\frac{\partial\log L}{\partial\beta}$.

Note the following relationships:

• MLE (Eq 2) is the partial derivative of Eq 1; it's one element of a vector.
• Jacobian (Eq 3) is the partial derivative of Eq 2; it's one element of a matrix.
• The jacobian is size $p\times p$ and $mle$ is a $p$-vector, since the $\beta$-vector is $\in\mathbb{R}^p$.

From there, you have all the values you need to get $\beta_1$ (a.k.a. $\beta_{n+1}$) from Eq4. Using this new $\beta$ value, compare the change in beta to the tolerance, and iterate through the equations with the new $\beta$ value if needed.

In terms of code:

1. define $mle(\beta_n)$ and $jac(\beta_n)$ using equations 2 and 3 above
2. Initialize $\Delta\beta=10^8$

while $err\geq\epsilon$ {

$\Delta\beta = [jac(\beta_n)]^{-1}*mle(\beta_n)$;

$err = ||\Delta\beta||$;

$\beta_{n+1} = \beta_{n} + \Delta\beta$;

}

Just a quick note: although the iteration above shows the inverse of the jacobian being used, for large matrices this becomes computationally expensive. You didn't specify what language you were using, but if it's something like MATLAB, use $jac(\beta_n)\backslash mle(\beta_n)$, as gaussian elimination saves loads of time.