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.
99,249 characters · 19 sections · 0 citation commands
Testing error distribution by kernelized Stein discrepancy in multivariate time series models
\baselineskip=0.7 true cm
{\it Keywords and phrases}: Consistent test; Kernelized Stein discrepancy; Multivariate time series model; Testing multivariate error distribution.
Consider a multivariate stationary time series $\{Y_t\}$ with $Y_t=(Y_{1t},...,Y_{dt})^\top\in\mathbb{R}^d$, and $Y_t$ admits the following specification
where $I_t = \{Y_t, Y_{t-1}, ...\}$ is the information set up to time $t$, $\theta_0\in \mathbb{R}^p $ is the true yet unknown model parameter, $\eta_t \in \mathbb{R}^d $ is a sequence of independent and identically distributed (i.i.d.) errors with zero mean and identity covariance matrix $\mathrm{I}_d$, $M(\cdot;\theta_{0})\in \mathbb{R}^d$ is a known measurable vector function indexed by $\theta_{0}$, and $C(\cdot;\theta_0)\in\mathbb{R}^{d\times d}$ is a known measurable symmetric and positive definite matrix function indexed by $\theta_{0}$. Let $\mathcal{F}_{t} := \sigma(I_t)$ be a sigma-field generated by $I_{t}$. Conditional on $\mathcal{F}_{t-1}$, $M(I_{t-1};\theta_{0})$ and $C(I_{t-1};\theta_0)$ in ((ref)) are the conditional mean vector and conditional covariance matrix of $Y_t$, respectively. The general specification in ((ref)) covers many often used multivariate models including, for example, the vector autoregressive and moving-average (VARMA) model, the multivariate generalized autoregressive conditional heteroskedasticity (MGARCH) model, and their variants and combinations. For surveys on the multivariate time series models, we refer to L\"{u}tkepohl (2005), Bauwens et al. (2006), Tsay (2013), and Francq and Zako\"{i}an (2019).
For model ((ref)), $\eta_t$ is assumed to have certain continuous probability density function (p.d.f.) $p_0(x)$ in a myriad of applications, which include the validity of capital asset pricing model (Berk, 1997), the optimal forecasts (Christoffersen and Diebold, 1997), the density forecasts (Diebold et al., 1998), the interval forecasts (Zhu and Li, 2015), the option pricing (Zhu and Ling, 2015), and the Value-at-Risk and Expected Shortfall calculations (Taylor, 2019). However, the true p.d.f. of $\eta_t$, denoted by $p(x)$, is generally unknown in practice, and the empirical researchers could make wrong conclusions if their assumed p.d.f. $p_0$ is different from $p$. Motivated by this, it is important to testing for the following hypotheses
Our considered hypotheses in ((ref)) are designed for the unobserved model error $\eta_t$, which nests the observed data (i.e., $\eta_t=Y_t$) as a special case. In this paper, we mainly focus on the testing for unobserved $\eta_t$, and the testing methodologies for univariate/multivariate observed time series can be found in Lobato and Velasco (2004), Bai and Ng (2005), Mecklin and Mundfrom (2004), Sz\'{e}kely and Rizzo (2005), and the references therein.
Since $\eta_t$ is unobserved, one need use the model residual $\widehat{\eta}_t$ to form valid tests for the hypotheses in ((ref)). When $Y_t$ is univariate (i.e., $d=1$), a number of different testing methods were proposed in the literature. Bontemps and Meddahi (2005) considered the robust moment tests for normality of $\eta_t$ by using the Hermite polynomials, and their idea was further extended in Bontemps and Meddahi (2012) to examine the general distribution of $\eta_t$. Although these robust moment tests are easy-to-implement with a chi-square limiting null distribution, they are inconsistent as only a finite number of moments of $\eta_t$ are considered for the testing purpose. To construct consistent tests, other strategies have been adopted. For the general model as in ((ref)), Bai (2003) developed Kolmogorov--Smirnov (KS) and Cram\'{e}r--von Mises (CvM) tests by measuring the distance between the empirical distribution of $\widehat{\eta}_t$ and the cumulative distribution of $\eta_t$. For the GARCH model, Horv\'{a}th and Zitikis (2006) gave a smooth-type test by measuring the distance between the kernel density estimator of $p$ and the assumed density $p_0$ in $L_{\nu}$-norm with $1<\nu<\infty$, and Klar et al. (2012) constructed an integrated test by measuring the distance between the empirical characteristic function of $\widehat{\eta}_t$ and the characteristic function of $p_0$. For the ARMA--GARCH model, Koul and Ling (2006) studied a weighted KS test based on a vector of certain weighted residual empirical processes.
When $Y_t$ is multivariate (i.e., $d>1$) and both $M(\cdot;\theta_{0})$ and $C(\cdot;\theta_0)$ are constants, most of earlier efforts were made to detect the normality of $\eta_t$. See, for example, Mardia (1974), Henze and Zirkler (1990), Doornik and Hansen (2008), and references therein. When $Y_t$ is multivariate but either $M(\cdot;\theta_{0})$ or $C(\cdot;\theta_0)$ is non-constant, only few testing methods were provided for the MGARCH model. For instance, Bai and Chen (2008) applied a similar idea as Bai (2003) to propose consistent KS tests for detecting the multivariate normal and $t_{\nu}$ distributions of $\eta_t$. Their tests are asymptotically distribution-free, however, they are not fully consistent and require the explicit form of the conditional cumulative distribution function of $Y_{it}$ (conditional on $(Y_{1t},...,Y_{i-1,t})$), which is neither available for other multivariate distributions, nor easily computable for the dimension $d>2$. Francq et al. (2017) developed KS and CvM tests to examine whether $\eta_t$ has the elliptic distribution by extending the idea of Henze et al. (2014). These KS and CvM tests are not asymptotically distribution-free, and their application scope could be narrowed down when detecting the exact distribution of $\eta_t$ is needed. Henze et al. (2019) constructed consistent tests for the normality of $\eta_t$ by using the identity
where $R_{\eta}(v)$ is the real part of $\varphi_{\eta}(v)$, $\varphi_{\eta}(v)$ is the characteristic function of $\eta_t$, and $M_{\eta}(v)$ is the moment generating function of $\eta_t$. Since the identity above only holds for multivariate normal distributions, their idea can not be extended to testing for other distributions. In economic and financial applications, many heavy-tailed or skewed error distributions could outperform the multivariate normal or $t_{\nu}$ distribution (see, e.g., Haas et al. (2004), Bauwens and Laurent (2005), De Luca et al. (2006), and references therein). Hence, it is necessary to construct a valid test for detecting the general multivariate distribution of $\eta_t$ in model ((ref)).
This paper is motivated to propose a new consistent test for $H_0$ based on the kernelized Stein discrepancy (KSD) in Liu et al. (2016). The KSD measures the distance between the (Stein) score functions of $p$ and $p_0$ under the norm induced by a kernel function. For the observed data $\eta_t$, Liu et al. (2016) constructed a test statistic for $H_0$, and established its asymptotics. However, when $\eta_t$ is replaced by $\widehat{\eta}_t$, we find that their results are not applicable any more due to the estimation effect in $\widehat{\eta}_t$. To handle the estimation effect, our new KSD-based test is constructed based on a subsample of $\widehat{\eta}_t$. Under certain conditions, we show that our test has no estimation effect, and establish its asymptotics under $H_0$ and $H_1$. Although the estimation effect is negligible in theory, it may still exist in finite samples especially when the sample size is small. To overcome this difficulty, we introduce a simple parametric bootstrap method to calculate the critical values of our test. Simulations show that our test performs well in the examined cases, even when no or few effective data $\{\widehat{\eta}_t\}$ are discarded by our subsampling technique. A real data analysis is further given to demonstrate the usefulness of our test.
The remaining paper is organized as follows. Section 2 introduces the KSD-based test statistic. Section 3 studies the asymptotics of the KSD-based test statistic and provides a parametric bootstrap method to calculate the critical values. Simulation results are reported in Section 4, and a real example is offered in Section 5. Concluding remarks are given in Section 6. Proofs are deferred into Appendices.
In this paper, we construct a new test for hypotheses in ((ref)) based on the kernelized Stein discrepancy (KSD) in Liu et al. (2016). Let $p(x)$ be the true p.d.f. of $\eta_t$ in ((ref)) with the support $\aleph\subseteq \mathbb{R}^d$. To introduce the KSD, we first need define the (Stein) score function of $p$ and the Stein class of $p$.
When $\aleph=\mathbb{R}^d$, by using integration by parts, the condition ((ref)) holds if $$\lim_{\|x\|\to\infty}f(x)p(x)=0,$$ which holds, for example, if $p(x)$ is bounded and $\lim_{\|x\|\to\infty}f(x)=0$.
Next, let $k(x,x')$ be an integrally strictly positive definite kernel function, that is, $$\int_{x\in\aleph}\int_{x'\in\aleph}g(x)k(x,x')g(x')dxdx'>0$$ for any function $g(x)$ satisfying $0<\|g\|_{2}^{2}<\infty$. With the kernel function $k$, we are ready to give the definition of KSD between the distributions of $p$ and $p_0$.
Clearly, the KSD $\mathbb{S}(p,p_0)$ measures the difference between the (Stein) score functions of $p$ and $p_0$ under a norm induced by the kernel function $k$. If $p$ and $p_0$ are continuous with $\|p\delta_{p_0,p}\|_{2}^{2}<\infty$, Liu et al. (2016) showed that $\mathbb{S}(p,p_0)\geq 0$ and
In view of the result ((ref)), we can detect the null hypothesis $H_0$ in ((ref)) by examining whether $\mathbb{S}(p,p_0)$ is significantly different from zero. However, a direct testing implementation based on ((ref)) is infeasible, since the score difference $\delta_{p_0,p}$ is unknown. To overcome this difficulty, we need an additional condition on the kernel function $k$.
When the kernel function $k$ is in the Stein class of $p$, Liu et al. (2016) found that the KSD in ((ref)) becomes
where
Now, the formula of $\mathbb{S}(p, p_0)$ in ((ref)) is tractable for the testing purpose, since it only depends on the score function $s_{p_0}$ and the kernel function $k$, both of which are known under $H_0$.
To form our test statistic, a sample counterpart of $\mathbb{S}(p, p_0)$ in ((ref)), based on the model residuals, is needed. Let $\theta = (\theta_1, ..., \theta_p)^\top \in \Theta \subset \mathbb{R}^p $ be the unknown parameter of model ((ref)), where $\Theta$ is compact parametric space. Assume that $\theta_0$ is an interior point of $\Theta$, and denote
By ((ref)), the model residual in ((ref)) can be computed as
where $\widehat{I}_{t}$, containing possible given initial values, is the truncated information set at time $t$, and $\widehat{\theta}_n$ is an estimator of $\theta_0$. With model residuals $\{\widehat{\eta}_{t}\}_{t=1}^{n}$, the KSD-based test statistic $\widehat{\mathbb{S}}$ as the estimator of $\mathbb{S}(p, p_0)$ in ((ref)) is given by
where $n_0=[K_0n^{1-\varepsilon}]$ for some $K_0, \varepsilon>0$. Clearly, $\widehat{\mathbb{S}}$ is a U-statistic with kernel function $u(\widehat{\eta}_i, \widehat{\eta}_j)$, and the calculation of $\widehat{\mathbb{S}}$ only requires the computation of $s_{p_0}$, $\nabla_{x}k(x,x')$, $\nabla_{x'}k(x,x')$, and $\nabla_{x,x'}k(x,x')$, which does not raise any computational burden even for a large dimension $d$. For the kernel function $k$, the often used one is the Gaussian kernel
where $\sigma>0$ is a fixed constant; in this case, we have
For the score function $s_{p_0}$, we show how to calculate it for some well-known distributions.
Unlike Liu et al. (2016), our test statistic $\widehat{\mathbb{S}}$ does not use the entire data $\{\widehat{\eta}_{t}\}_{t=1}^{n}$. This is because we have to sacrifice some part of $\{\widehat{\eta}_{t}\}_{t=1}^{n}$ to deal with the effect of estimation uncertainty caused by replacing $\theta_0$ via $\widehat{\theta}_{n}$ and the effect of unobserved initial values resulting from substituting $I_t$ by $\widehat{I}_t$. With the assist of bootstrap scheme, our numerical studies in Section 4 show that $\widehat{\mathbb{S}}$ can have good a size and power performance even when no or few data $\{\widehat{\eta}_{t}\}_{t=1}^{n}$ are discarded. Hence, it suggests that $\widehat{\mathbb{S}}$ can be used with $n_0=n$ or $n_0\approx n$ in practice, and the subsampling technique seems only theoretically relevant.
Denote $$g_t(\theta)=g(Y_t, I_{t-1}; \theta),\,\,\,\widehat g_{t}(\theta)=g(Y_{t},\widehat I_{t-1};\theta),\,\,\,\mbox{ and }\,\,\,\widehat{R}_{t}(\theta)= \widehat g_{t}(\theta)-g_t(\theta),$$ where $g(Y_t, I_{t-1}; \theta)$ is defined in ((ref)). In this subsection, we give some technical assumptions to study the asymptotics of $\widehat{\mathbb{S}}$.
A few remarks are in order related to the aforementioned assumptions. Assumptions (ref)--(ref) are regular in many time series applications. Assumption (ref) poses some moment conditions on the derivatives of $g_t(\theta)$ for the purpose of proof, and Assumption (ref) holds for most estimators such as the least squares estimator (LSE) for VARMA models and the quasi-maximum likelihood estimator (QMLE) for VARMA--GARCH models. Sufficient conditions to validate Assumptions (ref)--(ref) can be found in L\"{u}tkepohl (2005) for VARMA models, Comte and Lieberman (2003), Hafner and Preminger (2009), and Francq and Zako\"{i}an (2012) for MGARCH models, and Ling and McAleer (2003) for VARMA--MGARCH models. Assumption (ref) is a condition on the approximation error by replacing the information set $I_t$ by $\widehat{I}_{t}$, and it is used to show that the unobserved initial values have the negligible effect on the asymptotic theory. See also Hong and Lee (2005) and Escanciano (2006) for the similar conditions.
Assumption (ref) requires both $p$ and $p_0$ to have certain smooth conditions. The condition $\|p\delta_{p_0,p}\|_{2}^{2}<\infty$ is sufficient to prove the equivalence result ((ref)). As argued in Liu et al. (2016), this condition is mild. For example, it holds when $p$ is the density function of multivariate normal and $t_{\nu}$ distributions or $p$ has an exponentially decayed tail, but it may not hold when $p$ has a heavy tail. Note that the exclusion of heavy-tailed $p$ is also implied by Assumption (ref). Assumption (ref)(i) ensures the validity of ((ref)), and Assumption (ref)(ii) poses some boundedness conditions on $k$ and its derivatives. It is easy to check that Gaussian kernel in ((ref)) satisfies Assumption (ref) for any smooth density $p$ supported on $\aleph=\mathbb{R}^d$. Hence, we follow Liu et al. (2016) to use the Gaussian kernel in this paper.
According to Theorem 3.7 in Liu et al. (2016), the kernel function $u(x,x')$ is positive definite, and then by Mercer's theorem, $u(x,x')$ admits the expansion
where $\{l_{m}(\cdot)\}$ and $\{\lambda_m\}$ are the orthonormal eigenfunctions and eigenvalues of $u(x,x')$. We are ready to give the limiting null distribution of $\widehat{\mathbb{S}}$.
Our limiting null distribution in Theorem (ref) is the same as the one in Theorem 4.1 of Liu et al. (2016), since the effects of estimation uncertainty and unobserved initial values are asymptotically negligible by using the sub-sample technique with $\varepsilon>1/2$. When $\varepsilon\leq1/2$, how to establish the limiting null distribution of $\widehat{\mathbb{S}}$ is unclear at this stage, and we leave this topic for future study.
Although the effects of estimation uncertainty and unobserved initial values are asymptotically negligible in theory, they may exist in finite samples especially when $n$ is small. To redeem this drawback, we propose a simple parametric bootstrap method in Subsection 3.3 below to calculate the critical values of $n_0\widehat{\mathbb{S}}$. Owing to the use of bootstrap, our simulation studies will show that $\widehat{\mathbb{S}}$ has a good finite-sample performance even for very small value of $\varepsilon$, indicating that the condition $\varepsilon>1/2$ should not be an obstacle for applications.
Next, the behavior of $\widehat{\mathbb{S}}$ under $H_1$ is given in the following theorem.
Let $c_{\alpha}$ be the critical value of $n_0\widehat{\mathbb{S}}$ at the level $\alpha$. Then, the preceding theorem implies that under $H_1$, the power function $\Lambda_{n}:=P(n_0\widehat{\mathbb{S}}>c_{\alpha})$ converges to 1 as $n\to\infty$, and hence $\widehat{\mathbb{S}}$ can detect $H_1$ consistently.
To end this subsection, we discuss how the choice of $\sigma$ in ((ref)) affects the value of $\Lambda_{n}$. By ((ref)) in Appendix A.2, we can show that under $H_1$, for large $n$,
To further calculate $\Lambda_{n}$, we assume $n_0\approx n$ (as recommended for practical use) and $$\Big(\sqrt{n}\big(\mathbb{S}^{(0)}-\mathbb{S}(p,p_0)\big), \sqrt{n}(\widehat{\theta}_n-\theta_0)^\top\Big)^\top \xrightarrow{d} N(0,\Sigma_{\mathbb{S},\theta})$$ as $n\to\infty$, where $\Sigma_{\mathbb{S},\theta}\in\mathbb{R}^{(q+1)\times(q+1)}$ is the asymptotic covariance matrix. Then, by ((ref)) it is straightforward to see
where $\kappa=\sqrt{(1,s_1^\top)\Sigma_{\mathbb{S},\theta}(1,s_1^\top)^\top}$, and $s_1$ is the limit of $\mathbb{S}^{(1)}$ by the law of large numbers for U-statistics. From ((ref)), we know that $\sigma$ should be chosen such that $\mathbb{S}(p,p_0)/\kappa$ is maximized. However, this implementation can not be accomplished in an easy way, since an explicit form of $\kappa$ is not available. Therefore, it seems hard to choose $\sigma$ optimally. In practice, we can follow Liu et al. (2016) to choose $\sigma$ as the median of residual distance:
where $\chi_{ij}={\rVert \widehat{\eta}_i - \widehat{\eta}_j \rVert}^2$. Our simulation studies in Section 4 below show that $\widehat{\mathbb{S}}$ has a good finite-sample performance based on this choice of $\sigma$.
When $\eta_t$ is observed (i.e., $\eta_t=Y_t$), Liu et al. (2016) applied a Wild bootstrap method to calculate critical values for their test. However, when $\eta_t$ is unobserved as in our settings, their bootstrap scheme may not work, since it does not account for the effects of estimation uncertainty and unobserved initial values, which can affect our critical value $c_{\alpha}$ in the finite sample. In this paper, we apply the following parametric bootstrap method to calculate $c_{\alpha}$:
Step 1. Draw bootstrap i.i.d. errors $\eta_{t}^{*}\sim p_0$ and calculate the bootstrap data sample $$Y_{t}^{*}=M(I_{t-1}^{*}; \widehat{\theta}_{n})+C^{1/2}(I_{t-1}^{*}; \widehat{\theta}_{n})\eta_{t}^{*},$$ where $I_{t-1}^{*}$ is the bootstrap counterpart of $I_{t-1}$.
Step 2. Calculate the bootstrap estimator $\widehat{\theta}_{n}^{*}$ and the bootstrap residuals $$\widehat{\eta}_{t}^{*}=g(Y_{t}^{*},\widehat{I}_{t-1}^{*};\widehat{\theta}_{n}^{*}),$$ where $\widehat{I}_{t-1}^{*}$ is the bootstrap counterpart of $\widehat{I}_{t-1}$.
Step 3. Compute the bootstrap test statistic $n_0\widehat{\mathbb{S}}^{*}$, based on the bootstrap residuals.
Step 4. Repeat steps 1--3 $m$ times to get $\{n_0\widehat{\mathbb{S}}_{(1)}^{*},...,n_0\widehat{\mathbb{S}}_{(m)}^{*}\}$, whose empirical $\alpha$ upper quantile is taken as the critical value $c_{\alpha}$.
The validity of $c_{\alpha}$ under $H_0$ and $H_1$ can be justified by using the similar arguments as for Theorems (ref) and (ref), respectively, and hence we omit the details.
In this section, we carry out simulation experiments to assess the performance of our KSD-based test $\widehat{\mathbb{S}}$ in finite samples. For the purpose of comparison, some widely-used tests (see Appendix A.3 for their definitions and asymptotics) are also considered. The data generating processes (DGPs) considered below cover the dimension $d=2$ and $5$. In all simulations, we take the sample size $n=100$ or $500$, choose the number of repetitions $J=10,000$, and set the significance level $\alpha= 1\%$, 5%, or 10%. For $\widehat{\mathbb{S}}$, we use the Gaussian kernel in ((ref)) with $\sigma$ taken as in ((ref)), and choose $n_0=n$ such that no data $\{\widehat{\eta}_t\}_{t=1}^{n}$ are discarded. To reduce the computational burden in simulations, we follow Francq et al. (2017) to adopt the Warp-Speed method of Giacomini et al. (2013) for evaluating the bootstrap scheme proposed in Subsection 3.3. With the Warp-Speed method, rather than computing critical value $c_{\alpha}$ for each repetition sample, only one resample is generated for each repetition sample and the resampling test statistic $\widehat{\mathbb{S}}^{*}$ is computed for that sample. Then the critical value $c_{\alpha}$ is computed from the empirical distribution determined by the resampling repetitions $\{\widehat{\mathbb{S}}^{*}_{(i)}\}_{i=1}^{J}$.
We consider the DGP given by a constant mean and constant covariance model
where $M$ and $C$ are constant mean and constant covariance of $Y_t$, respectively, and they are chosen as $$M = \left(
\right), \,\,\, C^{1/2} = \left(
\right)$$ for $d = 2$, and $$M = \left(
\right), \,\,\, C^{1/2} = \left(
\right)$$ for $d = 5$. In model (\ref{E4.2.1}), the distribution of $\eta_t$ (i.e., the true distribution $p$) is $N_d(0,\mathrm{I}_d)$, $T_d(5)$, $SN_d(\gamma)$, or $ST_d(5,\xi)$, where the first three distributions are given in Examples \ref{exam_1}--\ref{exam_3}, and the fourth distribution $ST_d(\nu,\xi)$ is the multivariate skew-$t$ distribution in Bauwens and Laurent (2005) with mean zero, covariance matrix $\mathrm{I}_d$, $\nu>2$ being the degrees of freedom, and $\xi=(\xi_{1},...,\xi_{d})^\top\in\mathbb{R}^{d}$ being the asymmetry vector to control the skewness. In the sequel, we set $$ \gamma= \left\{
\right. \,\,\, \xi= \left\{
\right. $$
For the null distribution $p_0$ in ((ref)), we take it to be $N_d(0, \mathrm{I}_d)$, $T_d(5)$, or $SN_d(\gamma)$. When $p_0$ is $N_d(0, \mathrm{I}_d)$, we also consider Mardia's skewness test ($\widehat{\mathbb{T}}_{M,1}$), Mardia's kurtosis test ($\widehat{\mathbb{T}}_{M,2}$), Doornik--Hansen test ($\widehat{\mathbb{T}}_{DH}$), Henze--Zirkler test ($\widehat{\mathbb{T}}_{HZ}$), Bai--Chen tests ($\widehat{\mathbb{T}}_{BC,1}$, $\widehat{\mathbb{T}}_{BC,2}$, and $\widehat{\mathbb{T}}_{BC,3}$), and Henze--Jim\'{e}nez-Gamero--Meintanis test ($\widehat{\mathbb{T}}_{HJM}$). The first four tests $\widehat{\mathbb{T}}_{M,1}$, $\widehat{\mathbb{T}}_{M,2}$, $\widehat{\mathbb{T}}_{DH}$, and $\widehat{\mathbb{T}}_{HZ}$ are based on either the sample skewness or the sample kurtosis or both. The tests $\widehat{\mathbb{T}}_{BC,i}$, $i=1,2,3$, are based on the empirical distribution of the residuals, and the test $\widehat{\mathbb{T}}_{HJM}$ is based on the characteristic function of the residuals. Note that except $\widehat{\mathbb{T}}_{BC,i}$, all other tests work for $d>2$. When $p_0$ is $T_d(5)$ or $SN_d(\gamma)$, none of the competitive tests above is applicable, except that the tests $\widehat{\mathbb{T}}_{BC,i}$ can be used for the case of $T_2(5)$.
Tables (ref) and (ref) report the size and power of all examined tests for $d=2$ and $d=5$, respectively, where the size corresponds to the case of $p=p_0$. In calculation of $\widehat{\mathbb{S}}$, $\widehat{\mathbb{T}}_{BC,i}$, and $\widehat{\mathbb{T}}_{HJM}$, the residuals of model ((ref)) are computed by estimating $M$ and $C$ by the sample mean and sample covariance of $Y_t$, respectively. Note that since the tests $\widehat{\mathbb{T}}_{BC,i}$ are largely over-sized, we compute their size-adjusted power in the sequel. From Tables (ref) and (ref), our findings are as follows:
(1) Except for the tests $\widehat{\mathbb{T}}_{BC,i}$, all examined tests have an accurate size performance at three levels.
(2) When $p_0$ is $N_d(0, \mathrm{I}_d)$, $\widehat{\mathbb{S}}$ has a comparative power performance with any competitive test to detect the alternative hypotheses that $p$ are $SN_d(\gamma)$ and $ST_d(5,\xi)$. However, $\widehat{\mathbb{S}}$ has the best power performance to detect the alternative hypothesis that $p$ is $SN_{d}(\gamma)$, and the tests $\widehat{\mathbb{T}}_{M,2}$, $\widehat{\mathbb{T}}_{BC,i}$, and $\widehat{\mathbb{T}}_{HJM}$ have a much worse power performance in this case. The advantage of $\widehat{\mathbb{S}}$ is more obvious for the case $d=5$.
(3) When $p_0$ is $T_d(5)$, $\widehat{\mathbb{S}}$ has the satisfactory power performance especially for $n=500$, while the tests $\widehat{\mathbb{T}}_{BC,i}$ only exhibit the power to detect the alternative hypothesis that $p$ is $ST_d(5,\xi)$.
(4) When $p_0$ is $SN_d(\gamma)$, $\widehat{\mathbb{S}}$ is powerful to detect each examined alternative hypothesis, and its power to detect the heavy-tailed alternative distribution (e.g., $T_d(5)$ or $ST_d(5,\xi)$) is higher than that to detect the light-tailed alternative distribution (e.g., $N_d(0, \mathrm{I}_d)$).
Overall, our KSD-based test $\widehat{\mathbb{S}}$ exhibits the good size and power performance in all examined cases. All skewness- or kurtosis-based tests for normality generally perform well, except that $\widehat{\mathbb{T}}_{M,2}$ lacks the power to detect the alternative distribution $SN_{d}(\gamma)$. The tests $\widehat{\mathbb{T}}_{BC,i}$ have the over-sized problem in all examined cases, and their size-adjusted power in general is not satisfactory especially for the null distribution $T_{2}(5)$. The test $\widehat{\mathbb{T}}_{HJM}$ for the normality performs as good as $\widehat{\mathbb{S}}$, except that its power to detect the alternative distribution $SN_d(\gamma)$ is lower. Based on the aforementioned findings, it is reasonable to recommend $\widehat{\mathbb{S}}$ for use due to its generality and desirable power performance.
We consider the DGP given by a VAR(3) model
where $M$, $C^{1/2}$, and $\eta_t$ are chosen as in model ((ref)), and $$A_1 =
,\quad A_2 =
,\quad A_3 =
$$ for $d = 2$, and
for $d = 5$. As in Case 1, the null distribution $p_0$ is $N_d(0,\mathrm{I}_d)$, $T_d(5)$, or $SN_d(\gamma)$. For the VAR(3) model in ((ref)), the skewness- or kurtosis-based tests considered in Case 1 are not applicable any more. In this case, the tests $\widehat{\mathbb{T}}_{BC,i}$ work when $p_0$ is $N_d(0,\mathrm{I}_d)$ or $T_d(5)$ for $d=2$, and the test $\widehat{\mathbb{T}}_{HJM}$ works when $p_0$ is $N_d(0,\mathrm{I}_d)$ for $d=2$ and 5.
Tables (ref) and (ref) report the size and power of all examined tests for $d=2$ and $d=5$, respectively, where the size corresponds to the case of $p=p_0$. In calculation of $\widehat{\mathbb{S}}$, $\widehat{\mathbb{T}}_{BC,i}$, and $\widehat{\mathbb{T}}_{HJM}$, the residuals of model ((ref)) are computed by using the LSE to estimate the unknown parameters. From Tables (ref) and (ref), our findings are similar as those in Case 1.
We consider the DGP given by a CCC-GARCH(1, 1) model
where $\eta_t$ is chosen as in model ((ref)), and $C_t=\mbox{diag}\{\sigma_{1,t},...,\sigma_{d,t}\} \cdot R\cdot \mbox{diag}\{\sigma_{1,t},...,\sigma_{d,t}\}$ with
Here, the parameter matrices $R$, $W$, $B$, and $\Gamma$ are set to be
for $d=2$, and
for $d=5$.
As in Cases 1 and 2, the null distribution $p_0$ is $N_d(0,\mathrm{I}_d)$, $T_d(5)$, or $SN_d(\gamma)$. For the CCC-GARCH model in ((ref)), the competitive tests can only be chosen as in Case 2. Tables (ref) and (ref) report the size and power of all examined tests for $d=2$ and $d=5$, respectively, where the size corresponds to the case of $p=p_0$. In calculation of $\widehat{\mathbb{S}}$, $\widehat{\mathbb{T}}_{BC,i}$, and $\widehat{\mathbb{T}}_{HJM}$, the residuals of model ((ref)) are computed by using the QMLE to estimate the unknown parameters. Clearly, our findings from Tables (ref) and (ref) are similar to those in Case 1.
In our previous simulation studies, we take $n_0=n$ and $\sigma$ as in ((ref)) to compute our KSD-based test $\widehat{\mathbb{S}}$. In this subsection, we implement the sensitivity analysis on the choice of $n_0$ or $\sigma$ for $\widehat{\mathbb{S}}$, based on the DGP in ((ref)) with $p$ being $N_d(0, \mathrm{I}_d)$ and $p_0$ being $N_d(0, \mathrm{I}_d)$ (for the size study) or $T_{d}(5)$ (for the power study).
First, we consider the cases that $n_0$ is taken with the subsample ratio $n_0/n=0.8$, $0.9$, $0.95$, and $1$, while the value of $\sigma$ is chosen as in ((ref)). Fig\,(ref) plots the size and power of $\widehat{\mathbb{S}}$ across the subsample ratio $n_0/n$. From this figure, we can find that (1) $\widehat{\mathbb{S}}$ always has a good size performance; (2) when $n=100$, the power of $\widehat{\mathbb{S}}$ increases as the value of $n_0/n$ (or $n_0$) increases, and when $n=500$, the power of $\widehat{\mathbb{S}}$ reaches one in all examined cases. Therefore, as expected, we should recommend to use $n_0=n$ for $\widehat{\mathbb{S}}$, although this choice of $n_0$ is inconsistent to our theoretical setting.
Second, we consider the cases that $\sigma$ is set to be 0.5, 0.7, ..., 3.1, while the value of $n_0$ is taken as $n$. Fig\,(ref) plots the size and power of $\widehat{\mathbb{S}}$ across $\sigma$. From this figure, we can find that the size of $\widehat{\mathbb{S}}$ is always accurate for each examined $\sigma$, and the power of $\widehat{\mathbb{S}}$ for $\sigma\geq 1.1$ has only a marginal difference from that for the choice of $\sigma$ in ((ref)). These findings imply that $\widehat{\mathbb{S}}$ tends to have a stable size and power performance over the choice of $\sigma$.
In this section, we revisit a real example in Tsay (2005). This example considered a three-dimensional financial time series, which consists of the daily log returns (in percentage) of the S&P 500 index, the stock price of Cisco Systems, and the stock price of Intel Corporation from January 2, 1991 to December 31, 1999, with 2275 observations in total. We denote this multivariate time series by $Y_t=(Y_{1t},Y_{2t},Y_{3t})^\top$, and plot each entry of $Y_t$ in Fig\,(ref). Following Tsay (2005), $Y_t$ is fitted by a VAR(3)--CCC--GARCH(1, 1) model
with
For model ((ref)), after dropping the insignificant parameters, we follow Tsay (2005) to first estimate the VAR(3) model by using the LSE, and then estimate the CCC--GARCH(1, 1) model by using the QMLE, where the resulting estimators are given by
Next, we use our KSD-based test $\widehat{\mathbb{S}}$ to check the distribution of $\eta_t$. The null distribution $p_0(x)$ of interest is $N_{3}(0, \mathrm{I}_3)$, $T_{3}(\nu)$, or $SN_{3}(\gamma)$, where the degrees of freedom $\nu$ is $\nu_{MLE}$, 6, 7, 8, or 9, and the skewness vector $\gamma$ is $\gamma_{MLE}$. Here, $\nu_{MLE}=7.724$ is the maximum likelihood estimator (MLE) of $\nu$ based on $\eta_t\sim T_{3}(\nu)$, and $\gamma_{MLE}=(-0.181, -0.023, 0)$ is the MLE of $\gamma$ based on $\eta_t\sim SN_{3}(\gamma)$. To calculate $\widehat{\mathbb{S}}$, we choose $n_0=n$ and use the Gaussian kernel $k$ in ((ref)) with $\sigma$ taken as in ((ref)). The p-value of $\widehat{\mathbb{S}}$ is computed based on the parametric bootstrap in Subsection 3.3 with $m=1000$.
Table (ref) reports the p-values of $\widehat{\mathbb{S}}$ for all chosen null distributions $p_0(x)$. From this table, we can find that $\widehat{\mathbb{S}}$ gives the strong evidence to reject the null distributions $N_{3}(0, \mathrm{I}_3)$ and $SN_{3}(\gamma_{MLE})$, and on the contrary, it can not reject the null distributions $T_{3}(\nu_{MLE})$, $T_{3}(7)$, $T_{3}(8)$, and $T_{3}(9)$ at the significant level 5%. Since $\widehat{\mathbb{S}}$ has the largest p-value for the null distribution $T_{3}(8)$, it is reasonable to conclude that $\eta_t$ in model ((ref)) follows $T_{3}(8)$.
This paper constructed a new KSD-based test to detect the error distribution in multivariate time series models with general specifications. The KSD-based test is easy-to-implement as long as the (Stein) score function of the null distribution has an explicit form. Hence, it allows the null distribution of interest to be not only multivariate normal, but also multivariate $t_{\nu}$, skew-normal, and many others. Since most of the existing tests only deal with the multivariate normal null distribution, the KSD-based test can largely broaden the testing scope for practitioners. This progress driven by the KSD-based test is important in view of the fact that the non-normal distributed errors are often recommended in various economic and financial applications.
Furthermore, our extensive simulation studies found that the KSD-based test not only shows its generality advantage to deal with multivariate non-normal null distributions, but also exhibits the comparative power with the existing tests to handle the multivariate normal null distribution. Finally, we studied a 3-dimensional financial time series by a VAR(3)--CCC--GARCH(1, 1) model, and the results of KSD-based test indicated that the error of this model follows a 3-dimensional multivariate $t_{8}$ distribution.
\setcounter{equation}{0} \setcounter{section}{0}
To facilitate our proofs, we need a useful expansion of $\widehat{\mathbb{S}}$. First, we give some notation to present this expansion. Denote $\zeta_n=\widehat\theta_n-\theta_0$ and
where
Second, we define three U-statistics $\mathbb{S}^{(a)}$ (for $a=0, 1, 2$) as follows:
where
with $b_{rs}=\big(\nabla_{\theta_r} G_{ij}(\theta_0)^\top\big)H_{ij}\big(\nabla_{\theta_s} G_{ij}(\theta_0)\big).$
With these notation, by Taylor's expansion we have
where $H_{ij}^{\dagger}=H(\eta_{ij}^\dagger)$, $\eta_{ij}^\dagger$ lies between $\eta_{ij}$ and $\widehat\eta_{ij}$, and
Furthermore, by Taylor's expansion again we have
where $\theta^\dagger$ lies between $\theta_0$ and $\widehat\theta_n$, and
By ((ref)) and ((ref))--((ref)), it follows that
where the U-statistics $\mathbb{S}^{(a)}$ (for $a=0, 1, 2$) are defined in ((ref)), and the remainder term $\widehat{R}$ is defined by
with $R_{ij}=R_{ij}^{(1)}+R_{ij}^{(2)}$ and
From the expansion ((ref)), it is clear that the estimation effect has an impact on the limiting distribution of $\widehat{\mathbb{S}}$ through the linear term $\zeta_n^\top\mathbb{S}^{(1)}$, the quadratic term $\zeta_n^\top\mathbb{S}^{(2)}\zeta_n$, and the remainder term $\widehat{R}$, and that the effect of unobserved initial values is involved in the remainder term $\widehat{R}$ via $\bar{R}_{ij}^{(2)}$.
To prove Theorems (ref)--(ref), we need two technical lemmas to handle the effects of estimation uncertainty and unobserved initial values.
Proof of Lemma (ref). By Assumptions (ref), (ref), (ref) and (ref)(ii) and the law of large numbers for U-statistics, it is not hard to see that $\mathbb{S}^{(1)}=O_{p}(1)$ and $\mathbb{S}^{(2)}=O_{p}(1)$. Since $\sqrt{n}\zeta_n=O_p(1)$ by Assumption (ref), it follows that $n_0\zeta_n^\top\mathbb{S}^{(1)}=O_{p}(n_0/\sqrt{n})=O_{p}(1/n^{\varepsilon-1/2})$ and $n_0\zeta_n^\top\mathbb{S}^{(2)}\zeta_n=O_{p}(n_0/n)=O_{p}(1/n^{\varepsilon})$. Hence, the conclusions hold. \qed
Proof of Lemma (ref). For simplicity, we only show that
since the proof for $R^{(2)}_{ij}$ is similar and even simpler. By ((ref)), we can rewrite $R^{(1)}_{ij}$ as
Then, it follows that $2(n_0-1)^{-1}\sum_{n-n_0+1\le i<j\le n}R^{(1)}_{ij}=\sum_{a=1}^6 \Delta_a^{(1)}$, where $$\Delta_a^{(1)}=\dfrac{2}{n_0-1}\sum_{n-n_0+1\le i<j\le n} r^{(1)}_{a,ij}.$$
Let $K>0$ be a generic constant whose value may change from place to place. Next, we show that $\Delta_1^{(1)}=o_{p}(1)$. To facilitate it, we claim
where $o_p(1)$ holds uniformly in $i,j$. With loss of generality, we prove ((ref)) for $\nabla_{x,x}u(\hat\eta_i^\dagger,\hat\eta_j^\dagger)-\nabla_{x,x}u(\eta_i,\eta_j)$, the first block entry of $H_{ij}^\dagger-H_{ij}$. Denote $u(x,x'):=u_1(x,x')+u_2(x,x')+u_3(x,x')+u_4(x,x')$, where $u_1(x,x')=s_{p_0}(x)^\top k(x,x')s_{p_0}(x')$, $u_2(x,x')=s_{p_0}(x)^\top k_{x'}(x,x')$, $u_3(x,x')=k_x(x,x')^\top s_{p_0}(x')$, and $u_4(x,x')=\text{trace}\big(k_{xx'}(x,x')\big)$. Below, we first prove
Rewrite
Note that
where $o_p(1)$ in ((ref)) holds uniformly in $t$ due to the fact that $\sqrt{n}\|\widehat\theta_n-\theta_0\|=O_p(1)$ and $\nonumber n^{-1/2}\max_{1\le t\le n}\sup_{\theta}\left\lVert\nabla_{\theta} g_t(\theta)\right\rVert=o_p(1) $ by Assumptions (ref)(i) and (ref), and ((ref)) holds by Assumption (ref). Therefore, by ((ref))--((ref)) and Assumption (ref), the adding and subtracting arguments give us
Similarly, the same result holds for $\|T_b^{(r,s)}(\widehat\eta_i^\dagger,\widehat\eta_j^\dagger)-T_b^{(r,s)}(\eta_i,\eta_j)\|$ with $b=2,3,4.$ Hence, the result ((ref)) holds, and then we can show the same result for $u_{b}(\cdot,\cdot)$ with $b=2,3,4.$ Therefore, it entails that the result ((ref)) holds.
Note that $E\left\lVert\eta_{ij}\right\rVert^4\le K(E\left\lVert\eta_{i}\right\rVert^4+E\left\lVert\eta_{j}\right\rVert^4)<\infty$ by Assumption (ref), $E\|\bar R_{ij}^{(2)}\|^4\le K(E\|\bar R_{i}^{(2)}\|^4+E\{\bar R_{j}^{(2)}\|^4)$, and
by Assumption (ref). Hence, by ((ref)) and H\"older's inequality, we can show
implying that $\Delta_1^{(1)}=o_{p}(1)$.
Furthermore, by Taylor's expansion, Assumptions (ref)--(ref), and a similar argument as for ((ref)), it is straightforward to see
where $o_p(1)$ holds uniformly in $i,j$. By ((ref)) and the similar arguments as for ((ref)), we can show that $\Delta_a^{(1)}=o_p(1)$ for $2\leq a\leq 6$. Therefore, it follows that the result ((ref)) holds. This completes the proof. \qed
Proof of Theorem (ref). By ((ref)) and Lemmas (ref)--(ref), $n_0\widehat{\mathbb{S}}=n_0{\mathbb{S}}^{(0)}+o_{p}(1)$, and the result follows by Theorem 4.1(2) in Liu et al. (2016). This completes the proof. \qed
Proof of Theorem (ref). By ((ref)) and Lemmas (ref)(ii) and (ref),
Now, the conclusion holds since $\sqrt{n_0}\big(\widehat{\mathbb{S}}-\mathbb{S}(p,p_0)\big)=O_{p}(1)$ by Theorem 4.1(1) in Liu et al. (2016), $\sqrt{n}\zeta_n=O_{p}(1)$, $\mathbb{S}^{(1)}=O_{p}(1)$, and $\mathbb{S}(p,p_0)>0$ under $H_1$. This completes the proof. \qed
1. Mardia's tests. Consider the null hypothesis that
Mardia (1974) detected $H_0$ in ((ref)) by proposing the following two test statistics: $$\widehat{\mathbb{T}}_{M,1} = \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n m_{ij}^3 \quad \mbox{and} \quad \widehat{\mathbb{T}}_{M,2} = \frac{1}{n} \sum_{i=1}^n m_{ii}^2,$$ where $m_{ij} = (Y_i - \bar{Y})^\top S_{Y}^{-1}(Y_j - \bar{Y})$, and $\bar{Y}$ and $S_{Y}$ are the sample mean and variance of $\{Y_t\}_{t=1}^{n}$, respectively. The tests $\widehat{\mathbb{T}}_{1,M}$ and $\widehat{\mathbb{T}}_{2,M}$ make use of the multivariate extensions of skewness and kurtosis measures in Mardia (1970), and they have the following limiting null distributions $$(n/6)\widehat{\mathbb{T}}_{M,1}\xrightarrow{d} \chi^2_{d(d+1)(d+2)/6} \mbox{ and }\widehat{\mathbb{T}}_{M,2}\xrightarrow{d} N(d(d+2), 8d(d+2)/n).$$
2. Doornik--Hansen test. Let $s = m_3/m_2^{3/2}$ and $k = m_4/m_2^2$ be the original sample skewness and kurtosis, where $m_j = \frac{1}{n}\sum_{i=1}^n(Y_i - \bar{Y})^j$. Next, transform $s$ and $k$ into $z_1$ and $z_2$, respectively, where $$z_1 = \delta \log(y+\sqrt{y^2-1})\mbox{ and }z_2 =\sqrt{9\alpha}\Big(\frac{1}{9\alpha}-1+\sqrt[3]{\frac{\chi}{2\alpha}}\Big).$$ Here,
Based on $z_1$ and $z_2$, Doornik and Hansen (2008) proposed the test statistic $\widehat{\mathbb{T}}_{DH}:= z_1^2 + z_2^2$ to detect $H_0$ in ((ref)), where the limiting null distribution of $\widehat{\mathbb{T}}_{DH}$ is $\chi^2_{2}$.
3. Henze--Zirkler test. To detect $H_0$ in ((ref)), Henze and Zirkler (1990) proposed a test statistic given by $$\widehat{\mathbb{T}}_{HZ}:= \frac{1}{n}\sum_{i=1}^n\sum_{j=1}^n e^{-\frac{\beta^2}{2}D_{ij}} - 2(1+\beta^2)^{-\frac{d}{2}}\sum_{i=1}^n e^{-\frac{\beta^2}{2(1+\beta^2)}D_i} + n(1+2\beta^2)^{-\frac{d}{2}},$$ where $\beta = \frac{1}{\sqrt{2}}\big[\frac{n(2d+1)}{4}\big]^{\frac{1}{d+4}}, $ $D_{ij}=(Y_i - Y_j)^\top S_{Y}^{-1}(Y_i - Y_j)$ is the squared Mahalanobis distance between $Y_i$ and $Y_j$, and $D_i= (Y_i - \bar{Y})^\top S_{Y}^{-1}(Y_i-\bar{Y})$ is the squared distance of $Y_i$ to the centroid.
Under $H_0$ in ((ref)), the limiting null distribution of $\widehat{\mathbb{T}}_{HZ}$ is log-normal with mean $\mu_{HZ}$ and variance $\sigma_{HZ}^2$, where
with $\omega_\beta = (1+\beta^2)(1 + 3\beta^2)$. Note that Henze and Zirkler (1990) suggested that this test is proper for sample size $n \geq 20$.
4. Bai--Chen test. For model ((ref)), Bai and Chen (2008) tested the multivariate normal and $t_{\nu}$ distributions for $\eta_t$ by using the martingale transformation. Their testing method requires the explicit formula of $P_0(Y_{it}|Y_{1t,...,Y_{i-1,t}})$ for $i=1,...,d$, where $P_0$ is the c.d.f. of $\eta_t$ under $H_0$ in ((ref)). However, it is difficult to derive the explicit formula of $P_0(Y_{it}|Y_{1t,...,Y_{i-1,t}})$ for $d>2$, even when $P_0$ is the c.d.f. of multivariate normal or $t_{\nu}$. Below, we only consider the case of $d=2$ as in Bai and Chen (2008).
Partition $$M(I_{t-1};\theta) = \left[
\right ] and C(I_{t-1};\theta) = \left[
\right ].$$ Denote $\widehat{\mu}_{it}=\mu_{i}(\widehat{I}_{t-1};\widehat{\theta}_{n})$, $\widehat{\sigma}_{it}=\sigma_{i}(\widehat{I}_{t-1};\widehat{\theta}_n)$, and $\widehat{\sigma}_{ij,t}=\sigma_{ij}(\widehat{I}_{t-1};\widehat{\theta}_n)$. Define
where
with $C_k(s) = \int_s^1 \dot{g}_{k}(r)\dot{g}_{k}^\top(r)dr$, $\dot{g}_{k}(r)$ is the first derivative of $g_{k}(r)$, and $$\widehat{J}_{n,k}(r) = \frac{1}{\sqrt{n}} \sum_{t=1}^n[I(\widehat{U}_{kt} \leq r) - r] \mbox{ for }k = 1,2,\,\, \widehat{J}_{n,3}(r)=\frac{1}{\sqrt{2}}\Big[\widehat{J}_{n,1}(r)+\widehat{J}_{n,2}(r)\Big].$$ The choices of $\widehat{U}_{kt}$ and $g_{k}(r)$ are given as follows:
Note that $\widehat{\mathbb{T}}_{BC,i}$, $i=1,2,3$, can be computed by using a similar numerical method as in Appendix B of Bai (2003). Under $H_0$ in ((ref)), the limiting distributions of $\widehat{\mathbb{T}}_{BC,i}$ can be found in Corollary 3.2 of Bai and Chen (2008). Let $cv_{BC,i}=(cv_{i,0.01},cv_{i,0.05},cv_{i,0.1})$ be a vector containing the critical values of $\widehat{\mathbb{T}}_{BC,i}$ at levels 1%, 5% and 10%. By direction simulations, we have that $cv_{BC,1}=(2.211,2.469,2.993)$, $cv_{BC,2}=(3.443,3.792,4.504)$, and $cv_{BC,3}=(2.782,2.214,1.940)$.
5. Henze--Jim\'{e}nez-Gamero--Meintanis test. When $p_0$ in ((ref)) is multivariate normal, Henze et al. (2019) made use of the identity ((ref)) to propose a test statistic given by
where $\gamma_0>0$ is a fixed constant. As the simulation studies in Henze et al. (2019), we take $\gamma_0=1.5$ and use a similar parametric bootstrap as ours in Subsection 3.3 to compute the critical values of $\widehat{\mathbb{T}}_{HJM}$.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Arellano-Valle, R. B. and Azzalini, A. (2008). The centred parametrization for the multivariate skew-normal distribution. {\it Journal of Multivariate Analysis} {\bf 99}, 1362--1382.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bai, J. (2003). Testing parametric conditional distributions of dynamic models. {\it Review of Economics and Statistics} {\bf 85}, 531--549.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bai, J. and Chen, Z. (2008). Testing multivariate distributions in GARCH models. {\it Journal of Econometrics} {\bf 143}, 19--36.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bai, J. and Ng, S. (2005). Tests for skewness, kurtosis, and normality for time series data. {\it Journal of Business & Economic Statistics} {\bf 23}, 49--60.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bauwens, L. and Laurent, S. (2005). A new class of multivariate skew densities, with application to generalized autoregressive conditional heteroscedasticity models. {\it Journal of Business & Economic Statistics} {\bf 23}, 346--354.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bauwens, L., Laurent, S. and Rombouts, J. V. K. (2006). Multivariate GARCH models: a survey. {\it Journal of Applied Econometrics} {\bf 21}, 79--109.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Berk, J. (1997). Necessary conditions for the CAPM. {\it Journal of Economic Theory} {\bf 73}, 245--257.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bontemps, C. and Meddahi, N. (2005). Testing normality: a GMM approach. {\it Journal of Econometrics} {\bf 124}, 149--186.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Bontemps, C. and Meddahi, N. (2012). Testing distributional assumptions: A GMM approach. {\it Journal of Applied Econometrics} {\bf 27}, 978--1012.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Christoffersen, P. F. and Diebold, F. X. (1997). Optimal prediction under asymmetric loss. {\it Econometric Theory} {\bf 13}, 808--817.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Comte, F. and Lieberman, O. (2003). Asymptotic theory for multivariate GARCH processes. {\it Journal of Multivariate Analysis} {\bf 84}, 61--84.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt De Luca, G., Genton, M. G. and Loperfido, N. (2006). A multivariate skew-garch model. {\it Advances in Econometrics} {\bf 20}, 33--57.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Diebold, F. X., Gunther, T. A. and Tay, A. S. (1998). Evaluating density forecasts with applications to financial risk management. {\it International Economic Review} {\bf 39}, 863--883.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Doornik, J. A. and Hansen, H. (2008). An omnibus test for univariate and multivariate normality. {\it Oxford Bulletin of Economics and Statistics} {\bf 70}, 927--939.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Escanciano, J. C. (2006). Goodness-of-fit tests for linear and non-linear time series models. {\it Journal of the American Statistical Association} {\bf 101}, 531--541.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Francq, C., Jim\'{e}nez-Gamero, M. D. and Meintanis, S. G. (2017). Tests for conditional ellipticity in multivariate GARCH models. {\it Journal of Econometrics} {\bf 196}, 305--319.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Francq, C. and Zako\"{i}an, J.-M. (2012). QML estimation of a class of multivariate asymmetric GARCH models. {\it Econometric Theory} {\bf 28}, 179--206.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Francq, C. and Zako\"{i}an, J.-M. (2019). {\it GARCH Models: Structure, Statistical Inference and Financial Applications (2nd Edition)}. Wiley, Chichester, UK.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Giacomini, R., Politis, D. N. and White, H. (2013). A warp-speed method for conducting monte carlo experiments involving bootstrap estimators. {\it Econometric Theory} {\bf 29}, 567--589.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Haas, M., Mittnik, S. and Paolella, M. S. (2004). Mixed normal conditional heteroskedasticity. {\it Journal of Financial Econometrics} {\bf 2}, 211--250.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Hafner, C. M. and Preminger, A. (2009). On asymptotic theory for multivariate GARCH models. {\it Journal of Multivariate Analysis} {\bf 100}, 2044--2054.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Henze, N., Hl\'{a}vka, Z. and Meintanis, S. G. (2014). Testing for spherical symmetry via the empirical characteristic function. {\it Statisics} {\bf 48}, 1282--1296.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Henze, N., Jim\'{e}nez-Gamero, M. D. and Meintanis, S. G. (2019). Characterizations of multinormality and corresponding tests of fit, including for GARCH models. {\it Econometric Theory} {\bf 35}, 510--546.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Henze, N. and Zirkler, B. (1990). A class of invariant consistent tests for multivariate normality. {\it Communications in Statistics--Theory and Methods} {\bf 19}, 3595--3618.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Hong, Y. and Lee, Y. J. (2005). Generalized spectral tests for conditional mean models in time series with conditional heteroskedasticity of unknown form. {\it Review of Economic Studies} {\bf 72}, 499--541.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Horv\'{a}th, L. and Zitikis, R. (2006). Testing goodness of fit based on densities of GARCH innovations. {\it Econometric Theory} {\bf 22}, 457--482
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Khmaladze, E. V. (1982). Martingale approach in the theory of goodness-of-fit tests. {\it Theory of Probability & Its Applications} {\bf 26}, 240--257.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Klar, B., Lindner, F. and Meintanis, S. G. (2012). Specification tests for the error distribution in GARCH models. {\it Computational Statistics & Data Analysis} {\bf 56}, 3587--3598.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Koul, H. and Ling, S. (2006). Fitting an error distribution in some heteroscedastic time series models. {\it Annals of Statistics} {\bf 34}, 994--1012.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Ling, S. and McAleer, M. (2003). Asymptotic theory for a new vector ARMA-GARCH model. {\it Econometric Theory} {\bf 19}, 280--310.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Liu, Q., Lee, J. and Jordan, M. (2016). A kernelized Stein discrepancy for goodness-of-fit tests. In {\it International Conference on Machine Learning} (pp. 276--284).
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt L\"{u}tkepohl, H. (2005). {\it New Introduction to Multiple Time Series Analysis}. Springer, Berlin.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Lobato, I. N. and Velasco, C. (2004). A simple test of normality for time series. {\it Econometric Theory} {\bf 20}, 671--689.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. {\it Biometrika} {\bf 57}, 519--530.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Mardia, K. V. (1974). Applications of some measures of multivariate skewness and kurtosis for testing normality and robustness studies. {\it Sankhy A} {\bf 36}, 115--128.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Mecklin, C. J. and Mundfrom, D. J. (2004). An appraisal and bibliography of tests for multivariate normality. {\it International Statistical Review} {\bf 72}, 123--138.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Sz\'{e}kely, G. J. and Rizzo, M. L. (2005). A new test for multivariate normality. {\it Journal of Multivariate Analysis} {\bf 93}, 58--80.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Taylor, J. W. (2019). Forecasting value at risk and expected shortfall using a semiparametric approach based on the asymmetric Laplace distribution. {\it Journal of Business & Economic Statistics} {\bf 37}, 121--133.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Tsay, R. S. (2005). {\it Analysis of Financial Time Series}. John Wiley & Sons, Hoboken, NJ.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Tsay, R. S. (2013). {\it Multivariate Time Series Analysis: With R and Financial Applications}. John Wiley & Sons, Hoboken, NJ.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Zhu, K. and Li, W. K. (2015). A new Pearson-type QMLE for conditionally heteroskedastic models. {\it Journal of Business & Economic Statistics} {\bf 33}, 552--565.
\hbox to -29.0pt{\hfil} \hangindent=20.0pt \hangafter=1\hskip 10pt Zhu, K. and Ling, S. (2015). Model-based pricing for financial derivatives. {\it Journal of Econometrics} {\bf 187}, 447--457.