The initial function is log-likelihood of multinomial probit model with $J$ alternatives:
$ ln \ell=\sum_{i=1}^N\sum_{j=1}^{J-1} y_{ij} \cdot ln \Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})+ ({n_i-\sum_{j=1}^{J-1}y_{ij}}) \cdot ln[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]$
Then, the first derivative of
$\dfrac{\partial}{\partial\beta_{kj}}\sum_{i=1}^N \left\{y_{ij} \cdot ln \Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})+ y_{iJ} \cdot ln\left[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})\right]\right\}$
would be
$\sum_{i=1}^N \left\{y_{ij} \dfrac{x_{ik}\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})}{\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})}- y_{iJ} \dfrac{x_{ik}\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})}{1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})}\right\}$
Can somene help me with the second derivative? I am aware that for each $\beta_{kj}$ I need to differentiate wrt every other $\beta_{kj}$:
$\dfrac{\partial^2}{\partial\beta_{kj}\beta_{k'j'}}=\sum_{i=1}^N \left\{y_{ij} \dfrac{x_{ik}\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})}{\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})}- y_{iJ} \dfrac{x_{ik}\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})}{1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})}\right\}$
I know that there will be two cases where $j'=j$ and $j'\ne j$ and that derivative of the first fraction in the bracket when $j'=j$ would be
$-x_{ik'}x_{ik}y_{ij}\cdot\dfrac{(\sum_{k=1}^Kx_{ik}\beta_{jk})\cdot\phi(\sum_{k=1}^Kx_{ik}\beta_{jk})}{\Phi(\sum_{k=1}^Kx_{ik}\beta_{jk})}-x_{ik'}x_{ik}y_{ij}\cdot\dfrac{\phi^2(\sum_{k=1}^Kx_{ik}\beta_{jk})}{\Phi^2(\sum_{k=1}^Kx_{ik}\beta_{jk})}$.
Any help with the second derivative of the second fraction when $j'=j$ and both fractions when $j'$ not equal to $j$ would be very much appreciated.
* EDIT: 2 * this is my attempt of doing second derivative of the probit log-likelihood function:
The second derivative of the log-likelihood function should have the following form:
$\dfrac{\partial^2}{\partial\beta_{kj}\beta_{k'j'}}=\sum_{i=1}^N\left\{y_ij\dfrac{\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot\varphi'(\sum_{k=1}^Kx_{ik}\beta_{kj})-\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot\Phi'(\sum_{k=1}^Kx_{ik}\beta_{kj})}{\Phi^2(\sum_{k=1}^Kx_{ik}\beta_{kj})}+y_{iJ}\dfrac{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]\cdot\varphi'(\sum_{k=1}^Kx_{ik}\beta_{kj})-\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot[1-\sum_{j=1}^{J-1}\Phi'(\sum_{k=1}^Kx_{ik}\beta_{kj})]'}{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]^2} \right\}$
Hence, the derivatives of the numerator and denominator will depend whether $j'=j$ (Case 1) or $j'\ne j$ (Case 2).
Case 1 $j'=j$
$\varphi'(\sum_{k=1}^Kx_{ik}\beta_{kj})=-x_{ik'}\cdot\sum_{k=1}^Kx_{ik}\beta_{kj}\cdot\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})$
$\Phi'(\sum_{k=1}^Kx_{ik}\beta_{kj})=x_{ik'}\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})$
$[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]'=-x_{ik'}\cdot\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})$
2) case $j'\ne j$
$\varphi'(\sum_{k=1}^Kx_{ik}\beta_{kj})=0$
$\Phi'(\sum_{k=1}^Kx_{ik}\beta_{kj})=0$
$[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]'=-x_{ik'}\cdot\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj'})$
By substituting these terms in the main equation for the case (1) yields
$=-\sum_{i=1}^N\left\{y_{ij}\dfrac{(\sum_{k=1}^Kx_{ik}\beta_{jk})\cdot\phi(\sum_{k=1}^Kx_{ik}\beta_{jk})}{\Phi(\sum_{k=1}^Kx_{ik}\beta_{jk})}+y_{ij}\dfrac{\phi^2(\sum_{k=1}^Kx_{ik}\beta_{jk})}{\Phi^2(\sum_{k=1}^Kx_{ik}\beta_{jk})}+y_{iJ}\dfrac{(\sum_{k=1}^Kx_{ik}\beta_{jk})\cdot\phi(\sum_{k=1}^Kx_{ik}\beta_{jk})}{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]}+y_{iJ}\dfrac{\phi^2(\sum_{k=1}^Kx_{ik}\beta_{jk})}{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]^2}\right\}\cdot x_{ik}x_{ik'}$
By solving for case (2), first fraction is equal to $0$
$\sum_{i=1}^N\left\{y_{ij}\dfrac{\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot0-\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot0}{\Phi^2(\sum_{k=1}^Kx_{ik}\beta_{kj})}+y_{iJ}\dfrac{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]\cdot0-\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot(-x_{ik'})\cdot\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj'})}{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]^2} \right\}$
And the result is
$=-\sum_{i=1}^N\left\{y_{iJ}\dfrac{\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj})\cdot x_{ik'}\cdot\varphi(\sum_{k=1}^Kx_{ik}\beta_{kj'})}{[1-\sum_{j=1}^{J-1}\Phi(\sum_{k=1}^Kx_{ik}\beta_{kj})]^2} \right\}$
Is this correct?
With the second fraction (although the working is the same as with the first fraction),
$$y_{ij}\cdot \dfrac{x_{ik}\cdot\varphi(\sum\limits_{k=1}^Kx_{ik}\beta_{kj})}{1-\sum\limits_{j=1}^{J-1}\Phi(\sum\limits_{k=1}^Kx_{ik}\beta_{kj})}$$
just continue with the quotient rule as before. The denominator is just a sum of functions, and so presents no real problem.
As for the $j,j'$ issue, it makes no difference whether $j=j'$ or $j\ne j'$ - both cases can be differentiated in the same way. However, due to the nature of partial differentiation, if $j\ne j'$, then $\dfrac{\partial}{\partial b_{j}}b_{j'}=0$ as $j$ and $j'$ can be considered to be independent from each other.
And the 'J' in $y_{iJ}$ should be in lowercase.