Method of moments and MLE estimates for Lomax (Pareto Type 2)

Question

I have this dataset, on which I am supposed to fit Lomax distribution with MM and MLE. Lomax pdf is: $$f(x|\alpha, \lambda) = \frac{\alpha\lambda^\alpha}{\left(\lambda+x\right)^{\alpha+1}}$$ For MM, it is possible to show that: $$\hat{\alpha}=\frac{2\hat{\sigma}^2}{\hat{\sigma}^2-\bar{X}^2}$$ $$\hat{\lambda}= \bar{X}\frac{\hat{\sigma}^2+\bar{X}^2}{\hat{\sigma}^2-\bar{X}^2}$$ Where $\hat{\sigma}^2$ is the sample variance and $\bar{X}$ is sample mean. The estimates are:

df <- read.table('Theft.txt')
df <- df$V1

s <- var(df)
m <- mean(df)

alpha <- (2*s)/(s-m^2)
lambda <- m*((s+m^2)/(s-m^2))

> alpha
[1] 2.70862
> lambda
[1] 3451.911

For MLE, I have log-likelihood function: $$\ell(\alpha, \lambda|x) = n\log(\alpha)+\alpha n\log(\lambda) - (\alpha+1)\sum_{i=1}^{n}\log(\lambda+x_i)$$

and the implementation:

llik <- function(alpha, lambda,x){
  n<-length(x)
  res <- n*log(alpha)+n*alpha*log(lambda)-(alpha+1)*sum(log(x+lambda))
  return(-res)
}
mle1 <- mle(minuslogl = llik, start = list(alpha=alpha,lambda=lambda),
fixed = list(x=df), method = 'BFGS')
> mle1@coef
      alpha      lambda 
   2.860708 3451.907162

I used as starting values the MM estimates. The resulting coefficients are quite similar to MM, however after using flomax() function from Renext package, I am getting completely different estimates, with higher likelihood:

> flomax(df)
$estimate
      shape       scale 
   1.880468 1872.132104 

I have also done some simulations, in which both MM and MLE are really bad at estimating the 'real' parameters of Lomax. Why are these estimates this bad? Why is in my case MM so different from MLE? Why is mle() so sensitive to starting values?

Answer 1

The issue appears to be the greatly different scales of the two parameters and how that interacts with BFGS. When I try optim using BFGS on the raw data, I get similar results to mle above (not surprisingly):

x <- df / 1000

llik <- function(theta, x){
   alpha <- theta[1]
   lambda <- theta[2]
   n<-length(x)
   res <- n*log(alpha)+n*alpha*log(lambda)-(alpha+1)*sum(log(x+lambda))
   return(-res)
}

alpha <- 2.7
lambda <- 3450    
mle1 <- optim(c(alpha, lambda), llik, method="BFGS", x = 1000*x)

mle1$par
[1]    2.859574 3449.996428

But working with the rescaled data:

alpha <- 2.7
lambda <- 3.450
mle1 <- optim(c(alpha, lambda), llik, method="BFGS", x = x)

mle1$par
[1] 1.880470 1.872135

llik(c(mle1$par[1], 1000*mle1$par[2]), 1000*x)
[1] 1012.211

Using a different technique (Nelder-Mead) on the original data gives good results, although we really ought to rewrite the log likelihood function so as to not fail when negative values of the two parameters are passed:

alpha <- 2.7
lambda <- 3450

mle1 <- optim(c(alpha, lambda), llik, method="Nelder-Mead", x = 1000*x)
Warning messages:
1: In log(alpha) : NaNs produced
2: In log(alpha) : NaNs produced
3: In log(alpha) : NaNs produced
4: In log(alpha) : NaNs produced
5: In log(alpha) : NaNs produced
6: In log(alpha) : NaNs produced
7: In log(alpha) : NaNs produced

mle1$par
[1]    1.879401 1870.984994