EconBase
← Back to paper

Spectral Targeting Estimation of $λ$-GARCH models

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.

52,014 characters · 14 sections · 38 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Spectral Targeting Estimation of $-$GARCH models

\abstract{ This paper presents a novel estimator of orthogonal GARCH models, which combines (eigenvalue and -vector) targeting estimation with stepwise (univariate) estimation. We denote this the spectral targeting estimator. This two-step estimator is consistent under finite second order moments, while asymptotic normality holds under finite fourth order moments. The estimator is especially well suited for modelling larger portfolios: we compare the empirical performance of the spectral targeting estimator to that of the quasi maximum likelihood estimator for five portfolios of 25 assets. The spectral targeting estimator dominates in terms of computational complexity, being up to 57 times faster in estimation, while both estimators produce similar out-of-sample forecasts, indicating that the spectral targeting estimator is well suited for high-dimensional empirical applications. } \\ \\ Keywords: Asymptotic theory, Multivariate GARCH, Variance targeting, Two-step estimation. \\ JEL classifications: C32, C58.

Introduction

Multivariate conditionally heteroskedastic (MGARCH) models are a popular tool for risk management and dynamic portfolio allocation, where forecasts of conditional covariance matrices play an important role. As well known, MGARCH models suffer from the “curse of dimensionality”, making them difficult and time consuming to estimate for larger portfolios using quasi maximum likelihood (QML) techniques. Many practitioners and academics alike have therefore preferred using alternative estimation methods: Two popular choices are the variance targeting (VT) estimator and the equation-by-equation (EbE) estimator, see e.g. bauwens2006.

In the context of orthogonal GARCH models, such as the conditional eigenvalue GARCH ($\lambda-$GARCH) model of Hetland2019, we can combine the idea behind the two methods in what we denote the spectral targeting estimator (STE): By estimating the unconditional eigenvalues and -vectors using a sample moment estimator, the remainder of the parameters of the GARCH model may be estimated univariately in a stepwise manner, in which we target the unconditional eigenvalues and -vectors. This estimation procedure dramatically reduces the computational complexity of the optimization problem and speeds up numerical estimation compared to the QML estimator.

In this paper, we derive the large-sample properties of this two step estimator. Numerical illustrations show that the estimator is superior to the QML estimator in cross-sections larger than 10 financial assets, being up to 57 times faster in estimation, while the out-of-sample forecasts from the QML and ST estimator are similar in portfolios of 25 assets.

In general, asymptotic theory of QML estimation in MGARCH models is well-understood (see e.g. francq2019 (chapter 11) for a review of existing theory), whereas less attention has been paid to alternative estimation methods. Large sample properties of the two step VT estimator are considered in Pedersen2014 and francq2014 for the BEKK model engle1995 and (extended) CCC model (bollerslev1990 and jeantheau1998) respectively, while Francq2016 consider the two step EbE estimator for various MGARCH specifications. Both the VT and EbE estimators are two-step estimators, which are quite common in econometrics, see e.g. Newey1994. The EbEE and VTE both aim at making high(er) dimensional estimation feasible, and do so in two distinct ways: The EbEE estimates univariate volatility models in a first step, and subsequently a (conditional) correlation dynamic in a second step, whereas the VTE estimates the unconditional covariance matrix using a moment estimator, followed by a joint (profiled) estimation of the volatility and covariance dynamics. The ST estimator is related to both, as we recover sample eigenvalues and -vectors from the unconditional covariance matrix, and estimate univariate dynamics for “rotated” (orthogonalized) returns in a second step. The resulting estimator is well-behaved and easily implemented: Because the $\lambda-$GARCH model is specified using the spectral decomposition, the (profiled) log-likelihood, conditional on the initial estimator, can be rewritten as a sum of orthogonal univariate log-likelihood functions, making stepwise estimation feasible. This also means that the ST estimator is, in terms of asymptotic theory, equivalent to the VT estimator of the $\lambda-$GARCH. Furthermore, by recovering the (constant conditional) eigenvectors we avoid having to parameterize the eigenvectors under the restriction of orthonormality.

Consistency of the ST estimator follows under mild conditions (finite second order moments), while asymptotic normality requires finite fourth order moments, which may be violated empirically in financial data. Hence, the estimator is useful in the sense that it helps circumvent numerical issues often associated with (estimation of) MGARCH processes, and will produce consistent estimates under mild assumptions. However, it may not be suitable for inference on the parameters of the MGARCH process because of the moment requirement. Both of these moment conditions for consistency and asymptotic normality stem from the first step estimator of the unconditional eigenvalues and eigenvectors, which uses the sample covariance estimator. This is in contrast to the joint QML estimator of the $\lambda-$GARCH, which requires fractional moments for consistency and finite $2+\delta$ moments, for $\delta>0$, for asymptotic normality Hetland2019. \\ \\ The remainder of the paper proceeds as follows: Section (ref) introduces the $\lambda-$GARCH model and spectral targeting. Section (ref) presents the two-step estimator and Section (ref) presents novel asymptotic results and discuss practical considerations for implementation. Section (ref) investigates the empirical properties of the estimator compared to the QML estimator. Finally, Section (ref) concludes. All proofs are relegated to the appendices.

Notation

