Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
136,985 characters · 24 sections · 62 citation commands
Convolution-$t$ Distributions
Keywords:{ Multivariate heavy-tailed distributions. Convolutions of $t$-distributions, Voigt profile.}
JEL Classification:{ C01, C32, C46, C58}
Heavy-tailed distributions play a central role in the modeling of risk and extreme events in economics, finance, and beyond, and their multivariate properties are key to risk management, portfolio optimization, and the assessment of systemic risk, see e.g. EmbrechtsKluppelbergMikosch:1997, Harvey2013, and IbragimovIbragimovWalden:2015.
This paper introduces a new class of multivariate heavy-tailed distributions, convolution-$t$ distributions, that are convolutions of mutually independent multivariate $t$-distributions. The new distributions can accommodate heterogeneous marginal distributions and various forms of nonlinear dependencies. Moreover, the nonlinear dependencies can be used to identify factors of economic significance. Conveniently, the expression for the density function of convolution-$t$ distributions is simple. This facilitates straight forward estimation and inference. The framework includes the multivariate-$t$ distribution as a special case, and makes clear how it is restrictive in a number of ways. For instance, the multivariate-$t$ distribution implies all marginal distributions have identical shape (e.g. same kurtosis) and it implies a particular type of nonlinear dependence between elements, because they all share a single common Gamma-distributed (mixing) variable.
A convolution-$t$ distribution, denoted $\mathrm{CT}_{\boldsymbol{n},\boldsymbol{\nu}}(\mu,\Xi)$, is characterized by a location vector, $\mu$, a scale-rotation matrix, $\Xi$, and two $K$-dimensional vectors, $\bm{n}$ and $\bm{\nu}$, that specify cluster sizes and degrees of freedom, respectively. The matrix, $\Xi$, define key features of the distribution and the resulting non-linear dependencies identify the latent factor structure. We provide a novel identifying representation of $\Xi$, which can accommodate block structure and other sparse structures when such are needed. This makes the framework amenable to high-dimensional applications of the convolution-$t$ distributions.
The main contributions of this paper are as follows: First, we introduce a family of convolution-$t$ distributions and obtain their marginal density and cumulative distribution function by means of the characteristic function. Second, we characterize the identification problem and show that convolution-$t$ distributions can help unearth latent factor structures, which is not possible from the covariance structure alone. Third, we derive the score, hessian, and information matrix from the log-likelihood function. And the proof of consistency and asymptotic normality of maximum-likelihood estimators are well established. The asymptotic properties are confirmed in a simulation study. Fourth, we extend the moment-based approximation method by Patil:1965 to convolutions of any number of $t$-distributions. Fifth, we show that the convolution-$t$ distributions provide substantially empirical gain in an application with ten financial volatility series. Importantly, the convolution-$t$ framework adds insight about the latent factor structure, where nonlinear dependencies identify a common market factor, and a cluster structure that aligns with sector classifications.
Convolutions of $t$-distributions arise in classical problems, such as the Behrens-Fisher problem (statistics) and the Voigt profile (spectroscopy), where the latter is the convolution of a Gaussian distribution and a Cauchy distribution.\footnote{This convolution is closely related to a Mills ratio that appears in the Heckman model.} The literature has mainly focused on convolutions of two $t$-distributions, see e.g. Chapman:1950, Ghosh:1975, and PrudnikovBrychkovMarichev:1986.\footnote{Ruben:1960 expressed the density function as an integral that, in the general case, involves hypergeometric functions. RahmanSaleh:1974 and Dayal:1976 derived an expression for the distribution of the Behrens-Fisher statistic using Appell series.} BergVignat:2010 derived an expression for the density, which includes convolutions of two multivariate $t$-distributions. Their expression comprises an infinite sum with coefficients that are given from integrals.
The convolution-$t$ distributions offer a parametric approach to modeling multivariate variables with complex dependencies and heterogeneous marginal distributions, similar to copula-based methods, see e.g. Patton:Copula2012 and FanPatton:Copula2014. For instance, Fang2002 proposed the meta-elliptical distributions, which construct the joint density function by combining elliptical copula function with certain marginals, e.g. Student's $t$-distribution. The Meta-$t$ distribution of Demarta2007 is often used in high dimensional settings, see e.g. OhPatton2017JBES\nocite{OhPatton2017}\nocite{OhPatton2023}, OpschoorLucasBarraVanDick:2021, and CrealTsay2015. Another related strand of literature is that on multivariate stochastic scale mixture of Gaussian family distributions, see e.g. Eltoft2006, FinegoldDrton:2011, and Forbes2014. However, for the aforementioned two types of multivariate distributions, only the scaling (or covariance) matrix matters. The structure of convolution-$t$ distributions is differ from existing distributions in many ways, with an important one being the underlying cluster structure in convolution-$t$ distributions that, in conjunction with a novel scale-rotation matrix, $\Xi$, define particular nonlinear dependencies that can be used to identify clusters and latent factor structures in the data.
The rest of this paper is organized as follows. We introduce the convolution-$t$ distributions in Section (ref) and establish a number of its properties, including moments and marginal densities. We derive the score, Hessian, and Fisher information from the log-likelihood function in Section (ref), and characterize an interesting identification problem in this model. We consider the finite-sample properties of the maximum likelihood estimator (MLE) in Section (ref), and confirm that the analytical expressions for standard errors, based on the asymptotic expressions, are reliable. We introduce a standardized variant of convolution-$t$ distribution with finite variance in Section (ref), which makes it easier to interpret estimators in the empirical analysis. We derive a moment-based approximation methods of marginal convolution-$t$ distributions in Section (ref). In Section (ref), we apply the convolution-$t$ framework to model a vector of realized volatilities. We conclude in Section 8 and present all proofs in the Appendix. Some additional theoretical results and simulation results for non-standard situations are presented in the `Supplemental Material'.
The $n$-dimensional Student's $t$-distribution has density function,
where $\nu>0$ is the degrees of freedom, $\mu\in\mathbb{R}^{n}$ is the location parameter, and $\Sigma\in\mathbb{R}^{n\times n}$ is the symmetric and positive definite scale matrix. We write $X\sim t_{n,\nu}(\mu,\Sigma)$ to denote a random variable with this distribution and include the limited case as $\nu\rightarrow\infty$, such that $t_{\infty,n}(\mu,\Sigma)$ represents the Gaussian distribution, $X\sim N_{n}(\mu,\Sigma)$, with density,
If $X\sim t_{n,\nu}(\mu,\Sigma)$, then it is easy to verify that
for any $m\times n$ matrix $B$ with full row rank. Any marginal distributions of a multivariate $t$-distribution is a univariate $t$-distribution with the same degrees of freedom, $\nu$, as that of the multivariate distribution. This can be too restrictive and rules out marginal distributions with heterogeneity in the kurtosis. Recall that a multivariate $t$-distribution has the representation: $\sqrt{\frac{\nu}{\xi}}Z\sim t_{n,\nu}(\mu,\Sigma)$, where $Z\sim N_{n}(\mu,\Sigma)$ and $\xi\sim\mathrm{Gamma}(\tfrac{\nu}{2},2)$ are independent. This highlights another characteristic of the multivariate $t$-distribution, which is that common mixing variable $\xi$, induces very particular nonlinear dependencies. The convolution-$t$ distribution, which we introduce next, is a more versatile class of distributions.
Consider now the case where $X$ is composed of $K$ independent multivariate $t$-distributions, \[ X=\left(
\right),\qquadwhere\quad X_{k}\sim t_{n_{k},\nu_{k}}(0,I_{k}),\qquadfor\quad k=1,\ldots,K, \] such that the dimension of $X$ is $n=\sum_{k=1}^{K}n_{k}$. Because $X_{1},\ldots,X_{K}$ are mutually independent, it follows that the density function for $X$ is given by
where $f_{X_{k}}(x_{k})$ has the form in ((ref)) if $\nu_{k}<\infty$ and the form ((ref)) if $\nu_{k}=\infty$. We will use the following notation for convolutions of heterogeneous multivariate $t$-distributions (including Gaussian distributions).
An important characteristic of the convolution-$t$ distribution is that the multivariate density of $Y=\mu+\Xi X$ is very simple when $\Xi$ is an invertible matrix. In this case, we have
The simple expression makes the analysis of the log-likelihood function straight forward.
The properties of any convolution-$t$ distribution can be deduced from the case where $\Xi$ is invertible (which implies $m=n$). The case $m>\mathrm{rank}(\Xi)$ implies perfect collinearity in $Y$, such that the properties can be inferred from a lower-dimensional subvector. For the case $m<n$, we can introduce $m-n$ auxiliary $Y$-variables, which can be integrated out to obtain the $m$-dimensional distribution. We derived detailed results for the special case, $m=1$, which represents a marginal distribution of any convolution-$t$ distribution.
If $X\sim t_{n,\nu}(0,I_{n})$, then it is well known that $\Xi X$ and $\tilde{\Xi}X$ are observationally equivalent whenever $\Xi\Xi^{\prime}=\tilde{\Xi}\tilde{\Xi}^{\prime}$. This situation is different for convolution-$t$ distributions, where more structure in $\Xi$ can be identified. We illustrate this with a simple trivariate example, which also highlights some properties of convolution-$t$ distributions.
Let $n_{1}=1$ and $n_{2}=2$ such that $Y=\mu+\Xi X\in\mathbb{R}^{3}$, where $X_{1}\sim t_{\nu_{1}}(0,1)$ and $X_{2}\sim t_{\nu_{2}}(0,I_{2})$, and suppose that $\nu_{1},\nu_{2}>2$ such that the variances, $\mathrm{var}(\sqrt{\tfrac{\nu_{1}-2}{\nu_{1}}}X_{1})=1$ and $\mathrm{var}(\sqrt{\tfrac{\nu_{2}-2}{\nu_{2}}}X_{2})=I_{2}$ are well-defined. Thus with \[ \Xi=A\left[
\right], \] we have $\mathrm{var}(Y)=AA^{\prime}$. It is easy to verify that \[ A_{sym}=\left[
\right],\quadand\quad A_{asym}=\left[
\right], \] both result in the same covariance matrix,\footnote{The $A$-matrices are presented with approximate numerical values for readability. The exact values are $1/(3\sqrt{2})\approx0.236$, $1/\sqrt{2}\approx0.707$, and $2\sqrt{2}/3\approx0.943$.} which is given by \[ \mathrm{var}(Y)=\left[
\right],\quad{\rm with}\quad\rho=\tfrac{1}{2}. \] The symmetric $A$-matrix, $A_{\mathrm{sym}}$, is the symmetric square root of $\mathrm{var}(Y)$ and $A_{\mathrm{asym}}$ is an asymmetric matrix. These two $A$-matrices would result in the exact same distribution of $Y$, if $X$ had a multivariate $t$-distribution (including the Gaussian distribution). For convoluted $t$-distributions the two $A$-matrices lead to very different distributions for $Y$. In this example, we intentionally chose $A_{\mathrm{asym}}$ to have a symmetric $2\times2$ lower-right submatrix, because this emerges as an identifying assumption in our likelihood analysis.
The joint density function of $Y$ is simply given by \[ f_{Y}\left(y\right)=\frac{c_{1}c_{2}}{|\det\Xi|}\left(1+\tfrac{1}{\nu_{1}}x_{1}^{2}\right)^{-\frac{\nu_{1}+1}{2}}\left(1+\tfrac{1}{\nu_{2}}x_{2}^{\prime}x_{2}\right)^{-\frac{\nu_{2}+2}{2}},\qquad\text{for }y\in\mathbb{R}^{3}, \] where $\left(x_{1},x_{2}\right)^{\prime}=\Xi^{-1}y$ and $c_{i}=c(\nu_{i},n_{i})$ with $c(\nu,n)=(\nu\pi)^{-\frac{n}{2}}\Gamma\left(\tfrac{\nu+n}{2}\right)/\Gamma\left(\tfrac{\nu}{2}\right)$. While the covariance matrices for $(Y_{1},Y_{2})$, $(Y_{1},Y_{3})$, and $(Y_{2},Y_{3})$, are identical, this is not the case for the bivariate densities. We do have $f_{Y_{1},Y_{3}}=f_{Y_{1},Y_{2}}$, but \[ f_{Y_{1},Y_{2}}\left(y_{1},y_{2}\right)=\int_{-\infty}^{+\infty}f_{Y}\left(y\right){\rm d}y_{3},\qquad f_{Y_{2},Y_{3}}\left(y_{2},y_{3}\right)=\int_{-\infty}^{+\infty}f_{Y}\left(y\right){\rm d}y_{1}, \] are different and the bivariate distributions also depend on the two $A$-matrices.
Contour plots of the bivariate distributions of convolution-$t$ distributions are presented in Figure (ref), panels (c)-(f), for the case where $\nu_{1}=4$ and $\nu_{2}=8$. Panels (c) and (d) display the distributions of $\left(Y_{1},Y_{2}\right)$ and $\left(Y_{2},Y_{3}\right)$ for the symmetric $A$-matrix and panels (e) and (f) are for the asymmetric $A$-matrix. For comparison, we include the corresponding bivariate distributions when $X\sim N(0,I_{3})$ in panel (a) and for $X\sim t_{6}(0,I_{3})$, in panel (b).\footnote{We only display one bivariate distribution for each of the cases where $X\sim N(0,I_{3})$ and $X\sim t_{6}(0,I_{3})$, because they are identical for all pairs of variables and both $A$-matrices.}
The contour plots reveal many interesting features of convolution-$t$ distributions. First, the convolution-$t$ distributions can generate non-elliptical distributions, which is not possible with a multivariate $t$-distribution. The Gaussian and multivariate $t$ lead to particular forms of tail dependence, and the correlation, $\rho=1/2$, implies that their probability is more concentrated along the 45$^{\circ}$ line. Second, the convolution-$t$ distribution can produce very heterogeneous bivariate distributions, even if the $A$-matrix is symmetric, as can be seen from panels (c) and (d). The distribution in panel (d) has quasi-elliptical shape, however this is not always the case for a pair of variables in the same group, even if $A$ is symmetric. Third, the asymmetric $A$-matrix lead to contour plots with additional type shapes and different degrees of tail dependence, such as the very high level of tail dependence in panel (f). Many other bivariate distributions can be generated by varying the choice of asymmetric $A$-matrix, while preserving the covariance matrix.
Next, we characterize the marginal distribution of a convolution-$t$ random variable.
Consider the univariate convolution of $t$-distributions, $Y_{1}=\mu+\beta^{\prime}X\in\ensuremath{\mathbb{R}}$, where $\beta\in\mathbb{R}^{n}$. We seek to characterize the distribution of $Y_{1}$, that could represent an element of $Y$. We achieve this by means of the characteristic function for $Y_{1}$.\footnote{This approach as has a long history, see e.g. Gurland1948, Gil1951, Imhof1961, Bohman1970, Davies1973, Shephard:1991\nocite{Shephard1991Numerical}, and Waller1995. This approach is also commonly used in for derivative pricing when density function is unavailable in closed-form, see e.g. Heston1993, HestonNandi2000, and BakshiMadan2000.}
The characteristics function for $t_{\nu}(0,1)$ is given by
where \[ K_{\nu}(x)=\tfrac{1}{2}\int_{0}^{\infty}u^{\nu-1}e^{-x\frac{u+u^{-1}}{2}}\mathrm{d}u, \] is the modified Bessel function of the second kind.\footnote{The modified Bessel function can be expressed in many ways and is, confusingly, sometimes called the “modified Bessel function of the third kind”, e.g. Hurst:1995. It is also known as the Basset function, the Macdonald function, and the modified Hankel function.} The expression for $\phi_{\nu}(s)$, ((ref)), is due to Hurst:1995 and Joarder:1995, see Gaunt:2021 for a very elegant proof. We also include the standard normal, $N(0,1)$, in the analysis with the convention $\phi_{\infty}(s)\equiv e^{-s^{2}/2}$. For $\nu=1$ (the Cauchy distribution) the expression in ((ref)) simplifies to $\phi_{1}(s)=e^{-|s|}$. Convolution-$t$ distribution with odd degrees of freedom typically have much simpler expressions, because their characteristic function is simpler, see the example in Section (ref).
Now we turn to the interesting case where $X$ is composite and derive the characteristic function of $Y_{1}=\mu+\beta^{\prime}X=\mu+\sum_{k=1}^{K}\beta_{k}^{\prime}X_{k}$, which we use to obtain expressions for its density and cumulative distribution function.
By the Gil-Pelaez inversion theorem we now have the following expressions for the marginal density and cumulative distribution function for $Y_{1}$. These are given by
respective, where ${\rm Re}\left[z\right]$ and ${\rm Im}\left[z\right]$ denotes the real and imaginary part of $z\in\mathbb{C}$, respectively.
The main advantage of these expressions is that there is just a single variable, $s$, to be integrated out, regardless of the dimensions, $n$, and the underlying number of independent $t$-distributions, $K$. In this case, $Y_{1}$ is a linear combination of $n=\sum_{k}n_{k}$ variables, and the conventional approach to obtain its marginal distribution is to integrate our $n-1$ variables. This is computationally impractical unless $n$ is small. Thus, the expressions in ((ref)) are likely to be computationally advantageous for $n\geq3$, whereas the conventional approach is simpler when $n=2$, assuming that the densities of the underlying variables are readily available, as is the case for $t$-distributed random variables.
A key feature of the multivariate convolution-t distribution is that it can generate heterogeneous marginal distributions with different levels of heavy tails. The kurtosis of $Y_{1}$ is well-defined if $\nu_{\min}=\min_{k}\nu_{k}>4$, and, as shown below, the excess kurtosis is a linear combination of the excess kurtosis of $X_{1},\ldots,X_{K}$, where the weights are determined by vector $\beta$ and $\boldsymbol{\nu}$.
This shows that the convolution-$t$ distribution can have heterogeneous marginal distributions with different shapes. In Section (ref), we use the results of Theorem (ref) to approximate the distributions of $Y_{1}$, using a univariate $t$-distribution
where the first four moments (assumed to be finite) are matched by setting \[ \mu_{\star}=\mu,\quad\nu_{\star}=4+\tfrac{6}{\kappa_{Y_{1}}},\quad\sigma_{\star}^{2}=\omega^{\prime}\omega\tfrac{\nu_{\star}-2}{\nu_{\star}}. \]
The approximating distribution can be used to simplify calculating of several interesting quantities, such as the Value-at-Risk (VaR) and the Expected Shortfall (ES), given by \[ {\rm VaR}_{\alpha}^{\star}\left(Y_{1}\right)=F_{\star}^{-1}\left(\alpha\right),\quad\mathrm{ES}_{\alpha}^{\star}\left(Y_{1}\right)=\mu_{\star}-\sigma_{\star}\frac{\nu_{\star}+\left(F_{\star}^{-1}\left(\alpha\right)\right)^{2}}{\alpha\left(\nu_{\star}-1\right)}f_{Y_{1}}^{\star}\left(F_{\star}^{-1}\left(\alpha\right)\right), \] where $F_{\star}^{-1}\left(\alpha\right)$ is the $\alpha$-quantile of the univariate $t$-distribution $f^{\star}(y)$.
Because a convolution of symmetric distributions is symmetric, it follows that all odd moments (less than $\nu_{\min}$) are zero. General moments, including fractional moments, $r>0$, can be obtained with the method in Kawata1972. For instance, for $0<r<2$, we have \[ \mathbb{E}|Y_{1}|^{r}=\delta(r)\int_{0}^{\infty}s^{-(1+r)}\left(1-{\rm Re}[\varphi_{Y_{1}}(s)]\right)\mathrm{d}s,\qquad\delta(r)=\frac{r(1-r)}{\Gamma(2-r)}\frac{1}{\sin(\tfrac{\pi}{2}(1-r))}, \] with the convention $\delta(1)=2/\pi$. See Appendix (ref) for expressions of general moments $2<r<\nu_{\min}$.
An interesting convolution is that of a Gaussian distributed random variable and an independent Cauchy distributed random variable.\footnote{The convolution of a Gaussian random variable and a $t$-distribution with odd degrees of freedom is analyzed in Nason:2006, and Forchini:2008 generalized this analysis to $t$-distributions with any degrees of freedom.} Suppose that $Z\sim N(0,1)$ and $X\sim\text{\ensuremath{\mathrm{Cauchy}(0,1)}}$ are independent, and we seek the distribution of their convolution, $Y=Z+X$. The resulting distribution is known as the Voigt profile in the field of spectroscopy, and Kendall_D:1938 derived an expression involving the complementary error function.\footnote{Impressively, David G. Kendall published this result as a second year undergraduate student at Oxford University, see Bingham:1996.} We can apply Theorem (ref) to obtain the expression for the density in Kendall_D:1938.
The characteristic functions for $Z$ and $X$ are $\varphi_{Z}(s)=e^{-s^{2}/2}$ and $\varphi_{X}(s)=e^{-|s|}$, respectively, such that $\varphi_{Y}(s)=e^{-s^{2}/2-s}$, for $s\geq0$. By Euler's formula, the expression for the density in Theorem (ref) simplifies to \[ f_{Y}(y)=\frac{1}{\pi}\int_{0}^{\infty}\cos(sy)e^{-s^{2}/2-s}\mathrm{d}s, \] and with some additional simplifications we arrive at the following expression for $f_{Y}(y)$.
The analytical expression in ((ref)) is very accurate and fast to evaluate. The scaled complementary error function can be expressed as ${\rm erfcx}(z)=\exp(z^{2})[1-{\rm erf}(z)]$, where ${\rm erf}(z)=\frac{2}{\sqrt{\pi}}\int_{-\infty}^{z}e^{-t^{2}}\mathrm{d}t$ is the error function. These are standard non-elementary functions that are implemented in software, such as Julia, MATLAB, and R. Interestingly, the function, ${\rm erfcx}$, is also related to the Mills ratio that appears in other econometric problems, such as the Heckman model, see Heckman:1979. This follows from \[ {\rm erfcx}\left(\tfrac{u}{\sqrt{2}}\right)=e^{u^{2}/2}\frac{2}{\sqrt{\pi}}\int_{-\infty}^{u/\sqrt{2}}e^{-t^{2}}\mathrm{d}t=\frac{2}{\sqrt{\pi}}\frac{\tfrac{1}{\sqrt{2}}\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{u}e^{-s^{2}/2}\mathrm{d}s}{\frac{1}{\sqrt{2\pi}}e^{-u^{2}/2}}=\sqrt{\frac{2}{\pi}}\frac{1-\Phi(u)}{\varphi(u)}, \] where we here used $\Phi(u)$ and $\varphi(u)$ to denote the cumulative distribution function and the density, respectively, for a standard normal random variable.
The expression for the density of a convolution-$t$ distribution is typically complicated. An exception is a convolution of $t$-distributions with odd degrees of freedom, which has a much simpler expression. The reason is that the characteristic function for a $t$-distribution with odd degrees of freedom have the simple form, $\phi_{\nu}(s)=e^{-s_{\nu}}p(s_{\nu})$, where $s_{\nu}=\sqrt{\nu}|s|$ and $p(s_{\nu})$ is a polynomial of order $(\nu-1)/2\in\mathbb{N}$. The first few are given by, $\phi_{1}(s)=e^{-|s|}$, $\phi_{3}(s)=e^{-\sqrt{3}|s|}(1+\sqrt{3}|s|)$, $\phi_{5}(s)=e^{-s_{5}}(1+s_{5}+\tfrac{1}{3}s_{5}^{2})$, and $\phi_{7}(s)=e^{-s_{7}}(1+s_{7}+\tfrac{6}{15}s_{7}^{2}+\tfrac{1}{15}s_{7}^{3})$.\footnote{The highest-order term is $s_{\nu}^{m}/(\nu-2)!!$, where $m=(\nu-1)/2$. E.g. for $\nu=13$ this term is $s_{13}^{6}/10395$. Also note that $(\nu-2)!!=(2m-1)!!=2^{m}\Gamma(m+\tfrac{1}{2})/\Gamma(\tfrac{1}{2})$. } These simplifications explain that the most detailed results for convolutions involve odd degrees of freedom, see e.g. FisherHealy:1956, WalkerSaw:1978, FanBerger1990, NadarajahDey:2005, and Nason:2006.
For $\nu_{1}=1$ and $\nu_{2}=3$ we have $\varphi_{Y_{1}}(s)=(1+\sqrt{3}|s|)e^{-isy-(1+\sqrt{3})|s|}$ and from \[ \int_{-\infty}^{\infty}(1+\sqrt{3}|s|)e^{-isy-(1+\sqrt{3})|s|}ds=\frac{iy+2\sqrt{3}+1}{(iy+\sqrt{3}+1)^{2}}, \] we have that \[ f_{t_{1}+t_{3}}(y)=\frac{1}{\pi}\mathrm{Re}\left[\frac{iy+2\sqrt{3}+1}{(iy+\sqrt{3}+1)^{2}}\right]=\frac{1}{\pi}\frac{y^{2}+16+10\sqrt{3}}{\left(y^{2}+4+2\sqrt{3}\right)^{2}}. \] Similarly, we can obtain \[ f_{t_{1}+t_{5}}(y)=\frac{1}{\pi}\frac{y^{4}+2(11+3\sqrt{5})y^{2}+\tfrac{8}{3}(131+61\sqrt{5})}{\left(y^{2}+2(3+\sqrt{5})\right)^{3}}, \] and other convolutions of $t$-distributions with odd degrees of freedom lead to similar expressions that are ratios of polynomials in $y$. The expressions presented above are simpler than those in NadarajahDey:2005.\footnote{For instance, NadarajahDey:2005 has $f_{t_{1}+t_{3}}(y)=\frac{1}{\pi}(y^{2}+16+10\sqrt{3})(y^{2}+4+2\sqrt{3})^{2}(y^{2}+4-2\sqrt{3})^{4}/(4+8y^{2}+y^{4})^{2}$ and $f_{t_{1}+t_{5}}(y)=\frac{1}{\pi}(488\sqrt{5}+1048+(66+18\sqrt{5})y^{2}+3y^{4})(y^{2}+6+2\sqrt{5})^{3}(y^{2}+6-2\sqrt{5})^{6}/[5(16+12y^{2}+y^{4})^{8}]$, which can be verified to be identical to our expressions.}
We address identification, estimation, and inference of convolution-$t$ distributions in this section. Thus, consider $Y=\mu+\Xi X=\mu+\sum_{k=1}^{K}\Xi_{k}X_{k}$, where $\Xi=(\Xi_{1},\ldots,\Xi_{K})$ with $\Xi_{k}\in\mathbb{R}^{n\times n_{k}}$, and define $\Omega_{k}=\Xi_{k}\Xi_{k}^{\prime}$, for $k=1,\ldots,K$.
Throughout this section we make the following assumption about $\Xi\in\mathbb{R}^{n\times n}$.
It is helpful to consider the case, $K=1$, where $X$ has a multivariate $t$-distribution. In this case $\tilde{X}=Q^{\prime}X$ has the same distribution as $X$, for any orthonormal $n\times n$ matrix, $Q$. Thus, if we set $\tilde{\Xi}=\Xi Q$, then $\tilde{\Xi}X=\Xi QQ^{\prime}X$ has the same distribution as $\Xi X$. Given any $\Xi$ that satisfies Assumption (ref), then $Q=\Xi^{\prime}(\Xi\Xi^{\prime})^{-\frac{1}{2}}$ satisfies $Q^{\prime}Q=I_{n}$. With this choice we have $\tilde{\Xi}=\Xi Q=(\Xi\Xi^{\prime})^{\frac{1}{2}}$, which shows that any $\Xi$-matrix will have an observationally equivalent symmetric $\Xi$-matrix, in the special case where $K=1$. For $K>1$ the situation is different.
Identification is somewhat simplified if all degrees of freedom are unique, such that they can be ordered as $\nu_{1}<\cdots<\nu_{K}$. However, this rules out important convolutions of $t$-distributions with the same degrees of freedom. One exception is the Gaussian case, $\nu=\infty$, because we can always combing multiple Gaussian subvectors into a single subvector. Thus, the following assumption is without loss of generality.
That Assumption (ref) is without loss of generality follows from the a simple argument. Suppose $X_{i}\sim N(0,I_{n_{i}})$ and $X_{j}\sim N(0,I_{n_{j}})$ then $\sum_{k}\Xi_{k}X_{k}$ and $\sum_{k\neq i,j}\Xi_{k}X_{k}+\tilde{\Xi}_{m}\tilde{X}_{m}$ have the same distribution if we set $\tilde{\Xi}_{m}=(\Xi_{i},\Xi_{j})$ and $\tilde{X}_{m}\sim N(0,I_{n_{i}+n_{j}})$.
It is, however, not possible to combine $t$-distributed subvectors with identical degrees of freedom into a single subvector. Even if the subvectors are Cauchy distributed ($\nu_{k}=1)$, which is the only $t$-distribution that is a stable distribution (as is the Gaussian distribution). To see that there is not a way to merge two independent Cauchy distributions into a single distribution. Suppose that $X_{1}$ and $X_{2}$ are independent and Cauchy distributed. The question is if there is a way to express $\Xi_{1}X_{1}+\Xi_{2}X_{2}$ as $\tilde{\Xi}\tilde{X}$, such that their distributions are identical, where $\tilde{X}$ is an $n_{1}+n_{2}$ dimensional $t$-distribution. We have $\tilde{\Xi}\tilde{X}\sim\mathrm{Cauchy}(0,\tilde{\Xi}\tilde{\Xi}^{\prime})$ and from Bian1991 it follows that $\Xi_{1}X_{1}+\Xi_{2}X_{2}$ is Cauchy distributed if and only if $\Xi_{1}\Xi_{1}^{\prime}\propto\Xi_{2}\Xi_{2}^{\prime}$. However, this proportionality implies that $\Xi$ has reduced rank, which is a violation of Assumption (ref).
Theorem (ref) shows that it is the $n\times n$ scaling matrices, $\Omega_{k}=\Xi_{k}\Xi_{k}^{\prime}$, $k=1,\ldots,K$ that hold the key to identification.
We seek convenient identifying assumptions that will lead to a unique $\Xi_{k}\in\mathbb{R}^{n\times n_{k}}$, which satisfies $\Omega_{k}=\Xi_{k}\Xi_{k}^{\prime}$. The parameter-matrix, $\Xi_{k}$, is partially identified the subspace spanned by the reduced-rank matrix, $\Omega_{k}$, and that $\mathrm{tr}\{\Xi_{k}^{\prime}\Xi_{k}\}=\mathrm{tr}\{\Omega_{k}\}$. Suppose that $\Xi_{k}$ is a solution to $\Omega_{k}=\Xi_{k}\Xi_{k}^{\prime}$. For any orthonormal $n_{k}\times n_{k}$ matrix $Q$, (i.e. $QQ^{\prime}=I_{n_{k}}$) we observe that $\tilde{\Xi}_{k}=\Xi_{k}Q$ is also a solution, because $\Omega_{k}=\tilde{\Xi}_{k}\tilde{\Xi}_{k}^{\prime}$, and $\mathrm{tr}\{\tilde{\Xi}_{k}^{\prime}\tilde{\Xi}_{k}\}=\mathrm{tr}\{Q^{\prime}\Xi_{k}^{\prime}\Xi_{k}Q\}=\mathrm{tr}\{\Xi_{k}^{\prime}\Xi_{k}QQ^{\prime}\}=\mathrm{tr}\{\Xi_{k}^{\prime}\Xi_{k}\}$. The following theorem gives our preferred identification scheme for $\Xi$ matrix.
From the distribution of $Y$ we can identify the triplets, $(n_{k},\nu_{k},\Omega_{k})$, while the ordering of the $K$ triplets can arbitrary. For instance, suppose that $\Xi(\Omega_{1},\ldots,\Omega_{K})$ has the structure stated in Theorem (ref), then reordering of the triplet will, initially, lead to a different $\Xi$ matrix, $\tilde{\Xi}=(\Xi_{j_{1}},\ldots,\Xi_{j_{K}})$, where $j_{1},\ldots,j_{K}$ is a permutation of $1,\ldots,K$. However, this matrix can always be manipulated into, $\Xi(\Omega_{j_{1}},\ldots,\Omega_{j_{K}})$, which has the diagonal-symmetry structure in Theorem (ref).
Each permutation of the $K$ clusters define a particular $\Xi$-matrix with ((ref)), which are all observationally equivalent. Fortunately, there is only a finite number of them, and we will therefore identify $\Xi$ by the permutation that maximizes the trace, $\mathrm{tr}\{\Xi\}$. This choice will often be meaningful, because it sorts the clusters in a way that align with the observed variables. We illustrate this in Example 1. This identification strategy will almost always suffice, because permutations that produce $\Xi$-matrices with identical traces are rare (they have Lebesgue measure zero in $\mathbb{R}^{n\times n}$). Another simple identification scheme is available if the pairs, $(n_{k},\nu_{k})$, $k=1,\ldots,K$, are unique. Here we may simple pick a ordering of the clusters, which identifies $\Xi$ with the structure in Theorem (ref). This is an alternative to using the trace of $\Xi$. The two may also be combined. For instance, we may sort the clusters by $n_{k}$, and sort clusters with the same size by maximizing the trace.
It is the independent components in $X$ that enables us to identify additional parameters in $\Xi$. This is similar to independent component analysis (ICA), which is commonly used in signal processing and machine learning, where it is known as blind source separation, see e.g. Stone2004. However, our distributional framework is more general, because we do not require $n$ mutually independent components. To the contrary, elements of the same subvector, $X_{k}$, are dependent, and this dependence is informative about the underlying cluster structure.
There are other identification strategies than that of Theorem (ref). For instance, one could impose a Cholesky structure on $\Xi_{kk}$. We prefer the structure with symmetric diagonal blocks, because it is well suited for the parsimonious block structures we introduce below, as well as general symmetry restrictions that also reduce the number of free parameters.
This definition is taken from ArchakovHansen:CanonicalBlockMatrix. The case where the scaling matrix, $\Xi\Xi^{\prime}=\sum_{k=1}^{K}\Omega_{k}$, has a block matrix is interesting, because $\Xi$ may be assumed to have the same block structure, albeit not necessarily a symmetric block structure.
According to Theorem (ref), if $\Xi_{kk}$ is symmetric and positive definite, then the $\Xi_{k}$ matrix is unique and identified from $\Omega_{k}$. Thus, we can directly imposing the following structure of $\Xi$ matrix in estimation \[ \Xi=\left[\Xi_{1},\Xi_{2},\ldots,\Xi_{K}\right]=\left(
\right) \] where $\Xi_{kk}$ can by parameterized by $\Xi_{kk}=\exp(\gamma_{k})$, where $\gamma_{k}$ is an unrestricted symmetric matrix, and $\exp\left(\cdot\right)$ takes matrix exponential.
Note that if $\Omega_{k}$, $k=1,2,\ldots,K$, all have a block structure with identical block sizes, then $\Xi=\left[\Xi_{1},\Xi_{2},\ldots,\Xi_{K}\right]$ has the same block structure (but need not be symmetric). Block structures may be inferred from estimates of $\Omega=\sum_{k}\Omega_{k}$ that will have the same block structure, such as in our empirical analysis where $\hat{\Omega}$ is found to have an approximate block structures.
Next, we proceed to establish likelihood-based results that are used for inference about $\mu$, $\Xi$, and $\boldsymbol{\nu}$. We do not address inference about the integer-valued vector, $\boldsymbol{n}$, which is a non-standard problem, similar to lag-length selection in autoregressions and determining the cointegration rank in vector autoregressions.
We define the following three key parameter vectors with
where the $\tilde{\Xi}$ matrix is unrestricted (i.e. it need not have the identifying structure of Theorem (ref)). The corresponding score vectors are denoted by \[ \nabla_{\mu}=\tfrac{\partial\ell}{\partial\mu},\quad\nabla_{\tilde{\Xi}}=\tfrac{\partial\ell}{\partial{\rm vec}\left(\tilde{\Xi}\right)},\quad{\rm and}\quad\nabla_{\nu_{k}}=\tfrac{\partial\ell}{\partial\nu_{k}}, \] and, for later use, we define \[ W_{k}=\frac{\nu_{k}+n_{k}}{\nu_{k}+X_{k}^{\prime}X_{k}},\quad X_{k}=e_{k}^{\prime}X, \] where $e_{k}$ is a $n\times n_{k}$ matrix from identity matrix $I_{n}=\left(e_{1},e_{1},\ldots,e_{K}\right)$. We also introduce the notation $J_{k}\equiv e_{k}e_{k}^{\prime}$ and note that $I_{n}=\sum_{k=1}^{K}J_{k}$.
We have the following formulas for the score vectors, hessian matrix, and corresponding information matrix.
Note that in the expressions for the score and Hessian matrix, extreme values of $X_{k}$ are dampened by $W_{k}$. For instance, each of the terms $\nabla_{\mu}$, $\nabla_{\mu\mu^{\prime}}$, $\nabla_{\mu\nu_{k}}$, and $\nabla_{\nu_{k}\nu_{k}}$ are bounded, because $W_{k}\leq\frac{\nu_{k}+n_{k}}{\nu_{k}}$, $\left\Vert W_{k}X_{k}\right\Vert \leq\frac{1}{2}\sqrt{\frac{n_{k}}{\nu_{k}}}\left(\nu_{k}+n_{k}\right)$, and $W_{k}X_{k}^{\prime}X_{k}\leq\nu_{k}+n_{k}$, and it follows that all moments are finite for these terms. For the term, $\nabla_{\nu_{k}}$, we note that $\log\left(1+X_{k}^{\prime}X_{k}/\nu_{k}\right)$ has the same distribution as $\log(1/U)$, where $U\sim\mathrm{Beta}(\frac{\nu_{k}}{2},\frac{n_{k}}{2})$, for which all moments are finite,\footnote{$\mathbb{E}([\log(1/U)]^{m})=\psi^{(m-1)}(\nu_{k})-\psi^{(m-1)}(\nu_{k}+n_{k})$, where $\psi^{(n)}(x)$ is the $n$-th derivatives of $\log\Gamma(x)$, which are known as polygamma functions.} see Lemma (ref) for expressions of the first two moments. It can also be show that $\nabla_{\tilde{\Xi}}$ has finite moment when $K=1$, however this need not be the case in general. When $K\geq2$, there are interaction terms, such as $W_{k}X_{k}^{\prime}X_{l}$, which lacks $W_{l}$ to dampen the effect of large realizations of $X_{l}$. For the case, $K\geq2$, we show $\min_{k}\nu_{k}>1$ is needed to guarantee the existences of $\mathbb{E}\nabla_{\tilde{\Xi}}$, $\mathbb{E}\nabla_{\mu\tilde{\Xi}^{\prime}}$, and $\mathbb{E}\nabla_{\tilde{\Xi}\nu_{k}}$, and $\min_{k}\nu_{k}>2$ is needed for $\mathbb{E}\nabla_{\tilde{\Xi}\tilde{\Xi}^{\prime}}$ to be well-defined. The latter is intuitive since it involves the terms, $\frac{\nu_{k}}{\nu_{k}-2}$, $k=1,\ldots,K$.
Recall that $\Xi=\tilde{\Xi}P$, where $P=\mathrm{diag}(P_{11},\ldots,P_{KK})$ with $P_{kk}=\tilde{\Xi}_{kk}^{\prime}(\tilde{\Xi}_{kk}\tilde{\Xi}_{kk}^{\prime})^{-\frac{1}{2}}\in\mathbb{R}^{n_{k}\times n_{k}}$. From the results we established for the unstructured $\tilde{\Xi}$, we obtain the results for the identified parametrization, $\Xi$, by means of the Jacobian matrix of $\partial{\rm vec}(\Xi)$ with respect to $\partial{\rm vec}(\tilde{\Xi})$. Below, $K_{n,m}\in\mathbb{R}^{mn\times mn}$ denotes the commutation matrix.
Lemma (ref) is important for deriving the score and hessian matrix with respect to the identified $\Xi$. More specific, let $M_{\Xi}^{+}$ be the Moore-Penrose inverse of $M_{\tilde{\Xi}}$ matrix when evaluated at $\tilde{\Xi}=\Xi$, then we have
and $\nabla_{\Xi\mu^{\prime}}=M_{\Xi}^{+\prime}\nabla_{\tilde{\Xi}\nu\mu^{\prime}}$, $\nabla_{\Xi\nu^{\prime}}=M_{\Xi}^{+\prime}\nabla_{\tilde{\Xi}\nu^{\prime}}$. This also facilitate the computation of the two information matrices, given by \[ \mathcal{I}_{\Xi}=M_{\Xi}^{+\prime}\mathcal{I}_{\tilde{\Xi}}M_{\Xi}^{+},\quad\mathcal{J}_{\Xi}=M_{\Xi}^{+\prime}\mathcal{J}_{\tilde{\Xi}}M_{\Xi}^{+}. \] Note that the expression for $\mathcal{J}_{\Xi}$ does not requires $\frac{\partial{\rm vec}(M_{\Xi}^{+\prime})}{\partial{\rm vec}(\Xi)^{\prime}}$ to be evaluated because $\mathbb{E}[\nabla_{\tilde{\Xi}}]=0$. Also note that $\mathcal{I}_{\tilde{\Xi}}$ has reduced rank whereas $\mathcal{I}_{\Xi}$ is nonsingular, because $\Xi$ is identified.
Given a random sample of convolution-$t$ distributed random variables, $Y_{1},\ldots,Y_{T}$, the asymptotic properties of maximum likelihood estimators of $\theta=\{\mu,\Xi,\bm{\nu}\}$ are derived below.
When $K=1$, $\mathrm{CT}_{\boldsymbol{\nu},\boldsymbol{n}}(\mu,\Xi)$ simplifies to the multivariate $t$-distribution, for which many results exist and ZhuGalbraith:2010 established results for a univariate generalized $t$-distribution. Simulation-based evidence in the Supplemental Material suggests that the convergence rates for $\mu$, $\boldsymbol{\nu}$, and $\Xi_{kk}$ continues to be $\sqrt{T}$ in non-standard situations where $\min_{k}\nu_{k}\leq2$, whereas $\Xi$ coefficient in the off-diagonal blocks, $\Xi_{kl}$, $k\neq l$, have faster rates of convergence. The asymptotic normality of $\hat{\boldsymbol{\nu}}$ appears to be quite robust and also hold for small values of $\nu_{k}$, whereas the other parameters have more spiky limit distributions when $\min_{k}\nu_{k}\leq2$, that resemble Laplace distributions. This is particularly pronounced for the coefficients in the off-diagonal blocks of $\Xi$.
Note that the analytical expression of the Hessian matrix in Theorem (ref) and ((ref)) facilitate the computation of the sandwich-form of the asymptotic covariance matrix, $\frac{1}{T}\mathcal{\hat{J}}_{T}^{-1}\mathcal{\hat{I}}_{T}\mathcal{\hat{J}}_{T}^{-1}$, where $\mathcal{\hat{J}}_{T}=\frac{1}{T}\sum_{t=1}^{T}\nabla_{\hat{\theta}\hat{\theta}^{\prime},t}$ and $\mathcal{\hat{I}}_{T}=\frac{1}{T}\sum_{t=1}^{T}\left(\nabla_{\hat{\theta},t}\nabla_{\hat{\theta},t}^{\prime}\right)$. Note that $\hat{\theta}_{T}=\max_{\theta\in\Theta}\ell_{T}(\theta)$ implies $\frac{1}{T}\sum_{t=1}^{T}\nabla_{\tilde{\Xi},t}=0$ under our assumptions, such that the term $\frac{\partial{\rm vec}(M_{\Xi}^{+\prime})}{\partial{\rm vec}(\Xi)^{\prime}}$ is not needed in the computation of $\mathcal{\hat{J}}_{T}$.
We conduct a simulation study to verify the accuracy the asymptotic results for the maximum likelihood estimators. We consider the following trivariate system, given by \[ \mu=\left[
\right],\quad\Xi=\left[
\right]\quad\nu=\left[
\right] \] where the first two elements belong to one group, i.e. $n_{1}=1$ and $n_{2}=2$. And the structure of $\Xi_{11}$ and $\Xi_{22}$ are all symmetric and positive definite. Additional, because $n_{1}\neq n_{2}$, then according to Theorem (ref), the parameters $\theta=\left(\mu^{\prime},\nu^{\prime},{\rm vec}\left(\Xi\right)^{\prime}\right)^{\prime}$ are identified. We conduct a simulation study with $T=500,1000,2000,4000$, based on 50,000 Monte-Carlo simulations. We report the mean, standard deviation of the estimated parameters, $\alpha_{L}^{0.025}$ and $\alpha_{R}^{0.025}$ (defined below). The asymptotic standard deviations are also included for comparison, which are taken from the square root of the diagonal elements of the Cramèr-Rao bound, i.e., $\mathcal{I}_{\theta}^{-1}/T$. Note that we report the inverse of degree of freedoms, i.e. $1/\nu_{k}$, because its empirical distribution is more close to normal especially in small samples.
To examine the normality of asymptotic distribution, we reports the empirical left and right quantiles, defined as \[ \alpha_{L}^{0.025}=\frac{1}{T}\sum_{t=1}^{T}\left[\frac{\hat{\theta}_{t}-\theta}{{\rm aStd}\left(\theta\right)}<-1.96\right],\quad\alpha_{R}^{0.025}=\frac{1}{T}\sum_{t=1}^{T}\left[\frac{\hat{\theta}_{t}-\theta}{{\rm aStd}\left(\theta\right)}>1.96\right] \] and we will have $\alpha_{L}^{0.025}\rightarrow0.025$ and $\alpha_{R}^{0.025}\rightarrow0.025$ if $\hat{\theta}_{t}$ is approximately normally distributed. Table (ref) reports the simulation results in large and small samples, respectively. We can find that all maximum likelihood estimators $\hat{\theta}$ are consistent, and the empirical standard deviations match the asymptotic standard deviations very well. The empirical distributions of $\hat{\theta}$ converge to normal distribution when the sample size increases. In the Table (ref) in Appendix (ref), we also reports the estimation results in small samples.
The conventional multivariate $t$-distribution, $X\sim t_{n,\nu}(\mu,\Sigma)$ with $\nu>2$ has $\mathrm{var}(X)=\tfrac{\nu}{\nu-2}\Sigma$. Some theoretical results are more elegantly expressed in terms of standardized variables, $\tilde{X}=\mu+\sqrt{\tfrac{\nu-2}{\nu}}(X-\mu)$, which has $\mathbb{E}(\tilde{X})=\mu$, $\mathrm{var}(\tilde{X})=\Sigma$, and density
We will refer to this distribution as the standardized $t$-distribution and use the notation, $\tilde{X}\sim t_{n,\nu}^{\mathrm{std}}(\mu,\Sigma)$. Obviously, if $Y_{1}$ is a linear combination of independent $t$-distributed vectors, then it is also a linear combination of independent standardized $t$-distributed vectors if all $\nu_{k}>2$. Note that the standardized $t$-distribution disentangles the degrees of freedom parameter, $\nu$, from the variance. This is convenient in dynamic volatility models, such as the one we estimate in the empirical analysis. In our empirical analysis we will use the multivariate convolution-$t$ distribution, given by $Y=\mu+\Xi X$, where $X\in\mathbb{R}^{n}$ is composed of $K$ independent standardized multivariate $t$-distributions with $X_{k}\sim t_{n_{k},\nu_{k}}^{\mathrm{std}}(0,I_{k})$ and $n=\sum_{k}n_{k}$. This redefines $\Xi$ slightly, because elements are scaled by $\sqrt{\tfrac{\nu_{k}}{\nu_{k}-2}}$. We have the following notation for this standardized multivariate convolution-$t$ distribution \[ Y\sim\mathrm{CT}_{\boldsymbol{n},\boldsymbol{\nu}}^{\mathrm{std}}(\mu,\Xi), \] where $\boldsymbol{\nu}=(\nu_{1},\ldots,\nu_{K})^{\prime}$ with $\nu_{k}>2$ for all $k=1,\ldots,K$.
The marginal distribution of the univariate, $Y_{1}=\mu+\beta^{\prime}X$, where $X\sim\mathrm{CT}_{\boldsymbol{n},\boldsymbol{\nu}}^{\mathrm{std}}(0,I_{n}),$ and $\min_{k}\nu_{k}>4$, then convolution of $t$-distributions may be approximated by other distributions. In this section, we will approximate the distribution of $Y_{1}$ with a standardized $t$-distribution, which is either determined either by matching the first four moments (Method of Moments) or determined by minimizing the Kullback-Leibler divergence (KL divergence). For convolutions of two $t$-distributions, the method of moments approach was used in Patil:1965. Our result in Theorem (ref) makes it simple to generalize the results in Patil:1965 to convolutions of multiple $t$-distributions.
Let $g(\cdot)$ denote the true density function of $Y_{1}$, as derived in Theorem (ref), and let $f_{\mu,\sigma^{2},\nu}(y)$ denote the density of a standardized $t$-distribution, with mean $\mu$, variance, $\sigma^{2}$, and $\nu$ degrees of freedom. The first four moments of $Y_{1}$ were derived in Theorem (ref), and we will determine the $t$-distribution with the same first four moments. The third centralized moment is zero for both distributions, and we match the match the first, second, and fourth moments by selecting $\mu$, $\sigma^{2}$, and $\nu$, and this is achieved with $\mu_{\star}=\mu$, $\sigma_{\star}^{2}=\sum_{k=1}^{K}\beta_{k}^{\prime}\beta_{k}$, and $\nu_{\star}=4+\frac{6}{\kappa_{Y_{1}}},$ such that the moment-matching density is given by,
Alternatively we can determine $\mu$, $\sigma^{2}$, and $\nu$ by minimizing the Kullback-Leibler discrepancy between $f_{\mu,\sigma^{2},\nu}(y)$ and $g$. The Kullback-Leibler discrepancy is given by \[ \ensuremath{{\rm KLIC}(g,f)=\int_{-\infty}^{+\infty}\log\left(\frac{f(y)}{g(y)}\right)g(y)\mathrm{d}y}, \] and by solving \[ (\mu_{\ast},\sigma_{\ast}^{2},\nu_{\ast})=\arg\max_{\mu,\sigma^{2},\nu}{\rm KLIC}(g,f_{\mu,\sigma^{2},\nu}), \] we obtain the best approximating density, $f^{\ast}(y)=f_{\mu_{\ast},\sigma_{\ast}^{2},\nu_{\ast}}(y)$, in terms of the $\mathrm{KLIC}$.
In the left panels of Figure (ref), we consider the case where the true density is a convolution of two independent $t$-distributions, with the same degrees of freedom, and in the right panels we consider the case where the true density is a convolution of a normal distribution and an independent $t$-distribution. In the upper panels, (a) and (b), we report how well $f^{\star}$ (line) and $f^{\ast}$ (line) approximates the true density for $\nu>4$, in units of logarithm KLIC. In both cases it is evident that the approximating densities differ from the true density. Thus, a convolution $t$-distribution is not a $t$-distribution. The upper panels also show that the method of moments density is different from the KLIC minimizing density, especially for small values of $\nu$. And their differences decrease when $\nu$ becomes larger. However, one should note that the difference are small as the scale is in logarithmic units of KLIC.
In Panels (c) and (e) we plot the approximating densities (lines) and the true density (shaded area) of $Y=X_{1}+X_{2}$ with $X_{1},X_{2}\sim t_{\nu}^{{\rm std}}(0,1)$ for $\nu=7$ and $\nu=4.01$. Similarly in Panels (d) and (f) we plot the approximating densities (lines) and the true density (shaded area) of $Y=Z+X$ with $Z\sim N(0,1)$ and $X\sim t_{\nu}^{{\rm std}}$ with $\nu=7$ and $\nu=4.01$, respectively. Panels (c) and (d) show that the approximating densities (lines) match the true density very well when $\nu\geq7$. However, when $\nu$ become smaller and approaches $4$, Panel (e) and (f) shows that there exit significant differences between the approximating density from method of moments and true density. As for the $f^{\ast}$ (line) from minimizing KLIC, it still works in the sum of two $t$-distribution with same degree of freedom $\nu=4.01$, but it performs poor for the sum of a normal and $t$-distribution with $\nu=4.01$.
We note that the method of moments method is relatively poor, when the convolutions involve $t$-distributions with small degrees of freedom.
Realized measures of volatility are typically found to be very persistent and HAR model by Corsi:2009 (and variations thereof) are often used to model realized measures, see e.g. Andersen2007, CorsiFusariVecchia2013, and Bollerslev2016.
In this section, we model a vector of daily logarithmically transformed realized kernel estimators, see BNHLS-RK:2008\nocite{BNHLS-MRK:2011}. We use a HAR specification to model the conditional mean of the transformed realized measures and use convolution-$t$ distributions to model the errors in the model.
The sample spans the period from January 1, 2001 to December 31, 2020, which includes $T=4,972$ trading days and we include $n=10$ assets in stock markets in our analysis. These consist of S&P 500 index (SPX) presenting market portfolio, GE, HOM, and MMM, from the Industrials sector, JNJ, LLY, and MRK, from the Health Care sector, and AAPL, IBM, and INTC, from the Information Technology sector. We compute daily realized kernel estimator using intraday transaction data, where the high-frequency data of SPX were downloaded from TickData.com and the remaining nine individual stocks were obtained from the TAQ database. These were cleaned following the methodology in BNHLS-RKpractice:2009, and we used the Parzen kernel to compute the realized kernel estimates.
Let $Y_{i,t}=\log RK_{i,t}$ where $i=1,\ldots,n$, where $i$ indexes the stocks and $t=1,\ldots,T$ the day. The righthand side variables in the HAR model are defined by: \[ Y_{t}^{(d)}=Y_{t-1},\quad Y_{t}^{(w)}=\frac{1}{4}\sum_{i=2}^{5}Y_{t-i},\quad\text{and}\quad Y_{t}^{(m)}=\frac{1}{17}\sum_{i=6}^{22}Y_{t-i}, \] which correspond to daily, weekly, and monthly frequencies, respectively. Thees are defined such that no lagged $Y_{t-i}$ is used twice, which makes their coefficients easier to interpret. Our HAR model is given by
where $\xi\in\mathbb{R}^{n}$, $\beta_{d}$, $\beta_{w}$, and $\beta_{m}$ are diagonal matrices. From $V_{t}\sim\mathrm{CT}_{\boldsymbol{n},\boldsymbol{\nu}}^{\mathrm{std}}(0,\Xi)$ it follows that the conditional covariance of $Y_{t}$ is given by $\Sigma=\Xi\Xi^{\prime}$, and we can express each elements of $Y_{t}$ as an univariate HAR model
where $\beta_{d,i}$ is the $i$-th diagonal element of $\beta_{d}$, and $\beta_{w,i}$ and $\beta_{m,i}$ defined similarly.
Let $U_{t}=\Xi^{-1}V_{t}$ be the standardized residuals, and let $\mathcal{F}_{t}$ be the natural filtration, then the conditional log-likelihood function $Y_{t}|\mathcal{F}_{t-1}$ is,
The expression for the marginal density of convolution $t$-distributions in Theorem (ref) is particularly useful in our empirical analysis, because it facilitates a factorization of the joint density into marginal densities and the copula density by Sklar's theorem. This leads to the following decomposition of the log-likelihood, \[ \ell(Y_{t})=\sum_{i=1}^{n}\ell(Y_{t,i})+\log\left(c\left(Y_{t}\right)\right), \] where $c(Y_{t})$ denotes the copula density. This factorization can be used to compare different model specifications. For instance, if one specification has a larger value for the log-likelihood, then the factorization can be used to investigate if the gains are driven by gains in the marginal distributions, by gains in the copula density, or a combination of the two. It may reveal that some model features important for the marginal distribution, whereas other features are more specific to the copula.
Motivated by the diagonal-block structure of the information matrix shown in Theorem (ref), we adopt a two-sage estimation method. In the first stage, we estimate the regression coefficients, $\beta_{d,i}$, $\beta_{w,i}$, and $\beta_{m,i}$, by least squares for each of the ten assets separately and stack their residuals into the vector $\hat{V}_{t}$. In the second stage, we fit parametric distributions to $\hat{V}_{t}$, by maximizing the corresponding log-likelihood function, $\ell=\sum_{t=1}^{T}\ell_{t-1}\left(Y_{t}\right)$ in ((ref)).
The multivariate convolution-$t$ distribution, $\mathrm{CT}_{\boldsymbol{n},\boldsymbol{\nu}}^{\mathrm{std}}(0,\Xi)$, is characterized by the group assignments, $\bm{n}$, the vector of degrees of freedoms, $\bm{\nu}$, and a scale-rotation matrix, $\Xi$. In this paper, we consider three types of group assignments $\bm{n}$. The first time is the single cluster case, $K=1$ (i.e. with $\bm{n}=n=10$), which is simply a multivariate $t$-distribution (or Gaussian). The second case has $K=4$, with $\bm{n}=\left(1,3,3,3\right)$, which matches the sector classification of the 10 assets, and use “Cluster-$t$” to refer this structure. The third case, has $K=10$ independent components (i.e. no cluster structure) and $\bm{n}=\left(1,1,\ldots,1\right)$. We use “Hetero-$t$” as the label for this structure.
For the scale-rotation matrix, $\Xi$, we consider both the just-identified asymmetric structure and the over-identified symmetric structure. We also combine the symmetric and asymmetric $\Xi$-matrices with block structure, which has the label “Block”, whereas $\Xi$-matrices without an imposed block structure are labelled as “Just-Identified”. Imposing block structure on $\Xi$ is motivated by the approximate block structure seen in the empirical conditional covariance matrix of $Y_{t}$ (shown below). One should remember that, for the Gaussian case with $\bm{\nu}\rightarrow\infty$ and the conventional multivariate-$t$ case with $K=1$, only the covariance matrix $\Sigma=\Xi\Xi^{\prime}$ is identified.
We present the results for the HAR regression coefficients from first-stage estimation in Table (ref). The HAR coefficients are reported along with their corresponding standard errors in parentheses. For the logarithmic realized kernel estimators, $Y_{i,t}$, $i=1,\ldots,10,$ we report the estimated values of the model-implied expected values, $\mathbb{E}[Y_{t,i}]$, their persistence parameter, $\pi_{i}=\beta_{d,i}+\beta_{w,i}+\beta_{m,i}$, and the standard deviation of the residuals, $\sigma_{u_{i}}$. There are some variations in the levels of log-volatility, $\mathbb{E}[Y]$, across assets, as one would expect. The estimated HAR regression coefficients and conditional standard deviations, $\hat{\sigma}_{1},\ldots,\hat{\sigma}_{9}$, are very similar for all stocks, whereas the SPX is somewhat different. It has a higher weight on the first lag and less weight on distant lags as defined by $\beta_{m}$. The average log-volatility is also substantially smaller for SPX, implying that the volatilities for the individual stocks are substantially larger.\footnote{The log-differences of between 0.66 and 1.32 translate to variances that are between 92% and 275% larger or, alternatively, volatilities that are 39%-94% larger.}
From the first-stage residuals, $\hat{V}_{t}$, $t=1,\ldots,T$, we compute the empirical covariance matrix, $\hat{\Sigma}=\frac{1}{T}\sum_{t=1}^{T}\hat{V}_{t}\hat{V}_{t}^{\prime}$, which we present in Table (ref). The sector-based block structure is illustrated with the shaded regions. Interestingly, the covariance matrix has an approximate block structure that coincides with the sector classification of the stocks. The within-sector covariances are much higher than between-sector covariances, and we will explore this structure in Section (ref).
Next, we report the second stage estimation results for the just-identified $\Xi$ matrix, as defined by Theorem (ref). We maximize the log-likelihood for for the vectors of residuals, $\hat{V}_{t}$, $t=1,\ldots,T$, for each of the six distributional specifications. Table (ref) reports the estimates of the degrees of freedom vectors, $\boldsymbol{\nu}=(\nu_{1},\ldots,\nu_{K})^{\prime}$ while Table (ref) presents the estimated $\Xi$ matrix. In Table (ref), we also include the number of parameters $p$, the joint log-likelihood function $\ell$, the marginal log-likelihood $\ell_{m}$, the logarithmically transformed copula density, $\ell_{c}$, which is determined residually by $\ell_{c}=\ell-\ell_{m}$. We also report the Bayesian information criterion (BIC).
Several key observations emerge from Table (ref). There is strong evidence against the Gaussian specification. The multivariate-$t$ specification leads to a substantially better log-likelihood, resulting from large improvement in both the marginal distributions, $\ell_{m}$, and the log-copula density, $\ell_{c}$. The convolution-$t$ specifications, both the Cluster-$t$ and Hetero-$t$ specifications, lead to further improvements over the multivariate $t$-distribution. Once again we observe large improvements in the log-likelihood, with the largest gains originating from the log-copula density. There is some heterogenous values in the estimated degrees of freedom. For instance, the log-volatility related to the Health Care sector shows the heaviest tails. Both the Cluster-$t$ and the Hetero-$t$ distributions, suggest that a symmetric $\Xi$-matrix is not supported by the data. It is too restrictive to capture the complex nonlinear dependences, which the just-identified structure can accommodate. This is supported by the large improvements in the log-copula densities $\ell_{c}$, whereas marginal densities $\ell_{m}$ are slightly lower with the just-identified $\Xi$. In terms of the total log-likelihood the Hetero-$t$ specification is inferior to the Cluster-$t$ structure, despite having more free parameters.
The Hetero-$t$ specifications have ten degrees of freedom parameters whereas the Cluster-$t$ specifications have four. While the six additional parameters result in slightly better marginal densities, $\ell_{m}$, the copula densities, $\ell_{c}$, are actually worse for the Hetero-$t$ specifications. Recall that the Cluster-$t$ specification is not nested in the corresponding Hetero-$t$ specification, because the former imply non-linear dependencies between variables in the same cluster variable, $X_{k}$, whereas the Hetero-$t$ structure implies that all elements of $X$ are independent. The empirical results strongly indicate that the cluster structure captures important nonlinear dependences, because the log-copula densities are larger for the Cluster-$t$ specifications. Fifth, the empirical evidence also favors the (just identified) asymmetric structure for $\Xi$, over the (over identified) symmetric structure. Despite having 36 and 45 additional parameters, the BIC favors the asymmetric structure for both types of convolution-$t$ distributions.
Table (ref) reports the estimated $\Xi$ matrix for each of the convolution-$t$ specifications in Table (ref). The results for just-identified asymmetric $\Xi$-matrices are shown in Panels (a) and (b) and the results for over-identified symmetric $\Xi$-matrices are shown in Panels (c) and (d). Each of the four estimates of $\Xi$ have an approximate block structures, which we have highlighted with the (lighter) shaded regions. A darker shade is used to highlight the largest elements of $\hat{\Xi}$. The asymmetric structures in Panels (a) and (b) are particularly interesting, because they indicate a factor structure, where all assets load on a common (market) factor, where as the individual nine stocks also have large loadings on a single distinct variable. This factor structure is similar to the traditional CAPM model, but it is the structure implied by convolution-$t$ distributions, that enables us to identify this factor structure. Second, we also notice that the loading coefficients are very similar within each cluster, which motives the block-matrix structure on $\Xi$, which we explore below. Panel (c) and (c) present the estimates for over-identified symmetric $\Xi$ matrices. All coefficients are estimated to positive with the diagonal of $\hat{\Xi}$ being the dominant elements. For this specification, we can therefore largely link the degrees of freedom parameter of an element of $X$ to the corresponding elements of $Y$. The approximated block structure is also very much evident for the symmetric $\Xi$-matrices.
Next, we imposed a block structure on $\Xi$, which reduces the number of free parameters in $\Xi$. By imposing block structure, the free parameters in $\Xi$ would greatly decrease which only depends on the number of $K$, which is at most $K\left(K+1\right)$ if number of elements in each group is at least two.\footnote{When there are $\tilde{K}\leq K$ groups with only one element, this number become $K\left(K+1\right)-\tilde{K}$. The reason for the distinction between these two cases is that an $1\times1$ diagonal block has only one coefficient.} If additional asymmetry is imposed, then this number becomes $K\left(K+3\right)/2$.
Table (ref) presents the estimation results of degrees of freedom under block $\Xi$ matrix. When compared with the results without block structure in Table (ref), we find that: first, the estimates of degrees of freedom are very similar to these in Table (ref), and the findings about the log-likelihood values maintain. Second, importantly, the BIC values are much smaller than these with Just-Identified structure, which means the specification with block structure is not too restrictive and is preferred according to BIC. The estimated $\Xi$ matrix is provided in Table (ref), and we can find that the factor patterns are robust under block structure.
We have introduced the convolution-$t$ distributions, which is a versatile class of multivariate distributions for modeling heavy-tailed data with complex dependencies. We have characterized a several properties of these distributions and detailed results for the marginal distributions of convolution-$t$ distributions, such as expressions their densities and cumulative distribution functions. We obtained simple expressions for the first four moments and used these to develop an approximation method for convolution-$t$ distributions.
A attractive feature of this class of distributions is that their log-likelihood function have simple expressions. This makes estimation and inference relatively straight forward. We have analyzed the identification problem, which motivates a particular parametrization, and we established consistency and asymptotic normality of the maximum likelihood estimator in Theorem (ref),
Our analysis of the dynamic properties of ten realized volatility measures highlights the empirical relevance of the new class of distributions. The empirical results obtained with convolution-$t$ distributions provide a far more nuanced understanding of the nonlinear dependencies in these time series. The convolution-$t$ distributions improve the empirical fit, with improvements seen in both marginal distributions and their interdependencies, where the latter is characterized by the copula density. The specifications labelled Cluster-$t$ and Hetero-$t$ both offered substantial improvements in the empirical fit.
We identified interesting cluster structures that align with the sector classifications of stocks. There were identified from non-linear dependencies that are made possible by the convolution-$t$ distributions. We found heterogenous degrees of freedom with stocks in the Health Care sector exhibiting the heaviest tails, $\nu_{k}\approx5.96$, whereas stock in the Information Technology sector were about about $\nu_{k}\approx7.77$.
Conventional heavy-tailed distributions have also been found to be useful in multivariate GARCH models, see e.g. CrealKoopmanLucas:2011, and the convolution-$t$ distributions might prove to be useful in this context.