T-distribution parameters with QRM package I am fitting a t-distribution on some data I have using the fit.st function from the QRM package.
The function returns 2 set of parameters (par.ses and par.ests), both containing values for nu, mu, and sigma. The values for the 2 sets of parameters are quite different and I don't get which one to use.
I looked into the source code of the package function and from my understanding they solve the optimization problem through 2 different methods and return the 2 results regardless of the success of the first one. In that case I should take the par.ests results if converged==TRUE and the par.ses results if the convergence didn't happen.
Is my understanding correct?
Here is the code of the fit.st function in the package:
fit.st <- function(data, ...){
  if(is.timeSeries(data)) data <- series(data)
  mu <- mean(data)
  m2 <- mean((data - mu)^2)
  m4 <- mean((data - mu)^4)
  nu <- 4 + (6 * m2^2) / (m4 - 3 * m2^2)
  sigma <- sqrt((nu - 2) * m2 / nu)
  theta <- c(nu, mu, sigma)
  negloglik <- function(theta, y){
    - sum(log(dt((y - theta[2]) / abs(theta[3]), df = abs(theta[1]))) - log(abs(theta[3])))
  }
  optimfit <- optim(theta, fn = negloglik, y = data, ...)
  par.ests <- optimfit$par
  ifelse(optimfit$convergence == 0, converged <- TRUE, converged <- FALSE)
  par.ests[1] <- abs(par.ests[1])
  par.ests[2] <- abs(par.ests[2])
  nItheta <- hessian(negloglik, par.ests, y = data)
  asymp.cov <- solve(nItheta)
  loglh.max <- -negloglik(par.ests, y = data)
  par.ses <- sqrt(diag(asymp.cov))
  names(par.ests) <- c("nu", "mu", "sigma")
  names(par.ses) <- names(par.ests)
  dimnames(asymp.cov) <- list(names(par.ests), names(par.ests))
  list(converged = converged, par.ests = par.ests, par.ses = par.ses,
       asymp.cov = asymp.cov, ll.max = loglh.max)
}

and the link to the package documentation: https://rdrr.io/cran/QRM/src/R/Student.R

par.ests referes to the parameter estimates. par.ses refers to the parameter standard errors.
It is maybe not so obvious, since the function has no documentation. But you can see in these lines:
  nItheta <- hessian(negloglik, par.ests, y = data)
  asymp.cov <- solve(nItheta)
  loglh.max <- -negloglik(par.ests, y = data)
  par.ses <- sqrt(diag(asymp.cov))

that the function uses the hessian matrix to compute standard errors.