Some notation used throughout the paper. $\mathbb{R}$ denotes the real numbers, $\mathbb{N}$ the natural numbers, and $\mathbb{Z}$ the positive natural numbers. The absolute value of $a\in\mathbb{R}$ is denoted $|a|$. For $p,n\in \mathbb{Z}$, $I_p$ denotes the $(p\times p)$ identity matrix and $0_{n\times p}$ denotes a $n\times p$ matrix of zeros. The vector $\text{vec}(A)$ stacks the columns of the matrix A. We use the “diag” operator in two ways: If $W$ is a $p\times 1$ vector, $\text{diag}(W)$ returns a $p\times p$ diagonal matrix with $W$ on the diagonal, and if $A$ is a $p\times p$ matrix, $\text{diag}(A)$ returns the diagonal of $A$ as a $p\times 1$ vector. The trace of a square matrix is denoted $\text{tr}(A)$, and the determinant $\det(A)$. Furthermore, denote by $\rho(A)$ the spectral radius of any square matrix $A$, i.e. $\rho(A)=\max\{|\tilde\lambda_i|: \ \tilde\lambda_i \text{ is an eigenvalue of A}\}.$ We use $||\cdot||$ as a matrix norm. Let $\odot$ denote the Hadamard product, with $A^{\odot2}= A\odot A$, and $A\otimes B$ denotes the Kronecker product between A and B, and note that $A^{\otimes2}=A\otimes A$. Elements of matrices or vectors are denoted by lower case letters, e.g. $a_{ij}$ is the $(i,j)'th$ element of the matrix $A$. We use three kinds of convergence of random variables, $\overset{a.s.}{\rightarrow}$ denotes almost sure convergence, $\overset{p}{\rightarrow}$ denotes convergence in probability and $\overset{D}{\rightarrow}$ denotes convergence in distribution.

The $\lambda-$GARCH model

As in Hetland2019, we focus on the class of O-GARCH models originally introduced by alexander1997. The presented model has more general dynamics than the O-GARCH, allowing for eigenvalue-spillovers, and we denote this version of the model the Eigenvalue GARCH, or $\lambda-$GARCH for short.

Let $X_t$ be a $p\times1$ vector of asset returns,

align[align omitted — 48 chars of source]

where $t=1,\hdots, T$ and $Z_t$ is an $iid(0,I_p)$ sequence of random variables. $H_t^{1/2}=V\Lambda_t^{1/2}$ is the (asymmetric) matrix square root of the conditional covariance matrix, $H_t$ (following the literature on MGARCH models, see e.g. van2002 and lanne2007). , which is decomposed using the spectral theorem,

align[align omitted — 35 chars of source]

$V=

pmatrix[pmatrix omitted — 40 chars of source]

$ is an orthonormal matrix of eigenvectors, $VV'=I_p$, and $\Lambda_t$ is a diagonal matrix with time-varying eigenvalues, $\lambda_t$, on the diagonal,

align[align omitted — 48 chars of source]

The $p\times1$ vector of dynamic eigenvalues are assumed to follow a GARCH dynamic,

align[align omitted — 88 chars of source]

where $Y_t=V'X_t$ are “rotated” (or orthogonalized) returns: The orthonormal matrix $V$ rotates the returns $X_t$ to be linearly independent with conditional covariance $\Lambda_t$. To ensure that the covariance matrix is positive definite for all $t\in \mathbb{Z}$, we restrict $w_i>0$, $a_{ij}\geq0$, and $b_{ij}\geq0$ for $i,j=1,\hdots p$. Furthermore, to facilitate stepwise estimation, we restrict $B$ to be a diagonal matrix, letting the $i$'th lagged eigenvalue enter in equation $i$.

