The multivariate skewed variance gamma (VG) distribution is useful for modelling data with heavy-tails and high density around the location parameter. When the shape parameter is sufficiently small, the density function is unbounded at the location parameter. In this thesis, we proposed three modifications to appropriately bound the likelihood function so that the maximum is well-defined. These modified likelihoods are the capped, leave-one-out (LOO), and weighted LOO likelihoods. Moreover, we present expectation/conditional maximisation (ECM) algorithms to accurately estimate parameters of the VG distribution using its normal mean-variance mixture representation. Apart from parameter estimation, we also calculate standard errors (SEs) to assess the significance of the parameter estimates. However, this calculation requires the second order derivative of the log-likelihood function with respect to vector/matrices. We apply new matrix differentiation formulas to efficiently compute the observed and Fisher information matrices for the VG distribution. These SE calculations rely on asymptotic properties of the maximum likelihood estimator (MLE) which have been extensively studied under the smooth likelihood case. For the cusp/unbounded case, proving these asymptotic properties are a challenge as they do not satisfy the smoothness regularity condition. We numerically investigate these asymptotic properties for the location estimator when the likelihood function has cusp or unbounded points. We demonstrated its super-efficient rate of convergence and found the double generalised gamma distribution provides a good approximation to the asymptotic distribution of the location parameter. Lastly, the ECM algorithms are applied to the vector autoregressive moving average model with VG and Student's t innovations to capture serial correlation, leptokurtosis, skewness, and cross dependence of return data from high frequency stock indices and cryptocurrencies.