# Maximum likelihood of multivariate t-distributed variable with scaled covariance

I am trying to estimate the covariance of a iid multivariate t-distributed random variable, where I define the multivariate density as in the Statlect textbok, which is the same as the wikipedia page. The density in question is given by, $$f_{X}\left(x\right)=\frac{\Gamma\left(\frac{v+p}{2}\right)}{\Gamma\left(\frac{v}{2}\right)\sqrt{v^{p}\pi^{p}\left|\det \Sigma\right|}}\left(1+\frac{1}{v}x^{\prime}\Sigma^{-1}x\right)^{-\frac{v+p}{2}}$$ where $x$ is a realization of the $p\times1$ random vector $X$ with covariance $V[X]=\frac{v}{v-2}\Sigma$, and $v$ is the degrees of freedom parameter.

I'm interested in estimating the covariance of $Y=\sqrt{\frac{v-2}{v}}X$ which, for $v>2$, has covariance $$V[Y]=V\left[\sqrt{\frac{v-2}{v}}X\right]=\frac{v-2}{v}V[X]=\Sigma$$

From density tranformation it holds that the density of $Y$ is, $$f_Y(y)=\frac{1}{|\sqrt{\frac{v-2}{v}}|}f_X\left( y /\sqrt{\frac{v-2}{v}} \right)$$ which i have calculated to, $$f_{Y}\left(y\right)=\frac{\Gamma\left(\frac{v+p}{2}\right)}{\Gamma\left(\frac{v}{2}\right)\sqrt{(v-2)v^{p-1}\pi^{p}\left|\det\Sigma \right|}}\left(1+\frac{1}{v-2}y^{\prime}\Sigma^{-1}y\right)^{-\frac{v+p}{2}}$$ so my best guess at the log-likelihood is: $$L(\theta)=-\frac{1}{2}\times\sum_{i=1}^N\left[-2\delta(v)+p\log(\pi)+\log(v-2)+(p-1)\log(v)+\log(\det(\Sigma))+(v+p)\log\left(1+\frac{1}{v-2}y_i^{\prime}\Sigma^{-1}y_i \right) \right]$$ where $\delta(v)=\log\frac{\Gamma((v+p)/2)}{\Gamma(v/2)}$.

Are the derivations of the log-likelihood correct? Is this the likelihood of a variable $Y$ with covarians $\Sigma$, and can i expect the mariginals to be standard univariate t-distributions with variances $\text{diag}(\Sigma)$?

I have implemented this likelihood in R and tried to estimate on simulated data, but when I for example simulate with v=5, I get an estimate of around 50 to 60, like in the example below.

I'm quite confident about both the log-likelihood, as well as the 'R'-code, but the results don't add up, so clearly one or the other (or both) is worng! Can you help me figure out what I misunderstand.

Thank you in advance!

The simulation in R as an illustration:

N <-10000
v <- 5
set.seed(123)
sigma <- matrix(c(1,0.5,
0.5,1),2,2,byrow=T)

U <- chol(sigma)
Z <- matrix(rt(N*2,v),N,2)
Y <- Z%*%U*sqrt((v-2)/v)

# Check covariance
var(Y)
#          [,1]      [,2]
#[1,] 0.9764366 0.4701139
#[2,] 0.4701139 0.9515428

# Defining the (negative) log-likelihood function
lik.t <- function(theta,data){
# Number of obs
N <- dim(data)[1]
# Nummber of dimensions (set to 2)
p <- 2

# Empty vector for likelihood values
lik <- rep(NA,N)

# Lower Choleski parametrization
L <- matrix(c(theta[1:2],0,theta[3]),2,2)
sigma <- L%*%t(L)

# Ensuring v>2
v <- exp(theta[4])+2

delta <- log(gamma((v+p)/2))-log(gamma(v/2))
for(i in 1:N){
lik[i] <- 1/2*(-2*delta + p*log(pi) + log(v-2) + (p-1)*log(v) + log(det(sigma)) +
(v+p)*log(1 + 1/(v-2)*t(data[i,])%*%solve(sigma)%*%data[i,]))
}
sum(lik)
}

# Colecting tru parameters as starting value for BFGS
foo <- chol(sigma)[upper.tri(sigma,T)]
theta0 <- c(foo,v)

# Miniminzing minus the log-likelihood
out <- optim(theta0,lik.t,data=Y,method="BFGS",control=list(trace=1,REPORT=1))

#Printing results
L.hat <- matrix(0,2,2)
L.hat[lower.tri(L.hat,T)] <- out$par[1:3] sigma.hat <- L.hat%*%t(L.hat) v.hat <- exp(out$par[4])+2

sigma.hat  # Estimated covaraince
#         [,1]      [,2]
#[1,] 0.9347044 0.4534481
#[2,] 0.4534481 0.9088294
v.hat      # Estimated degrees of freedom
#[1] 51.53331

• Are you asking about the statistical aspect of this, or are you looking for code check? Note that the latter is off-topic here, but we can migrate it for you. – gung - Reinstate Monica Nov 25 '14 at 3:58
• I'm not looking for a code check. I'm would like to ask if the likelihood function is correctly specified for the given problem and if what I'm writing is true, in a statistical sence.But maybe i need to clarify :-) – Duffau Nov 25 '14 at 4:10

## 1 Answer

The EM algorithm is typically used to find MLEs of the parameters from an iid multivariate t sample. Instead of writting a long answer as to how to implement the algorithm I refer you to McLachlan and Krishnan "The EM algorithm and Extensions", second edition. They show the MLE procedure for both known and unknown degrees of freedom, as well as ways to accelerate the algorithm. If you don't have access to the book you can also look at this paper

https://arxiv.org/pdf/1707.01130.pdf

It turns out that estimating the degrees of freedom using ML results in an unbounded score function. The above paper finds estimators that minimize the MLq, which results in a bounded score. They also have a great review of the EM algorithm for the usual ML procedure.

Once you have estimators for both the degrees of freedom and the scale parameter, finding an estimate of the covariance is straight forward.

• When I asked this question originally I was working on a vector-GARCH project, where I was seeking to replace a multivariate normal assumption with a multivariate t- distribution. I remember that in the univariate GARCH literature this is done routinely, where naive ML techniques are used to both fit parameters of the GARCH process and the degrees of freedom simultaneously. Do you have a clue to why that succeeds? Is this also an issue with univariate iid t-variate? Thanks for the input! I never really let this go :-) – Duffau Jan 22 '19 at 18:44
• Hey @Duffau, I am not sure about the univariate case, but I suggest you look into the book I told you. Formulating a problem using multivariate t errors is straight forward once you have formulated it using multivariate normal errors. Whatever estimation procedure you have for the multivariate normal case is very similar to the M step of the EM algorithm. Then the E step just estimated the latent variables, which are usually gamma(nu/2,nu/2), where nu are the degrees of freedom. Just be careful with the intercept and the sums. – Carlos Llosa Jan 30 '19 at 18:42