By Lemma (ref) and (ref) in Appendix (ref), the stochastic process $\{X_t\}_{t\in\mathbb{Z}}$ can be initiated from the invariant distribution such that it is covariance stationary if and only if $\rho(A+B)<1$. If this is the case, the unconditional covariance matrix, $H = V(X_t) = E[X_tX_t']$, exists almost surely and is given by,

align[align omitted — 98 chars of source]

where $\lambda=E[\lambda_t]$ is the vector of unconditional eigenvalues.

To obtain the covariance targeting (alternatively eigenvalue targeting) $\lambda-$GARCH, we re-parameterize the model by substituting (ref) into (ref),

align[align omitted — 78 chars of source]

This implies that the $i$th rotated return is driven by an augmented GARCH(1,1) with spill-overs from the other squared rotated returns,

align[align omitted — 159 chars of source]

where $w_i = (1-b_i)\lambda_i-\sum_{j=1}^pa_{ij}\lambda_j$ for $i=1,\hdots,p$, and $y_{i,t}=V_i'X_t$. This specification is motivated by generality: it seems restrictive to assume that the conditional variance of a component is not influenced by the past of other components, and allowing for spill-overs between assets may improve the model fit and out-of-sample performance.

Spectral targeting estimation

While theory for classical joint QMLE of O-GARCH type models have been considered in Hetland2019, we consider spectral targeting estimation (STE). The stepwise estimation procedure examined in this paper makes estimation and inference for the $\lambda-$GARCH feasible, even in large systems, as long as the time series dimension dominates the cross-sectional dimension (ledoit2004,ledoit2012).

Define $\upsilon = \text{vec}(V)$, i.e. the vector of stacked eigenvectors, such that

align[align omitted — 124 chars of source]

where $\gamma$ contain the eigenvalues and -vectors of the unconditional covariance matrix, $H$. Hence, $\gamma$ denote “static” and $\kappa^{(i)}$ the “dynamic” parameters of equation $i$, such that $\theta^{(i)}=[\gamma', \kappa^{(i)\prime}]'$ is the vector of parameters associated with the $i$th rotated return, $i=1,\hdots,p$, of size $e=p^2+2p+1$. Likewise, define the parameter space $\Theta^{(i)} := \mathcal{L} \times \mathcal{V} \times \mathcal{K}^{(i)} \subset \mathbb{R}^{p}_{++}\times\mathbb{R}^{p^2}\times \mathbb{R}^{p+1}$ which is restricted such that $\rho(A+B)<1$, $W$ and $\lambda$ are element-wise strictly positive and eigenvectors are orthonormal, $VV'=I_p$, such that $H$ is positive definite and symmetric. The vector of all the parameters in the model is

align*[align* omitted — 79 chars of source]

which has $p(p+1)/2+p^2+p$ elements. To emphasize the dependence on the parameters in $\theta^{(i)}$, we restate the model for the $i'$th rotated return as,

align*[align* omitted — 204 chars of source]

which also explicitly states that the conditional eigenvalues are a non-linear function of the eigenvectors in $\gamma$, and linear in the dynamic parameters in $\kappa^{(i)}$. Furthermore, $H_t(\theta)=V\Lambda_t(\theta)V'$, such that the (constant conditional) eigenvectors only depend on $\gamma$, whereas the diagonal matrix of conditional eigenvalues depend on the full vector of parameters, $\theta$.

The STE consists of two steps: In the first step, we estimate $\gamma$ using a sample estimator. In the second step, the dynamic parameters of the model are estimated by univariate QMLE for each equation in (ref)-(ref) for $i=1,\hdots,p$. This procedure yields the STE for equation $i$, denoted $\theta^{(i)}$, and based on the joint vector of parameters, $\theta$, the sequence of filtrated conditional covariance matrices, $H_t(\theta)$, can be recovered for $t=1,\hdots,T$.

The moment estimator

The first step of the STE utilizes the (strong) law of large numbers for strictly stationary and ergodic processes, and we estimate $H$ by the sample covariance matrix,

align[align omitted — 107 chars of source]

If $X_t$ is covariance stationary and ergodic, $\hat H$ is a strongly consistent estimator for $H$ by the ergodic theorem. From $\hat H$ it is possible to recover the estimated eigenvalues, $\hat\lambda$, and estimated eigenvectors, $\hat V$, by solving the two equations,

align[align omitted — 107 chars of source]

and under Assumption (ref) below $\hat \lambda$ and $\hat \upsilon$ are strongly consistent estimators of $\lambda$ and $\upsilon$ respectively by the continuous mapping theorem.

In applications, these two equations are solved using iterative procedures and for $H$ symmetric and positive definite, all eigenvalues are almost surely strictly positive. Notice however, that the eigenvalue decomposition is not unique: the spectrum of $H$ is unique only up to the ordering, and while the eigenspace of $H$ is unique the eigenvectors are not. Furthermore, eigenvalues may not be unique. We discuss this further in remark (ref).

remark[Alternative first step estimator] Instead of estimating the eigenvalues and -vectors implicitly using the moment estimator of $H$, we can estimate them directly using an approach similar to that proposed by fan2008 and boswijk2011, wherein $V$ is specified using rotation matrices, \begin{align*} V(\phi)=\prod_{1\leq i< j \leq p}U_{ij}(\phi_{ij}), \end{align*} with $U_{ij}(\phi_{ij})$ a $p$-dimensional identity matrix apart from four elements: $(i,i)$ and $(j,j)$ are $\cos(\phi_{ij})$, $(i,j)$ and $(j,i)$ are $\sin(\phi_{ij})$ and $-\sin(\phi_{ij})$ respectively. $\phi$ is a $p(p-1)/2$ vector containing the rotation parameters, $\phi_{ij}$. This parameterization ensures that $V(\phi)V'(\phi)=I_p$. The eigenvectors and eigenvalues can then be estimated by numerically solving the minimization problem, \begin{align*} \underset{[\phi',\lambda']'\in\mathcal{C}}{\arg\min} \ \ C_T(\phi, \lambda) \end{align*} where $\mathcal{C}$ is an appropriate parameter space and $C_T(\phi,\lambda)$ is a cost function, e.g. the Gaussian log-likelihood, \begin{align*} C_T(\phi,\lambda) = \frac{1}{T} \sum_{t=1}^T \left (\log\det(\Lambda)+X_t'V(\phi)\Lambda^{-1}V'(\phi)X_t\right). \end{align*} The asymptotic theory for this estimator can be derived with relative ease, see e.g. Hetland2019 who parameterize the joint QMLE of the $\lambda-$GARCH in a similar fashion. One should however, keep in mind that the rotation parameters in $\phi$ are not uniquely identified unless we impose restrictions on the parameter space. A sufficient condition is $\phi_{ij}\in(0,\pi/2)$.

The alternative first step estimator outlined in remark (ref) requires numerical optimization of a cost function, and may therefore run into numerical problems as $p$ increase, such as failure of a Newton-type optimization procedure to converge, or the possibility of ending up in a local maximum -- problems similar to those of the joint QML estimator. We therefore choose to work with the sample moment estimator as it has a closed form solution and is the preferred first step estimator in the variance targeting literature.

The profiled maximum likelihood estimator

In the second step of the STE, we consider the profiled quasi log-likelihood function based on the multivariate Gaussian distribution. The joint Gaussian log-likelihood of the model, conditional on a fixed $X_0$ and $H_0$, is,

align[align omitted — 361 chars of source]

using $H_t^{(-1)}(\theta)=V\Lambda_t^{-1}(\theta)V'$, $\log\det(H_t(\theta))=\sum_{i=1}^p\log(\lambda_{i,t}(\gamma,\kappa^{(i)}))$. and $Y_t(\gamma)=V'X_t$. That is, because the rotated returns are orthogonal, the log-likelihood function can be decomposed as the sum of $p$ univariate log-likelihood functions, each of which depend on $\theta^{(i)}=[\gamma', \kappa^{(i)'}]'$,

