This code was kindly recommended to me in my original question. It returned the same parameter estimates as the software called CRAFT by Aon Benfield. I have also managed to replicate it for the Weibull distribution.
I was wondering if anyone could help me replicate if for Pareto? Based on CRAFT, I would expect a shape of ~0.59 and scale of ~55.15 (b~93 and q~1.7). Finally, I would also like to do the same but using the SSQ method instead of MLE?

My data = c(28.744, 385.714, 20.595, 99.350,31.864, 77.713, 264.408, 21.204, 31.937, 0.900, 18.762, 173.276, 23.707) 

constrained_mle <- function(x, p, q) {
  lnL <- function(shape) {
    # Solving q = qgamma(p, shape)*scale for the necssary scale
    scale <- q/qgamma(p, shape)
    sum(dgamma(x, shape, scale=scale, log=TRUE))
  }
  res <- optimise(lnL, lower=0, upper=1e+3, maximum=TRUE, tol=1e-8)
  scale <- q/qgamma(p, res$maximum)
  c(scale=scale, shape=res$maximum)
}
par <- constrained_mle(x, .95, 500.912)
par
#>       scale       shape 
#> 231.0561574   0.6038756
# checking that the solution is correct
qgamma(.95, scale=par[1], shape=par[2])
#> [1] 500.912 ```