Determine Maximum Likelihood Estimate (MLE) of loglogistic distribution

Question

I am given two data sets containing dates and losses (in some currency). I have to determine the maximum likelihood estimates of the parameters of loglogistic distribution.

I googled and found a paper that gives me the "MMLE" for sigma in (2.4).

1 Is that good enough to substitute for MLE? It seems to follow the same procedure of getting the MLE that I recall from introductory statistics.

2 To estimate sigma, it looks like I will have to censor s values from sample. How do I determine a value for s?

Exact instructions:

"Directions:

1. For each loss data set [given elsewhere],
   determine the maximum likelihood estimates of the parameters of each of the following models:
   a. Loglogistic distribution

...

2. For each model in no. 1, perform the following goodness-of-fit tests:
   a. Kolmogorov-Smirnov test b. Anderson-Darling test c. Chi-square goodness-of-fit test

3. The following must be included in your report:
   a) Word file (hard copy): Complete solution: discussion of the likelihood (or loglikelihood) function for each model, discussion of the goodness-of-fit tests for each model Computation of VaR and TVaR b) Excel file (by email): Implementation of the formulas obtained in 3.a "

Answer 1

If you are looking to fit a log-logistic distribution to your data, it is fairly straightforward to do so. In the example below, I am using the function dllog to get at the density of the log-logistic for a given set of values of the shape and scale parameter, but it is no trouble to write the PDF code yourself as well.

---

(Log-)Likelihood & MLE

The density of log-logistic distributed a random variable has the probability density function [PDF]: $$ f(X_i; \alpha, \beta) = \dfrac{\left(\tfrac{\beta}{\alpha}\right)\left(\tfrac{X_i}{\alpha}\right)^{\beta - 1}}{\left(1+ \left(\tfrac{X_i}{\alpha}\right)^{\beta}\right)^2} $$ where $\alpha$ and $\beta$ are the scale and shape parameters respectively.

For a given sample of data $X_1, \ldots, X_N$, this implies that the log-likelihood of the sample is: $$ \ell_N(\alpha, \beta \mid X_1, \ldots, X_N) = \sum_{i=1}^N \log f(X_i; \alpha, \beta) $$

The MLE of the parameters given the sample of the data, is given by the maximizer of the log-likelihood: $$ \hat{\alpha}_{MLE}, \hat{\beta}_{MLE} = \arg\max_{\alpha, \beta}\ell_N(\alpha, \beta \mid X_1, \ldots, X_N) $$

---

Computing & optimizing the log-likelihood

In the code below:
0. I have used the function rllog to generate a random sample from a log-logistic distribution with parameters c(5, 6).
1. The function fnLLLL computes the (negative) log-likelihood of the data.
2. The function fnLLLL uses the function dllog from the FAdist package to compute the PDF of the log-logistic distribution, $f$.
3. optim computes the values of $\alpha$, and $\beta$ that minimize the negative log-likelihood, and the values c(2, 3) are the intial values for the optimizer. Those optimized values are $5.132758$ & $5.654340$, and the optimized value of the negative log-likelihood function is $9239.179$.

# simulate some log-logistic data
library(FAdist)
vY = rllog(n = 1000, shape = 5, scale = 6)

# log-likelihood function
fnLLLL = function(vParams, vData) {
  # uses the density function of the log-logistic function from FAdist
  return(-sum(log(dllog(vData, shape = vParams[1], scale = vParams[2]))))
}

# optimize it
optim(c(2, 3), fnLLLL, vData = vY)

This gives:

> optim(c(2, 3), fnLLLL, vData = vY)
$par
[1] 5.132758 5.654340

$value
[1] 9239.179

$counts
function gradient 
      57       NA 

$convergence
[1] 0

$message
NULL

Answer 2

I asked above for the parametrization, as in Actuarial Science (at least in the US) the loglogistic (Fisk) is usually parameterized: $$ \begin{align} f(x) &= \frac{\gamma\left(\frac{x}{\theta}\right)^\gamma}{x + \left[1 + \left(\frac{x}{\theta}\right)^\gamma\right]^2}\\ F(x) &= \frac{\left(\frac{x}{\theta}\right)^\gamma}{1 + \left(\frac{x}{\theta}\right)^\gamma} \end{align} $$

Therefore: $$ LL(x) = \ln\gamma + \gamma\left(\ln x - \ln\theta\right) - \ln x - 2\ln\left(1 + \left(\frac{x}{\theta}\right)^\gamma\right) $$

As above, create an error function for this (well, negative log likelihood) similar to:

NLL_LGLG <- function(pars, X) {
  gamma <- pars[[1]]
  theta <- pars[[2]]
  LL <- log(gamma) + gamma * (log(X) - log(theta)) - log(X) -
        2 * log(1 + (X / theta) ^ gamma)
  return(-sum(LL))
}

You can use deriv3 if you want to use gradient based methods like L_BFGS (or calculate it by hand if speed is of the essence), but plugging the above and your data into optim or nloptr should be what you want. If you don't want to roll your own logliklihoods, and this parameterization is the one you want, you can find dllogis in the actuar package on CRAN.