align[align omitted — 246 chars of source]

where $y_{i,t}(\gamma)= V_i'X_t$ and $\lambda_{i,t}(\gamma,\kappa^{(i)})$ is given in (ref). Conditional on $\gamma$, each of the $i$ univariate log-likelihood functions are orthogonal and do not depend on $\kappa^{(j)}$ for $j\neq i$. The parameters of the model can therefore be estimated sequentially, and we define the STE of $\kappa^{(i)}$ as,

align[align omitted — 143 chars of source]

and the two-step procedure yields the STE of $\theta$,

align*[align* omitted — 99 chars of source]

Similar to (quasi) maximum likelihood estimation of multivariate GARCH models, we use the Gaussian log-likelihood function, but we do not assume that the vector of innovations $Z_t$ are Gaussian, only that they are centered with unit variance: Even if the innovations are drawn from a different distribution, the results in Theorem (ref) and (ref) below still hold, as long as the assumptions are satisfied.

Compared to joint QMLE, which estimates all $\frac{3}{2}(p^2+p)$ parameters jointly, the STE procedure vastly reduces the number of parameters estimated in each step: In the first step $p(p+1)/2$ parameters are estimated by method of moments and in the second step $p+1$ parameters are estimated for each rotated return, making the estimation procedure suitable in high-dimensional systems and less vulnerable to numerical problems.

Large-sample properties of sequential variance targeting estimation

In this section we establish consistency and asymptotic normality of the STE and discuss practical considerations for implementation. A novelty of the asymptotic theory presented here is that we parameterize the moment estimator in terms of the unconditional eigenvalues and vectors, rather than the vectorized covariance matrix. In doing so, we apply the mean-value theorem on the eigenvectors, which otherwise do not have a closed form solution as a function of the unconditional covariance matrix. This, in conjunction with the continuous mapping theorem, allows us to study the asymptotic behavior of both the first step estimator, $\gamma$, and the joint parameter vector of the $i$'th rotated return, $\theta^{(i)}$.

The two-step estimator is consistent under finite second order moments, and it has a limiting Gaussian distribution under the assumption of finite fourth order moments. Both of these moment conditions stem from the first step moment estimator, and are more strict that the moment conditions for the joint QML estimator (for which we need $E||X_t||^{2+\delta}<\infty$, $\delta>0$, see Theorem 3.3 in Hetland2019). These results are novel and extend the existing literature on targeting and stepwise estimation, see e.g. francq2014, Pedersen2014 and Francq2016. All proofs are relegated to Appendix (ref).

Before discussing the asymptotic properties in detail, we make the following assumptions. First, we assume that the process is covariance stationary and ergodic.

assumptionThe process $\{X_t\}_{t\in \mathbb{Z}}$ is strictly stationary, ergodic and has finite second order moments.

Furthermore, we need the following assumption on the algebraic multiplicity of the eigenvalues.

assumptionThe characteristic polynomial of the unconditional covariance matrix, $H$, has an algebraic multiplicity of 1.

We also assume that the dynamic parameters of the model are identified and that the true parameter vector is a subset of the parameter space.

assumptionThe true parameter vector $\theta_0^{(i)}\in\Theta^{(i)}$, with $\Theta^{(i)}$ compact.
assumptionFor $\kappa^{(i)}\in \mathcal{K}^{(i)}$, if $\kappa^{(i)}\neq\kappa_0^{(i)}$, then $\lambda_{i,t}(\gamma_0,\kappa^{(i)})\neq \lambda_{i,t}(\gamma_0,\kappa_0^{(i)})$.

These assumptions lead us to the following theorem on strong consistency of the ST estimator.

theoremUnder Assumptions (ref)-(ref), as $T\rightarrow\infty$, the ST estimator is consistent, \begin{align*} \hat\theta^{(i)}\overset{a.s}{\rightarrow}\theta_0^{(i)}. \end{align*}

Assumption (ref) is in line with the literature for variance-targeting estimation, both in the univariate and multivariate case, see e.g. Pedersen2014 or francq2011,francq2014, and is needed to ensure that the moment estimator converge to a well-defined unconditional covariance matrix for $T\rightarrow \infty$.

Assumption (ref) is novel in the (variance) targeting literature and is needed for the first step estimator: We assume that all the unconditional eigenvalues are simple, i.e. that the characteristic polynomial of the unconditional covariance matrix has an algebraic multiplicity of one. This is needed for two reasons: First, in the case of repeated eigenvalues, the associated eigenvectors are not uniquely determined and the parameters in the first step estimator are not uniquely identified, and hence the first step estimator is not consistent. Second, it is a requirement for $\lambda$ and $\upsilon$ to be continuously differentiable (Theorem 1, magnus1985), which is needed to apply the mean-value theorem when considering the asymptotic distribution of the estimator.

Assumptions (ref)-(ref) are standard for multivariate GARCH models, see e.g. COMTE200361 or Hafner2009a. Moreover, the normalization imposed on the first step estimator ensures that the eigenvalues and -vectors of the first step estimator are uniquely identified. Primitive conditions for the identification of the second step estimator can be found in e.g. francq2019 (chapter 10) and are also treated in Hetland2019.

remark[Identification of the first step estimator] The eigenvalue decomposition used in the first step estimator is not uniquely defined: the ordering of the eigenvalues is not fixed and the sign of the eigenvectors is unidentified. We can, however, without a loss of generality, sort the eigenvalues in non-decreasing order and normalize the eigenvectors such that the first non-zero element of each eigenvector is positive. These two normalizations, along with Assumption (ref), ensure that the eigenvalue decomposition is unique. Note, however, that any equivalent normalizations also suffice.

Next, we show that the estimator is asymptotically normal. To do so, we need two additional assumptions on existence of moments and the true parameter vector.

assumptionThe process $\{X_t\}_{t\in \mathbb{Z}}$ has finite fourth order moments, $E||X_t||^4<\infty$.
assumption$\theta_0^{(i)}$ is in the interior of $\Theta^{(i)}$.

This leads us to the next theorem on asymptotic normality of the estimator for the $i$'th rotated return,

theoremUnder Assumptions (ref)-(ref), for $T\rightarrow\infty$ \begin{align*} \sqrt{T}\left(\hat\theta^{(i)}-\theta_0^{(i)}\right)\overset{D}{\rightarrow}N\left(0,\Sigma_0^{(i)}\right), \end{align*} where $\Sigma^{(i)}$ is the asymptotic covariance matrix, given by, \begin{align} \underset{e\times e}{\Sigma_0^{(i)}} & = \begin{pmatrix} I_{p(p+1)} & 0_{p(p+1)\times p+1} \\ -(J_0^{(i)})^{-1}K_0^{(i)} & -(J_0^{(i)})^{-1} \end{pmatrix} \Omega_0^{(i)} \begin{pmatrix} I_{p(p+1)} & -(J_0^{(i)})^{-1}(K_0^{(i)})' \\ 0_{p+1\times p(p+1)} & -(J_0^{(i)})^{-1} \end{pmatrix}. \end{align} where ${J_0^{(i)}}$ and ${K_0^{(i)}}$ are defined in (ref) and ${\Omega_0^{(i)}}$ is given in (ref).

Assumption (ref) is required to ensure that the first step estimator, $\sqrt{T}(\hat\gamma-\gamma_0)$, converges to a Gaussian distribution with a finite variance. This assumption is common in the variance targeting literature and is also needed when reparameterizing the moment estimator in terms of the spectral decomposition. In fact, the moment requirement is not needed in the probability analysis of the profiled log-likelihood function, but it simplifies the exposition. Lemma (ref) in Appendix (ref) can be used to check the moment condition in Assumption (ref). Based on the simulations included in Appendix (ref), the moment conditions for consistency and asymptotic normality are sufficient and necessary. Assumption (ref) is standard in the literature, and is a technical requirement to ensure that the mean-value theorem can be applied on the optimality condition for the profiled log-likelihood functions.

In the derivation of the asymptotic distribution of the first step (and consequently the second step) estimator, we restate the moment estimator as the average of the conditional eigenvalues. However, as $\text{vec}(\hat H)=V^{\otimes2}\text{vec}(\frac{1}{T}\sum_{t=1}^T\Lambda_{0,t}^{1/2}Z_tZ_t'\Lambda_{0,t}^{1/2})$ is in terms of $\Lambda_t$, and not the vectorized eigenvalues, $\lambda_t$, we restate the dynamics of the conditional eigenvalues in (ref) as a (restricted) BEKK$(p^2,1,1,1)$ model for $\Lambda_t$. This parametrization is present in $\Omega_{0}^{(i)}$ in (ref). In doing so, $\hat\gamma-\gamma_0$ is a martingale difference, allowing us to use a central limit theorem on $\sqrt{T}

pmatrix[pmatrix omitted — 102 chars of source]

'$ jointly to show normality and find the expression for $\Omega^{(i)}_0$. The proof of joint normality of $\sqrt{T}(\hat\theta^{(i)}-\theta^{(i)}_0)$ applies the mean-value theorem on the optimality condition of the second step estimator, stacked with the moment estimator from step one, and lemmata (ref)-(ref) in Appendix (ref) verify that the mean-value theorem can be applied.

remark[Fixed initial values] Assumption (ref) assumes that the process is strictly stationary, implying that the process $\{X_t\}_{t\in\mathbb{Z}}$ is initiated in the invariant distribution or in the infinite past. In practice the observed process is initiated in some fixed values, $X_0$ and $H_0$, which by definition makes the process non-stationary. However, the appendix verifies that the choice of initial values are asymptotically irrelevant for both consistency and asymptotic normality of the estimator, see lemmata (ref) and (ref).

The model as presented in (ref)-(ref) restricts the $B$ matrix to facilitate stepwise estimation. However, in many applications, practitioners prefer the “diagonal” specification, in which both $A$ and $B$ are diagonal, this restriction is also feasible in terms of estimation, as discussed in the following remark.

remark[Diagonal model] The estimator is still consistent and asymptotically normal in the case of a diagonal $A$ matrix, such that the vector of dynamic parameters is $\kappa^{(i)}=[a_i,b_i]'$. Hence, by Theorem (ref), Assumption (ref) and Lemma (ref) implies that for $A$ and $B$ diagonal matrices, a sufficient and necessary condition for finite fourth order moments is, $\max(b_i^2+E[z_{i,t}^4]a_{i}+2a_ib_i)<1$ for $i=1,\hdots,p$.

In applications, the asymptotic variance matrix for the stepwise estimator may be approximated using plug-in sample estimators. That is,

align[align omitted — 1,056 chars of source]

where $(\hat\lambda_i I_p-\hat H)^+$ denotes the Penrose-Moore pseudo-inverse of $\hat\lambda_i I_p-\hat H$ for $i=1,\hdots,p$ and $D$ is defined in Lemma (ref). The expressions in (ref)-(ref) converge almost surely to their population counterparts due to the ergodic theorem. Note that we may substitute the estimators by their true value, as we have established strong consistency of $\theta^{(i)}$. This makes estimation of the asymptotic covariance matrix only slightly more cumbersome than that of the well-known “sandwich” covariance matrix estimator known from joint QMLE.

Once we know the asymptotic distribution of $\hat\theta^{(i)}$, it is possible to derive the asymptotic distribution of the intercept term, $w_i$, in the original model in (ref) using the delta method.

corollary[Limiting distribution of the intercept] Under Assumption (ref)-(ref), for $T\rightarrow \infty$, \begin{align*} & \sqrt{T}\left(\hat w_i-w_{0,i}\right)\overset{D}{\rightarrow}N(0,\Phi_0\Sigma_0^{(i)}\Phi_0'), \end{align*} for $i=1,\hdots,p$ with $\Sigma_0^{(i)}$ given in (ref) and, \begin{align} \underset{1\times e}{\Phi_0} = \begin{pmatrix} 1-b_{0,i}\mathbb{I}\{\lambda_{0,i}\}-\sum_{j=1}^p a_{0,ij}, & 0_{p^2\times 1}, & -\lambda_0, & -\lambda_{0,i} \end{pmatrix}, \end{align} where $\mathbb{I}\{\lambda_i\}$ is a $p\times 1$ vector of zeros, apart from a $1$ in the $i$'th row.

Note that we present the asymptotic theory in terms of $\theta^{(i)}$ rather than, $\theta$, following Francq2016 and their notation for an equation-by-equation estimator for various MGARCH models. The theorems listed above could easily be restated in terms of $\theta$, as the asymptotic results hold simultaneously due to the orthogonality of the conditional univariate log-likelihood functions, but we refrain from doing so for two reasons: First, the present formulation is coherent with the step-wise approach of the estimator. Second, the present formulation makes it straightforward to parallelize estimation, computing the asymptotic variance matrix in each iteration, which also speeding up the estimation procedure.

As already emphasized, the STE reduces the risk of numerical issues in estimation compared to the QML estimator, and in the context of the $\lambda-$GARCH model, the ST estimator is closely related to the variance targeting estimator: Because the profiled log-likelihood consists of $p$ orthogonal terms, the STE and VTE of the $\lambda-$GARCH are theoretically equivalent, and in practice is expected to produce similar estimates and standard errors. Note however, that we expect the STE to have a smaller computational burden, as it minimizes the log-likelihood function over a smaller parameter space.

remark[Variance targeting estimator of the $\lambda-$GARCH] An alternative to the stepwise estimation of the $\lambda-$GARCH is the variance targeting (VT) estimator, in which all $\kappa=[\kappa^{(1)\prime},\hdots,\kappa^{(p)\prime}]'$ are estimated jointly. This estimator still relies on the first step estimator of $\hat\gamma$ in (ref), and we denote the full VT estimator $\tilde\theta=[\hat\gamma',\tilde\kappa']'$. Because of the orthogonal structure of the log-likelihood function, consistency and asymptotic normality of the VT estimator of the $\lambda-$GARCH can be derived using similar techniques as in Appendix (ref) and (ref).

Empirical illustrations

In the following we compare the empirical performance of the ST estimator to that of the joint QML estimator. First, we consider the relative efficiency of the two estimators in a simulation setting for many different portfolio sizes. This exercise lets us compare (empirical) efficiency of the STE against the QMLE. Second, we consider the out-of-sample performance of the the two estimation methods in a recursive value-at-risk application for portfolios of $p=25$ assets. The empirical fit is assessed using the likelihood ratio tests of christoffersen1998. In both these exercises, we also consider the computational complexity (i.e. time spent on estimating the model) of the two methods. Finally, we briefly summarize the results.

Relative efficiency: STE vs. QMLE

We now compare the relative efficiency and the time complexity of the STE against the joint QMLE. This is done for the diagonal model of dimension $p$, where we simulate a data-generating process $N=399$ times with $A_0=0.05I_p$, $B_0=0.85I_p$, such that the process has finite fourth order moments. The unconditional eigenvalues are specified as $\lambda_{0,i}=(p+1-i)/10$ for $i=1,\hdots,p$ and the eigenvectors are constructed using rotation matrices with all rotation parameters $\phi_{0,i}=0.5$ for $i=1,\hdots,p(p-1)/2$ (see Remark (ref)). The innovations, $Z_t$, are drawn $iid$ from a standard normal distribution, and each path of the simulated process has $T=2000$ observations. The model has $p(p-1)/2+3p$ parameters, and the STE procedure estimates $p(p+1)/2$ parameters in the first step, and the remaining $2p$ parameters sequentially for each rotated return. The QMLE on the other hand estimates all $p(p-1)/2+3p$ parameters simultaneously.

In comparing the two estimators, we employ the same methodology as Francq2016 who use the quadratic form $T(\hat\vartheta-\vartheta_0)'\mathcal{I}(\hat\vartheta-\vartheta_0)$ as a measure of accuracy of an estimator $\hat\vartheta$, where $\mathcal{I}$ is the (numerically) approximated information matrix and the parameter vector is constructed identically for both estimators with $\vartheta = [\text{vec}(H)',\text{diag}(A)',\text{diag}(B)']'$. Because $\mathcal{I}$ is computationally demanding to compute in higher dimensions, we instead use the simulated information matrix, which is obtained as $\mathcal{I}=var(\hat\vartheta^n_{QMLE}-\vartheta_0)\approx \frac{1}{N}\sum_{n=1}^N(\hat\vartheta^n_{QMLE}-\vartheta_0)(\hat\vartheta^n_{QMLE}-\vartheta_0)'$ for $n=1,\hdots,N$, where $\hat\vartheta^n_{QMLE}$ is the QMLE parameter vector for the $n$'th simulated path. The relative efficiency is then computed as,

align*[align* omitted — 188 chars of source]

where $\bar\vartheta_{STE}=\frac{1}{N}\sum_{n=1}^N\hat\vartheta^n_{STE}$ with $\bar\vartheta_{QMLE}$ defined analogously. By this definition, if $RE<1$, the ST estimator is relatively more efficient than the QML estimator.

table[table omitted — 1,193 chars of source]

The (average) computation times and the relative efficiency for the two estimators are contained in table (ref). For the larger systems, $p>20$, the computation time for QMLE is very big, on average 50 minutes for $p=50$ and 260 minutes for $p=100$, whereas STE remains feasible in all but the $p=500$ case, in which the computation time is roughly one hour. Considering the relative efficiency of the two estimators, the QMLE performs favourably for $p\leq10$, after which its performance deteriorate drastically compared to the STE. For portfolios larger than 10 assets, the STE is preferred.

Here, estimations are initiated in $\vartheta_{init}=\vartheta_0-0.025$. However, one could argue that initiating the both estimation procedures in $\vartheta_{init}\neq\vartheta_{0}$ gives the joint QMLE a disadvantage, as it performs numerical optimization over a much larger parameter space. We therefore repeat the exercise, initiating in $\vartheta_{init}=\vartheta_0$. This yields almost identical results (available upon request) and leads to the same conclusion, namely that the STE is relatively more efficient than joint QMLE for systems larger than $p>10$ assets, and that it always has a lower computational complexity than joint QMLE.

An application in risk management

We now turn our attention to the empirical performance of the STE of the $\lambda$-GARCH, and compare it to the joint QML estimator. The out-of-sample performance is assessed by considering the conditional $5\%$ Value-at-Risk (VaR) for five different medium-sized portfolios consisting of $p=25$ assets from the SP100 index.

Methodology and data

We consider the out-of-sample performance by considering the conditional $5\%$ value-at-risk at 1 and 5-day horizons for five different portfolios. The first of the five portfolios is equally weighted while the weights of the remaining portfolios are drawn randomly such that the second and third portfolios are long-only, with the third portfolio $50\%$ geared. The fourth and fifth portfolios are long-short portfolios. The constituents of the portfolios are drawn randomly from the SP100 index and can, along with their weighting, be found in Appendix (ref).

Each of the three estimators is fitted on a (rolling window) sample of $T=1200$ daily observations, with the initial sample starting on December $28$th 2010 and ending on December $29$th 2015. The out-of-sample consists of 3 years of data from December $30$th 2015 to December $31$st 2018, leading to $\tau=756$ out-of-sample observations for the $1$-day forecast and $\tau=189$ observations for the $5$-day (non-overlapping) forecasts. The out-of-sample forecasts are computed using a filtered historical simulation in which we draw innovations $iid$ with replacement from the standardized residuals, $\hat Z_t$, see e.g. christoffersen2009.

Recall that the conditional VaR at risk level $\alpha$ for the $h$-period return of portfolio $i$, denoted VaR$^i_{t,h}(\alpha)$ is defined as,

align[align omitted — 68 chars of source]

where $P_t$ is the conditional distribution of the ex ante $h$-period return of portfolio $i$, $R_{t+h|h}^i$ . Define the (unconditional) “hit” variable for portfolio $i$ as,

align[align omitted — 79 chars of source]

such that the unconditional coverage for portfolio $i$ is $\pi_i = \frac{1}{\tau}\sum_{n=1}^\tau I^i_{n}$. Similarly, we define the conditional hit variable as $I^i_{t|t-1} = \textbf{1}\{I_{t}^i=1 | I_{t-1}^i=1\}$, denoting two hits in a row.

When assessing the adequacy of the VaR forecasts we consider the three likelihood ratio (LR) tests proposed by christoffersen1998. The first LR test examines the hypothesis that the unconditional coverage is correct, $E[I_t^i]=\alpha$, but fails to account for potential clustering in the VaR hits. This is rectified by the second test, in which $I_{t|t-1}^i$ follows a two-state Markov chain, and we test the hypothesis of independence between hits. However, this test does not test for correct coverage, and as a consequence, we also consider the third test of correct conditional coverage, which lets $I_{t|t-1}^i$ follow the two-state Markov chain, and tests it against the null of independence between hits and correct coverage. The tests are denoted $LR_{uc}$, $LR_{ind}$ and $LR_{cc}$ respectively.

table[table omitted — 2,496 chars of source]

Out-of-sample results

The results of the out-of-sample exercise is given in table (ref). Importantly, the STE procedure is roughly $57$ times faster than the QMLE. We note that the estimated $\lambda-$GARCH (on average) has finite second order moments but not fourth order moments. Intuitively, this means that both estimators are consistent, but only the QML estimator has a limiting Gaussian distribution.

The two estimation methods have a similar performance based unconditional coverage and the LR-tests: In general, the unconditional coverage is slightly different from the hypothesized $5\%$ and most of the LR-tests do not reject. Similar results are found for the 1 and 5-day $1\%$ VaR (not reported here). The rejected LR-tests relate to the equally weighted portfolio $P_1$.

In general, the VaR estimates produced by the two estimation methods are similar, but not identical: Consider figure (ref) which plots the estimated VaR for the two estimation methods along with the realized return of portfolio $P_4$. As shown, the VaR estimates are, for the majority of the sample, very similar, but the QML estimator sometimes produce more extreme VaR estimates than the STE. We note, however, that while the two VaR estimates at times differ, the unconditional and conditional hit sequences are almost identical, and based on the LR-test in table (ref), none of the estimation methods seem to dominate the other empirically. We therefore conclude that the estimation procedures seem to yield similar results, with the STE having the clear advantage that it is much faster in practice.

figure[figure omitted — 263 chars of source]

Brief summary of numerical exercises

The simulation evidence in (ref) indicates that not only is the STE relatively more efficient than QMLE in cross-sections of more than $p>10$ assets, it is also much more time efficient. This is verified in by the empirical study in Section (ref). One potential explanation is that the $\lambda-$GARCH is a non-linear function of the parameters in $\gamma$ through $Y_{t-1}$. By using a stepwise estimator, in which $\gamma$ is estimated using a closed form estimator, we mitigate the potential issues due to non-linearity, which seem to cause issues for large $p$ in the QML estimator.

In regards to the asymptotic results in Theorem (ref)-(ref), the simulation study in Appendix (ref) suggests that estimator is consistent in the case of finite second order moments of $X_t$. Furthermore, the simulations indicate that the asymptotic normality of the STE holds when $X_t$ has finite fourth moments. Hence, the moment requirements in Assumption (ref) and (ref) appear to be sufficient and necessary for consistency and asymptotic normality of the estimator.

Extensions and Concluding remarks

We have derived asymptotic properties of the spectral targeting estimator (STE) for the $\lambda-$GARCH, an extended version of the multivariate orthogonal GARCH (O-GARCH). The two-step estimator is consistent under finite second order moments, while it has a limiting Gaussian distribution when fourth order moments are finite. Simulations indicate that these moment conditions are sufficient and necessary. Moreover, we compare the empirical performance of the STE to that of the quasi maximum likelihood estimator (QMLE) for five portfolios of 25 assets. The STE dominates QMLE in terms of computational complexity, being up to 57 times faster in estimation, while both estimators produce similar out-of-sample forecasts. Finally, simulations indicate that in portfolios of more than 10 assets, the stepwise estimator is relatively more efficient than QMLE. The STE is therefore well suited for practitioners as it alleviates numerical problems and speeds up numerical optimization, while being easy to implement.

We note that while the STE delivered promising results in this exposition, the first step (sample) estimator may not be well-behaved when the ratio $p/T$ approaches one. This is discussed in e.g. ledoit2004,ledoit2012, who derive shrinkage estimators for the sample covariance matrix, minimizing the estimation error. An extension could therefore consider the asymptotic analysis of a spectral targeting estimator where the first step estimator is based on shrinkage. Another extension would be to consider spectral targeting estimation with infinite fourth order moments, in a similar fashion to the exposition in Pedersen2014b who consider the variance targeting estimator.

\printbibliography