Discussion Papers In Economics
And Business
A Penalised OLS Framework for High-Dimensional
Multivariate Stochastic Volatility Models

Benjamin Poignard
Manabu Asai

Discussion Paper 20-02

January 2020

Graduate School of Economics
Osaka University, Toyonaka, Osaka 560-0043, JAPANA Penalised OLS Framework for High-Dimensional
Multivariate Stochastic Volatility Models∗
Benjamin Poignard†
Graduate School of Economics, Osaka University and Riken-AIP, Japan
Manabu Asai‡
Faculty of Economics, Soka University, Japan

Abstract
Although multivariate stochastic volatility (MSV) models usually produce more accurate
forecasts compared to multivariate GARCH models, their estimation techniques such as Monte
Carlo likelihood or Bayesian Markov Chain Monte Carlo are computationally demanding andthus suffer from the so-called “curse of dimensionality”: using such methods, the applications
are typically restricted to low-dimensional vectors. In this paper, we propose a fast estimation
approach for MSV models based on a penalised ordinary least squares framework. Specifying
the MSV model as a multivariate state-space model, we propose a two-step penalised procedure
for estimating the latter using a broad range of potentially non-convex penalty functions. In
the first step, we approximate an EGARCH type dynamic using a penalised AR process witha sufficiently large number of lags, providing a sparse estimator. Conditionally on this firststep estimator, we estimate the state vector based on a AR type dynamic. This two-stepprocedure relies on OLS based loss functions and thus easily accommodates high-dimensional
vectors. We provide the large sample properties of the two-step estimator together with the socalled support recovery of the first step estimator. The empirical performances of our method
are illustrated through in-sample simulations and out-of-sample variance-covariance matrix
forecasts, where we consider as competitors commonly used MGARCH models.

Keywords: Forecasting; Multivariate Stochastic Volatility; Oracle Property; Penalised M-estimation.
JEL Classification: C13, C32.
∗

The authors are most grateful to Yoshi Baba for very helpful comments and suggestions. The second author
acknowledges the financial support of the Japan Ministry of Education, Culture, Sports, Science and Technology,
Japan Society for the Promotion of Science, and the Australian Academy of Science. The second author acknowledges the financial support of the Japan Society for the Promotion of Science.
†
E-mail address: bpoignard@econ.osaka-u.ac.jp‡
E-mail address: m-asai@soka.ac.jp1

1

Introduction

Over these past decades, various covariance models have been being developed for describing
dynamic structures for multivariate economic and financial time series. Within the Multivariate
GARCH (MGARCH) family, the dynamic conditional correlation (DCC) model of Engle (2002)
and Tse and Tsui (2002), the BEKK model of Baba et al. (1985) and Engle and Kroner (1995),
and their variants are commonly used: see the survey of Bauwens, Laurent, and Rombouts (2006)
for instance. In the multivariate stochastic volatility (MSV) family, the MSV model of Harvey,
Ruiz, and Shephard (1994) was extended, among others, by the factor model of Chib, Nardari,and Shephard (2006) and the dynamic correlation model of Asai and McAleer (2009): see, e.g.,Ghysels, Harvey, and Renault (1996), Asai, McAleer, and Yu (2006), and Chib, Omori, and Asai
(2009) for various univariate and multivariate stochastic volatility models. Based on a thoroughempirical analysis, Chib, Nardari, and Shephard (2006) highlighted that the MSV models usuallyoutperform MGARCH based models, such as, e.g, the DCC and the BEKK models in terms ofout-of-sample forecasts.
For estimating MSV models, Harvey, Ruiz, and Shephard (1994) derived a state space form
based on the vector of the log of squared returns. Using the corresponding state space form, they
performed a Kalman based filtering technique to evaluate the quasi log likelihood function. In
the recent literature, a commonly used method is the Bayesian MCMC one, as described, e.g., in
Chib, Omori, and Asai (2009) and Kastner, Frühwirth-Schnatter, and Lopes (2017), among others.
An alternative estimation approach is the Monte Carlo Likelihood (MCL) method suggested by
Durbin and Koopman (1997, 2001) and applied by Asai, Caprion, and McAleer (2015) and Asaiand McAleer (2009). Due to the severe cost in terms of computations, empirical applicationsin the literature are typically limited to low-dimensional random vectors when the MCMC and2

MCL approaches are carried out. MGARCH specifications also suffer from the so-called “curseof dimensionality” since the complexity is in general of order O(p2 ), where p corresponds to the
problem dimension as the specification of a general multivariate dynamic model often induces an
explosion of the number of free parameters, inducing practical problems of inference and possibly
overfitting. Moreover, tricky conditions are required on the model parameters to satisfy the
positive-definiteness of the variance-covariance process.
Another key hurdle of the aforementioned methods is the high non-linearity of the models,
which requires the use of likelihood based estimation approaches. Therefore, strongly reduced
versions of such multivariate models are most often considered as soon as p is larger than fouror five, typically. Another approach is given by factor modelling, which aims at reducing themodel complexity. Among others, Fan et al. (2008) emphasised the relevance of factor models forhigh-dimensional precision matrix estimation. However, this approach requires the identificationof the corresponding factors. An “expert” approach is based on some priors regarding the leading
underlying factors. Otherwise, latent unobserved factors induce particular estimation issues andtheir number is questionable.
The objective of this paper consists in modelling high-dimensional variance-covariance matrices
within the multivariate stochastic volatility framework in a flexible manner and breaking thecurse of dimensionality without relying on standard MCMC or MCL based procedures. To doso, we introduce a vector autoregressive and moving-average (VARMA) representation for the
MSV model of Harvey, Ruiz, and Shephard (1994) and apply an OLS-based two-step estimation
approach extending the idea of Hannan and Rissanen (1982) and Hannan and Kavalieris (1984).
Applying their approach to the MSV model requires an OLS estimation of a large dimensional
VAR model with a sufficiently large number of lags in the first step. Nonetheless, for the purpose3

of parsimony and to avoid overfitting, we have to enforce the nullity of possibly numerous model
coefficients. The OLS objective function is particularly adapted for regularisation/penalisation
procedures and fast closed form algorithms can be applied. Our study shares a similar spirit
with Poignard and Fermanian (2019), who provided a framework for high-dimensional variancecovariance within the MGARCH family: they derived some parameterizations to directly generate
positive-definite covariance matrices based on Multivariate ARCH processes, which enables a linear
representation with respect to the parameters and thus the use of a penalised OLS criterion at theestimation step. But our work differs from theirs in two main respects: our analysis lies within
the multivariate stochastic volatility family; we consider a general penalisation framework, which
includes a broad range of different penalty functions.
The main contributions of our method are as follows: using a parsimonious OLS framework, we can directly generate positive-definite variance-covariance matrices without relying on
MCMC/MCL like methods and manage high-dimensional matrix processes; the large sample
properties of the two-step estimator are provided together with the conditions for satisfying theso-called oracle property for the first step estimator in the sense of Fan and Li (2001), which
ensures the correct identification of the underlying set of nonzero coefficients. In the first step, we
consider a general penalised M-estimation framework, which encompasses potentially non-convex
penalty functions.
The remainder of the paper is organised as follows. In Section 2, we describe the frameworkand the new forecasting procedure based on a regularised OLS estimation framework. Section 3
contains the large sample properties of the regularised two-step OLS estimator. Section 4 reports
simulation based experiment results for in-sample estimates of covariance matrix together with
out-of-sample forecasting results based on a real financial portfolio. Finally, Section 5 concludes

4the paper. All proofs and intermediary results are in the Appendix.
Notations. Throughout this paper, we denote the cardinality of a set E by card(E). For a
Pp

vector v ∈ Rd , the `p norm is kvkp =k=1 |v k |p 1/p for p > 0, and kvk∞ = max|v i |.
i

Let thesubset A ⊆ {1, · · · , d}, then v A ∈ Rcard(A) is the vector v restricted to A. Mm×n (R) denotes
the space of m × n matrices with coefficients in R. For a matrix A, kAks is the spectral norm.
We write A0 (resp. v 0 ) to denote the transpose of the matrix A (resp. the vector v). We writevec(A) to denote the vectorization operator that stacks the columns of A on top of one anotherinto a vector. We denote by vech(A) the p(p + 1)/2 vector that stacks the columns of the lower
triangular part of the square and symmetric matrix A. The Ip matrix is the p-dimensional identity
matrix. For a function f : Rd → R, we denote by ∇f the gradient or subgradient of f and ∇2 f
the Hessian of f . We denote by (∇2 f )AA the Hessian of f restricted to the block A. We write Ac
to denote the complement of the set A.

2

Penalised OLS framework for MSV

2.1

Framework

We consider a p-dimensional vectorial stochastic process (yt )t=1,··· ,T and we denote by θ thevector of its model parameters. We then consider a Multivariate Stochastic Volatility (MSV)
decomposition given as
yt = Dt εt ,εt ∼ iid(0, Γ),

ht+1 = µ + Φ(ht − µ) + ηt ,(1)
ηt ∼ NRp (0, Ση ),(2)where Γ is a p×p correlation matrix, εt = (ε1t , . . . , εpt )0 is a p×1 random vector, which is independently and identically distributed (i.i.d.), centered with variance-covariance Γ, ht = (h1t , . . . , hpt )0
is a p × 1 vector of log-volatility, Dt = diag {exp(h1t /2), . . . , exp(hpt /2)} is a diagonal matrix of5

volatility, µ = (µ1 , . . . , µp )0 is a p × 1 vector, Φ is a p × p matrix, and Ση is a p × p covariancematrix of ηt . The MSV model (1) and (2) reduces to the MSV model of Harvey, Ruiz, and Shepard(1994), should we assume Φ diagonal and εit following the t distribution.
2 ), . . . , log(y 2 ))0 . As discussed in Harvey, Ruiz, and Shepard
Then we define xt = (log(y1tpt(1994), the MSV model can be formulated as a state space model, specified as
xt = c + αt + ζt ,(3)

αt+1 = Φαt + ηt ,(4)

where c = (c1 , . . . , cp )0 , ζt = (ζ1t , . . . , ζpt )0 , and αt = ht − µ with ci = µi + E[log(ε2it )] andζit = log(ε2it ) − E[log(ε2it )]. Assuming a t distribution for εit , Harvey, Ruiz, and Shepard (1994)
specified the covariance matrix of ζt as Σζ . Note that E[ζt ] = 0 by definition. Based on the state
space form, these authors then suggested a quasi maximum likelihood estimation of the MSV
model using the Kalman filter. Alternative methods were proposed such as the Bayesian MCMC
technique of Chib, Nardari, and Shephard (2006) or the Monte Carlo Likelihood (MCL) method of
Durbin and Koopman (1997, 2001). A significant drawback of these methods is the computational
cost and thus the curse of dimensionality: most of the applications are restricted to small vector
sizes and/or reduced forms are fostered.
In this paper, we aim at tackling this issue for MSV models using a penalised OLS estimation
method. Although the MSV model (1) and (2) might be a basic model, the following advantages with respect to MGARCH models can be highlighted: (i) relatively stable estimates andforecasts for variance-covariance matrices; (ii) simpler restrictions for stationarity conditions; (iii)
no intricate matrix parameterization and/or parameter restrictions to generate positive-definite
matrices.

62.2

Our proposed approach

In this section, we specify our proposed sparse OLS based MSV model - named as “penalised
OLS-MSV” - and discuss how the latter can manage high-dimensional stochastic vectors. Our
approach can be summarized as follows: derive a VARMA representation of xt , and considera regression-based estimator of (c, Φ), following the ideas of Hannan and Rissanen (1982) and
Hannan and Kavalieris (1984); conditionally on the VARMA estimators, use an ad-hoc estimatorfor Σζ such that the corresponding estimator is positive-definite; finally, obtain the estimator ofΓ.
More precisely, the procedure is as follows. Since xt is the sum of a VAR(1) process and an
i.i.d. noise by (3), the discussion of Granger and Morris (1976) suggests that xt has a VARMA(1,1)
representation. By equations (3) and (4), we obtain
xt = (I − Φ)c + Φxt−1 + (ζt + ηt−1 ) − Φζt−1 ,which has an alternative representation as
xt = (I − Φ)c + Φxt−1 + ut + Ξut−1 ,(5)E[ut ] = 0,

Var(ut ) = Σu ,E[ut u0s ] = 0 for t 6= s,

where Ξ and Σu are obtained by matching moments of wt = (ζt +ηt−1 )−Φζt−1 and wt∗ = ut +Ξut−1 .
0 ] = E[w ∗ w ∗0 ], the relationship between (Ξ, Σ )By considering E[wt wt0 ] = E[wt∗ wt∗0 ] and E[wt wt−1ut t−1

and other parameters is given as follows:Ση + Σζ + ΦΣζ Φ0 = Σu + ΞΣu Ξ0 ,−ΦΣζ = ΞΣu .

(6)
(7)

Given the values of Φ, Ξ, and Σu we may obtain Ση and Σζ by solving equations (6) and (7).

7Under suitable stationarity conditions, xt has an AR(∞) representation:
xt =

∞
X

Ψi xt−i + ut .

(8)

i=1

Based on a penalised OLS estimation, we can obtain a estimator of ût in the first step. Indeed, we
empirically need to specify m sufficiently large as a surrogate of ∞ in the summation in (8). Thus,
for the sake of parsimony and to avoid the overfitting issue, we propose to enforce the nullity ofpossibly numerous model coefficients in the Φi ’s.
In the second step, we calculate the OLS estimator of (ĉ∗ , Φ̂, Ξ̂) by regressing xt on a constant,xt−1 , and ût−1 .
For the third step, we start from the decomposition of the unconditional variance-covariancematrix of xt , which is given byΣx = Σα + Σζ ,(9)where Σx = E[(xt − c)(xt − c)0 ], Σζ = E[ζt ζt0 ], and Σα = E[αt αt0 ] withvec(Σα ) = [Ip2 − (Φ ⊗ Φ)]−1 vec(Ση ).
Denote the sample covariance matrix of xt and ût as Sx and Sû , respectively. An estimator of Σζ
is given as
Sζ = −i1 h −1Φ̂ Ξ̂Sû + Sû Ξ̂0 Φ̂0−1 ,2by equation (7). As there is no guarantee for Sζ and Sx − Sζ to be positive definite, we consider
ado hoc estimators for Σζ and Σα based on decomposition (9).
In the fourth step, we estimate Γ by a correlation matrix of yt .
We summarize our procedure as follows:8

Step 1. Define Ψ1:m = [Ψ1 · · · Ψm ] ∈ Mp×pm (R). Approximate (5) by
xt =

m
X

Ψi xt−i + ut ,

(10)

i=1

and consider the regularised OLS estimator
"
#
T
m
X
1 X
λT
2
||xt −
Ψ̂1:m = arg min
Ψi xt−i ||2 + p( , vec(Ψ1:m )) ,
2T
T
Ψ1:m
t=1

(11)

i=1

where p( λTT , ·) : Rd → R is a penalty function applied to each component of Ψ1:m , where λT
is the regularisation parameter, which depends on the sample size, and enforce a particular
type of sparse structure in the solution Ψ̂1:m . The parameter dimension in the first step isdenoted as d, which in Step 1 is d = mp2 . In this paper, as it will be detailed in Section 3,
we consider the SCAD due to Fan and Li (2001), the MCP due to Zhang (2010), the Lasso
of Tibshirani (1996) and the Bridge of Fu (1998).
Step 2. Define ût = xt −Pm

i=1 Ψ̂i xt−i . Conditionally on Ψ̂1:m , we consider the regression

xt = c∗ + Φxt−1 + Ξût−1 + vt ,where the parameters are (c∗ , Φ, Ξ). Thus, the second step objective function is
"
#
T
 2
1 X
∗
∗
(ĉ , Φ̂, Ξ̂)|Ψ̂1:m = arg min||xt − c + Φxt−1 + Ξût−1 ||2 ,2T(c∗ ,Φ,Ξ)t=1

such that we can obtain the estimator of c by ĉ = (I − Φ̂)−1 ĉ∗ . In this step, the second step
parameter dimension is p(1 + 2p).
Step 3. Consider the estimator of Σζ and Σα as follows:Σ̂α = (1 − r)Sx ,Σ̂ζ = rSx ,(12)

where r is a constant which satisfies 0 < r < 1. This ad hoc method aims at treating
the positive definiteness of the estimators and dealing with the high-dimensionality issue9of Σ̂ζ . While we consider a naive decomposition based on equation (9) for the former, we
set r = (π 2 /2)(p−1

Pp−1 in (12).
i=1 Sx,ii )Here, π 2 /2 is the value of E[ζit2 ] when εit follows

the standard normal distribution. Based on this specification, the average of the diagonal
elements of Σ̂ζ becomes π 2 /2. Using such approach, we are able to accurately estimateΣζ and Σα and importantly the computational cost is negligible, compared to alternative
estimators (e.g. GMM type method) that would require a numerical optimisation.
Step 4. Estimate Γ by a correlation matrix of yt .
When the tuning parameter λT shrinks to zero, Steps 1 and 2 reduce to the standard OLS
estimation for low-dimensional VARMA models, considered by Hannan and Rissanen (1982) and
Hannan and Kavalieris (1984). Step 1 actually corresponds to a multivariate version of the AR(∞)
representation of a log-GARCH model, which is a special case of the exponential GARCH model of
Nelson (1991). Thus, for significantly large m chosen ex-ante, (10) would be a relevant approximation of a log-GARCH type process. Our approach avoids the use of computationally demanding
methods such as MCMC and MCL. Although Harvey, Ruiz, and Shephard (1994) applied theKalman filter, its computational cost is non-negligible for larger p, since the cost increases withthe speed of O(T p2 ) for storing covariance matrices of p × 1 state vector for all t = 1, . . . , T .
For the estimators in Step 3, we may improve them by considering moment-matching methods
using equations (6), (7), and (9) with restrictions on the positive-definiteness of the estimatorsof Σζ and Σα . However, we use the above fast method described in Step 3 without the need ofa numerical optimization procedure. Finally, the fourth step can easily be adapted to a sparse
correlation matrix setting, especially when the size p/T is not negligible.
We now introduce our setting for generating the volatility process. For a low-dimensional case,
we can calculate the minimum mean square linear estimator (MMSLE) of αt based on the full10

sample x = (x01 , . . . , x0T )0 (T p×1) by the state space smoothing algorithm. In the high-dimensionalcase, we consider the multivariate version of the approach of Harvey (1998) with the vector formof (3), as follows:x = c† + α + ζwhere c† = (ιT ⊗ c) (T p × 1), α = (α10 , . . . , αT0 )0 (T p × 1), and ζ = (ζ10 , . . . , ζT0 )0 (T p × 1). By themodel structure, the covariance matrix of x is given by
Vx = Vα + Vζ ,
where



Σα
ΦΣα
Φ2 Σα
..
.

Σα Φ0
Σα
ΦΣα
..
.

Σα (Φ0 )2
Σα Φ0
Σα
..
.





Vα =  ΦT −2 Σα ΦT −3 Σα ΦT −4 ΣαΦT −1 Σα ΦT −2 Σα ΦT −3 Σα· · · Σα (Φ0 )T −2 Σα (Φ0 )T −1· · · Σα (Φ0 )T −3 Σα (Φ0 )T −2 · · · Σα (Φ0 )T −4 Σα (Φ0 )T −3 

,
..
..
..

.
.
.

0

···
Σα
Σα Φ
···
ΦΣα
Σα

and Vζ = (IT ⊗ Σζ ). Then the MMSLE isα̃ = Vα Vx−1 (x − c† ) + c† .

(13)

As in Harvey (1998), we consider an estimator of the covariance matrix, Ht = Dt ΓDt , such thatthe sample variance of the standardized variable of yit equals to one. We consider the estimatorasH̃t = D̃t Γ̂D̃t ,

(14)

where
o
n
D̃t = diag d˜1t , . . . , d˜pt ,d˜it = d¯i exp (x̃it /2) ,

v
u
T
u
X
2 exp (−x̃ ),¯yit
di = tT −1
it
t=1

for i = 1, . . . , p. The standardized variables are defined as z̃it = yit /d˜it , which implies by definition
T −1

PT

2
t=1 z̃it = 1. We call our proposed parameterisation “penalised OLS-MSV”.

11

2.3

Volatility forecasting

We now provide the forecasts for variance-covariance based on our proposed method. The MMSLE
for the lth-step-ahead forecast of αT is given by:α̂T +l = Rl Vx−1 (x − c† ) + c,

(15)

where
h
i
Rl = ΦT +l−1 Σα ΦT +l−2 Σα · · · Φl Σα .
The lth-step-ahead forecast of the covariance matrix is given by:Ĥt = D̂t Γ̂D̂t ,

(16)

where
n
o
D̂T +l = diag dˆ1,T +l , . . . , dˆp,T +l ,dˆi,T +l = d¯i exp (x̂i,T +l /2) ,for i = 1, . . . , p.

3

Asymptotic properties

3.1

First step: penalised estimator Ψ̂1:m

In this section, we focus on the large sample properties of the first step penalised estimator Ψ̂1:m .
We consider a loss function GT from RpT × Θ1 to R. The value GT (y; θ), with θ = vec(Ψ1:m ) ∈Θ1 ⊂ Rd and d = mp2 , evaluates the quality of the “fit” for the realizations of yt for everyt = 1, · · · , T and under Pθ . GT (y; θ) is the empirical loss associated to a continuous function` : RpT × Θ1 → R, idest
T
m
T
T
X
1X
1X
1X
1
1
2
2
GT (y; θ) :=kxt −Ψi xt−i k2 =kxt − Ψ1:m Zm,t−1 k2 :=`(ys , s ≤ t; θ),
T
2
T
2
T
t=1

i=1

t=1

12

t=1

2 ), Ψwhere xt corresponds to the vector of continuous transforms of log(yit1:m = (Ψ1 , · · · , Ψm ) ∈Mp×pm (R) and Zm,t−1 = (x0t−1 , · · · , x0t−m )0 ∈ Rpm . Then, the problem of interest isdXλTp( , |θi |)},
θ̂ = arg min {GT (y; θ) +
T
θ∈Θ1

(17)

i=1

where p( λTT , .) : Θ1 → R is a regulariser/penalty, λT is the regularisation parameter, which
depends on the sample size, and enforce a particular type of sparse structure in the solution.
The function θ → E[`(ys , s ≤ t; θ)] is supposed to be uniquely minimized at θ = θ0 so thatE[∇θ GT (y; θ0 )] = 0.


d
P
λT
The penalisation is performed through the term p T , |θi | , which is a coordinate-separable
i=1

penalty. In this paper, we consider the following set of penalties p(λ, θ), θ ∈ R:
Lasso
: λ|θ|,
Bridge : λ|θ|q , q ∈ (0, 1),
R |θ|
MCP
: sgn(θ)λ 0 (1 − z/(λbmcp ))+ dz,

for |θ| ≤ λ,
λ|θ|,
SCAD

:

2

2

scad λ|θ|+λ )− (θ −2b, for λ ≤ |θ| ≤ bscad λ,(2(bscad −1))(bscad + 1)λ2 /2,
for |θ| > bscad λ,and λ ≥ 0 is the regularisation parameter.
To carry out a sound asymptotic theory, we make the following assumptions.
Assumption 1. card(A) = k0 < d with A = supp(θ) := {i|θ0,i 6= 0}.
Assumption 2. The parameter set Θ1 ⊂ Rd is compact and convex.
Assumption 3. (yt ) is a strictly stationary, non-anticipative and ergodic process.
Assumption 4. For any θ ∈ Θ1 , there exists a measurable function g(.) such that |`(ys , s ≤t; θ)| ≤ g(ys , s ≤ t) with E[g(ys , s ≤ t)] < ∞.
Assumption 5. H = E[∇2θθ0 `(ys , s ≤ t; θ0 )] and M = E[∇θ `(ys , s ≤ t; θ0 )∇θ0 `(ys , s ≤ t; θ0 )] exist
and are positive-definite matrices.
13

These assumptions deserve a few comments. First, assumption 1 corresponds to the sparsity
assumption. As it is unknown, we rely on penalised M-estimation to recover this set. Assumption
2 is standard when analysing large sample properties. For the sake of clarity of our arguments,
we assume 3, which refers to the probabilistic property of the process. Note that it can be relaxedin favor for mixing conditions. Should we consider e.g. strongly mixing time series, where (yt )twould be a strongly mixing process with mixing coefficient α(.) satisfying α(ς) ≤ κρς with ς > 0
and 0 < ρ < 1, then we would rely on exponential bound adapted to the strongly mixing caseto control each element of the Taylor expansions we propose to analyse. To do so, Theorem 2 ofMerlevède, Peligrad and Rio (2009), could be used under strongly mixing and bounded randomvariables. Central limit theorem for strongly mixing processes, such as e.g. Theorem 5.2 of White(2001), could also be applied for that case. We highlight that relaxing such assumption for mixing
processes would not alter the convergence rates we propose to derive. Assumption 4 corresponds toa domination condition on the likelihood function, which will be used when proving the consistencyof the penalised estimator. Finally, assumption 5 is a standard regularity condition.
Theorem 1. Under assumptions 1-4, if λTT → λ0 , then for any compact B ⊂ Θ1 such that
θ0 ∈ B,dXPθ̂ −→ arg min{Gpen
(y;
x)}
:=
arg
min
{G
(y;
x)
+
p(λ0 , |xi |)} = θ0∗ ,
∞
∞
T →∞
x∈B
x∈B
i=1

with G∞ (y; x) = E[`(ys , s ≤ t; x)] and for any scalar x
Lasso : p(λ0 , |x|)= λ0 |x|,
Bridge : p(λ0 , |x|) = λ0 |x|q ,MCP : p(λ0 , |x|)= bmcp λ20 /21{|x|>bmcp λ0 } − (bmcp λ0 − |x|)2 /(2bmcp )1{|x|≤bmcp λ0 } ,
for |x| ≤ λ0 ,
λ0 |x|,2
2
SCAD : p(λ0 , |x|) =−(x − 2bscad λ0 |x| + λ0 )/(2(bscad − 1)), for λ0 ≤ |x| ≤ bscad λ0 ,(bscad + 1)λ20 /2,
for |x| > bscad λ0 ,Hence if λT = o(T ), then θ̂ is a consistent estimator.
14

Remark. The penalised estimator does not converge to θ0 when λT = O(T ) due to the bias
originating from the penalisation. Part of the proof consists in proving a uniform convergence ofthe penalised criterion. To do so, we rely on Lemma 2.4 of Newey and McFadden (1994) to derivea uniform law of large numbers of the parameter dependent loss function.

We now propose an asymptotic probabilistic bound for the penalised estimator. First, we
assume the following conditions on the penalty function.
Assumption 6. p( λTT , |.|) is twice continuously differentiable except at the origin. We define
A1,T = max|∇θi p(
i∈A

λT
λT
, |θ0,i |)|, A2,T = max|∇2θi θi p( , |θ0,i |)|,i∈A
T
T

so that A2,T → 0.
Remark. The condition on the second derivative implies that the penalty has less influence than
the non-penalised loss function in the regularised problem. Moreover, for the penalties of interest,
the scaling of (λT , T ) determines this rate.
Theorem 2. Under assumption 1-6, the sequence of penalised estimators θ̂ satisfieskθ̂ − θ0 k = Op (T −1/2 +pcard(A)A1,T ).

Remark. This result highlights that for a suitable choice of λT , we would obtain a√
T -consistentθ̂. This probability bound is similar to the so-called oracle bounds as it depends on the true sparse
support card(A).
√We now derive the asymptotic distribution for the rate λT = O( T ) for the Lasso, SCAD and
MCP and λT = O(T q/2 ) in the Bridge case.

15

Theorem 3. Under assumptions 1-6, suppose λT = o(T ), if the regularisation rate of the Lasso,√
SCAD and MCP satisfies λT = O( T ) and the Bridge regularisation rate satisfies λT = O(T q/2 ),if lim ∇x p( λTT , x) = λTT for the SCAD and MCP, then
x→0+

√

d

T (θ̂ − θ0 ) −→ arg min {F∞ (u)},
T →∞
u∈Rd

provided F∞ (.) is the random function in Rd where1
F∞ (u) = u0 W + u0 Hu + φ(λ0 , u, θ0 ),2
where W ∼ NRd (0, M) with M := M(θ0 ) = E[∇θ `(ys , s ≤ t; θ0 )∇θ0 `(ys , s ≤ t; θ0 )], H = E[∇2θθ0 `(ys , s ≤t; θ0 )] andLasso : φ(λ0 , u, θ0 )= λ0Bridge : φ(λ0 , u, θ0 ) = λ0
MCP : φ(λ0 , u, θ0 )SCAD : φ(λ0 , u, θ0 )

= λ0
= λ0

d
P


ui sgn(θ0,i )1θ0,i 6=0 + |ui |1θ0,i =0 ,

k=1
d
P

|ui |q 1θ0,i =0

k=1
d
P

|ui |1θ0,i =0 ,

k=1
d
P

|ui |1θ0,i =0 .

k=1Remark. The following comments can be noticed.
(i) To derive such distributions, we rely on specific theoretical results depending on the penalty
function and thus on the (non-)convexity of the latter. These results are reported in Appendix B.
√(ii) Theorem 3 establishes the T -consistency of the penalised estimator. However, for λT =√O( T ) in the Lasso case, the term in 1θ0,i 6=0 implies that the true active set A can not
be recovered with high probability (See Proposition 1 of Zou, 2006). To fix this issue, Zou(2006) proposed the adaptive Lasso, which consists in penalising the coefficients differently
through the introduction of stochastic weights. These weights are explicit functions of a16

first step estimator, which is√

T -consistent. As a consequence, these weights alter theconvergence rate of the penalisation parameter λT .
We now turn to the oracle property. It is well-known since Zou (2006) that the Lasso does
not satisfy this property. The only method to fix this issue is to specify adaptive weights in the
penalty function to penalise each coefficient differently. The key advantage of non-convex penalties(SCAD, MCP, Bridge) is that they actually allow for satisfying this property without the need ofthese stochastic weights.
Theorem 4. Suppose λT = o(T ), for θA satisfying kθA −θ0,A k = Op (T −1/2 ), suppose assumptionsT
1-6, suppose the MCP and SCAD regularisation rates satisfy Tλ1/2→ ∞ and

lim inf lim inf
T →∞

x→0+

λT
T
∇x p( , x) > 0,
λT
T

√
T
→ ∞, 0 < q < 1 and λT = O( T ),and suppose the Bridge satisfies the regularisation rate Tλq/2√
then the T -consistent local estimator θ̂ defined in Theorem 2 satisfies
lim P(Â = A) = 1, and

T →∞

√


d
−1
T θ̂ − θ0 A −→ NRk0 (0, H−1
AA MAA HAA ).
T →∞Remark. Theorem 4 deserves a few comments.
(i) This result establishes the conditions to satisfy the oracle property. Contrary to Theorem√
3, where the rate for the SCAD/MCP is λT = O( T ) and for the Bridge is λT = O(T q/2 ),√we now require λT / T → ∞ for those ones and λT /T q/2 → ∞ for the latter. Note thatthe adaptive Lasso is left aside as we report sparse estimation methods that do not requirea two-step estimation. More details can be found in Zou (2006) on the adaptive Lasso andits properties.

17

(ii) The non-random quantities originating from the expansion of the penalty functions
0bT,A= ∇θ1 p( λTT , |θ0,1 |)sgn(θ0,1 ), · · · , ∇θk0 p( λTT , |θ0,k0 |)sgn(θ0,k0 ) ,ST,AA = diag(∇2θi θi p( λTT , |θ0,i |), i = 1, · · · , k0 ),vanish when T is large enough due to the assumed convergence rates of the penalisation
parameter λT for the Bridge, SCAD and MCP. A detailed discussion is provided in the
proof.

3.2

Second step estimator (ĉ∗ , Φ̂, Ξ̂)We now focus on the large sample properties of the second step estimator γ̂ = (ĉ∗ , vec(Φ̂)0 , vec(Ξ̂))0 ,which is of size p(1 + 2p). Conditionally on θ̂, we consider a second step loss function LT fromRpT × Θ2 to R with Θ2 ⊂ Rp(1+2p) and LT (y; θ̂, γ) is the empirical loss associated to a continuousfunction f : RpT × Θ1 × Θ2 → R, idest
LT (y; θ̂, γ) =

T

T

t=1

t=1

1 X1
1X
kxt − α − Φxt−1 − Ξût−1 k22 :=f (ys , s ≤ t; θ̂, γ).
T
2T

The problem of interest isγ̂ = arg min {LT (y; θ̂, γ)}.
γ∈Θ2To derive the large sample properties of the second step estimator, we make the following assumptions.
Assumption 7. The second step parameter set Θ2 is compact.
Assumption 8. For any θ ∈ Θ1 and γ ∈ Θ2 , there exists a measurable function k(.) such that|f (ys , s ≤ t; θ, γ)| ≤ k(ys , s ≤ t) with E[k(ys , s ≤ t)] < ∞.
Assumption 9. For any (vec(Φ)0 , vec(Ξ)0 , vec(Ψ1:m )0 )0 ∈ Θ1 × Θ2 , the pm × pm matrixΦ + Ξ −ΞΨ1 −ΞΨ2 · · · −ΞΨm
 Ip
0
0
···
0 


 0

I
0
·
·
·
0
p
Λ=
,
 ..
.. 
..
 .
.
···
···
. 
0

···

···
18

Ip

0

satisfies kΛks < 1.
Assumption 9 is key to control for the first step estimator and allows for establishingsup |LT (y; θ̂, γ) − LT (y; θ0∗ , γ)| = op (1),γ∈Θ2where θ̂ is a sequence in Θ1 that tends to θ0∗ in probability.
Theorem 5. Let β̂ = (θ̂0 , γ̂ 0 )0 be a sequence of penalised OLS based and non-penalised OLS based0estimators. Then, under the assumptions of Theorem 1 and assumptions 7-9, β̂ −→ (θ0∗ , γ00 )0 .
P

T →∞We now turn to the asymptotic distribution of the second step estimator, conditionally on θ̂.
Theorem 6. Under the assumptions of Theorem 4 together with its conditions on the penalisation
rate on λT and the penalty function’s behaviour, under assumptions 7-9, then√T(θ̂ − θ0 )Aγ̂ − γ0d−→ NRdim (0, J−1 IJ−1 ),T →∞with dim = k0 + p(1 + 2p) and the variance-covariance is composed withE[∇2θθ0 `(ys , s ≤ t; θ0 )]AA0E[∇2θγ 0 f (ys , s ≤ t; θ0 , γ0 )]A• E[∇2γγ 0 f (ys , s ≤ t; θ0 , γ0 )]


J =
and

I = E





∇θ `(ys , s ≤ t; θ0 )A ∇θ0 `(ys , s ≤ t; θ0 )A∇θ `(ys , s ≤ t; θ0 )A ∇γ f (ys , s ≤ t; θ0 , γ0 ).
∇γ f (ys , s ≤ t; θ0 , γ0 )∇θ0 `(ys , s ≤ t; θ0 )A ∇γ f (ys , s ≤ t; θ0 , γ0 )∇γ 0 f (ys , s ≤ t; θ0 , γ0 )Remark. Theorem 6 deserves some comments.
(i) The notation E[∇2θγ 0 f (ys , s ≤ t; θ0 , γ0 )]A• means that the expectation is restricted only toblock A when evaluating the partial derivative with respect to θ ∈ Θ1 .
(ii) The effect of the two-step estimation can be measured through the asymptotic variancecovariance matrix. As a by-product, simple calculations provide the asymptotic variances
19

−1
−1
V(.) of θ̂ and γ̂: with obvious notations on the block locations in J and I, V(θ̂) = J11
I11 J11
,

and
−1
−1
−1
−1
−1
−1
V(γ̂) = J22
I22 J22− ΥI12 J22− J22
I21 Υ0 + ΥI11 Υ0 , Υ := J22
J21 J11.

4

Empirical analysis

4.1

Simulation experiment

In this section, we empirically investigate the ability of the proposed penalisation method to better
capture complex variance-covariance processes. We simulate the p-dimensional stochastic process(t ) based on two data generating processes: the multivariate ARCH and the BEKK processes.
For the multivariate ARCH with q ∗ lags - M-ARCH(q ∗ ) in the rest of the paper - case, we consider
1/2 t = Ht ηt ,q∗P(Ip ⊗ 0t−k )Ak (Ip ⊗ t−k ), Ht = Ω +k=1where q ∗ is the number of lagged matrices being functions of t−k and the p2 ×p2 square matrices Aksatisfy the stationarity conditions of Theorem 2 of Boussama (2006) together with the positivity
condition given by Gouriéroux (1997). We generate the diagonal elements of Ak from a uniform
distribution U([0.01, 0.05]) and the off-diagonal ones from U([−0.01, 0.01]) under the ordering
constraint ∀k ≥ 2, ∀i, j, |Ak,ij | ≤ |Ak−1,ij |. As for the matrix Ω, the diagonal and off-diagonal
elements are simulated from U([0.1; 0.2]) and U([−0.01, 0.01]) respectively. As for the BEKK
process, the data generating process is based on(

1/2t = Ht ηt ,Ht = Ω + At−1 0t−1 A0 + BHt−1 B 0 ,
where A, B are p×p matrices, satisfying the stationarity constraint kDp+ { A⊗A + B⊗B }Dp ks <1, where Dp is the duplication matrix and Dp+ the elemination matrix (see subsection 11.3 ”Stationarity of VEC and BEKK Models” of Francq and Zakoı̈an (2010) for the stationarity condition
20and remark 11.1 for the definition of the latter matrices). The entries of A and B are generated
from the uniform distribution U([−0.8, 0.8]). The matrix Ω is generated as in the M-ARCH(q ∗ )case. Unlike the M-ARCH(q ∗ ) case, the BEKK dynamic includes an autoregressive component
through B, which motivated the use of larger lags when estimating our proposed parameterizations. In both proposed dynamics, we initialize the observations (k , · · · , 1 ) with centered andunit variance multivariate Gaussian distribution, where k = q ∗ in the M-ARCH model and k = 1
in the BEKK. Then conditionally on the past k observations, we generate Ht and thus t according
to a centered multivariate Gaussian distribution with variance-covariance Ht .
We consider the problem sizes, p = 15, 50, and T = 2000 observations for each of them. For
the M-ARCH(q ∗ )-based data generating process, we considered q ∗ = 2 when p = 15 and q ∗ = 1
when p = 50. Then we propose to compare the true variance-covariance processes - BEKK and
M-ARCH(q ∗ ) - and the estimated ones through our proposed MSV model and the scalar DCC
corresponding to process (25) with scalar matrix parameters together with the constant correlation
model (CCC). The estimation of the DCC model is based on the classic two-step Gaussian QMLE,
where the marginal conditional volatility processes are specified as GARCH(1,1) and a correlation
targeting procedure is applied in the second step, providing an estimated trajectory Ĥtdcc . The
CCC is estimated thanks to a joint estimation of the GARCH(1,1) parameters and correlation
parameters through a Gaussian QML, which provides an estimated process Ĥtccc . More details on
the DCC and CCC can be found in Appendix D.
Regarding our proposed variance-covariance dynamic, the penalised OLS-MSV, denoted asĤtols,al for the adaptive Lasso OLS-MSV, Ĥtols,br for the Bridge OLS-MSV, Ĥtols,scad for the SCAD
OLS-MSV, Ĥtols,mcp for the MCP OLS-MSV, and the non-penalised version of the OLS-MSV
denoted as Ĥtols . We considered the same technique as in Zou (2006) for the adaptive Lasso,21

where we selected γ = 1.5 as the power entering in the stochastic weights and the first step
estimator is the unpenalised OLS estimator. In the M-ARCH(q ∗ ) case, we set the number of lagsin Step 1 in (10) as m = 10 when p = 15 and set as m = 5 for a dimension p = 50. In the BEKK
case, due to the autoregressive nature of the latter dynamic, we selected m = 30 when p = 15 andm = 15 when the dimension is p = 50.
We compare the true variance-covariance process and the estimated variance-covariance processes through the aforementioned models. To do so, we specify a matrix distance, namely thep
Frobenius norm, defined as ||A − B||F := Trace((A − B)0 (A − B)). We compute the previous
norm for each t and for A = Ht andB ∈ {Ĥtdcc , Ĥtccc , Ĥtols , Ĥtols,br , Ĥtols,al , Ĥtols,scad , Ĥtols,mcp }.
We take the average of those quantities over T = 2000 periods of time. Since we repeat thisexperiment 100 times, this provides an average gap for all those simulations.
By a cross-validation (CV) procedure - see e.g. Hastie and al. (2015, Chap. 2) - , we selected
the regularisation parameter and emphasize that the standard CV developed for i.i.d. data can
not be used in our time series framework. To fix this issue, we used the hv-CV procedure devisedby Racine (2000), which consists in leaving a gap between the test sample and the training sample,
on both sides of the test sample.
Clearly, the relevance of the penalisation procedure reported in Table 1 for the M-ARCH(q ∗ )based DGP and Table 2 for the BEKK based DGP increases with the size p, for any penalisation
method. Moreover, DCC/CCC models are always beaten by the penalised specifications, whateverp, whereas, interestingly, it is not the case when compared with the non-penalised MSV. Our
results emphasize the clear gain when considering penalisation, especially in the case m, p large as22

in the BEKK case: the autoregressive nature of such process fostered the use of a large numberof lags, that is a large m.

4.2

Application to real data

To assess the relevance of the proposed penalised method, we propose a real data experiment.
To do so, we compare the forecasting performances of the covariance matrices Ht for a portfolio of daily financial returns composed of the MSCI stock index for the following 23 countries:
Australia, Austria, Belgium, Canada, Denmark, Finland, France, Germany, Greece, Hong Kong,
Ireland, Italy, Japan, Netherlands, New Zealand, Norway, Portugal, Singapore, Spain, Sweden,
Switzerland, the United-Kingdom, the United-States. We focus on direct out-of-sample evaluation
methods, which allow for pairwise comparisons. They test whether some of the variance-covariancemodels provide better forecasts in terms of portfolio volatility behavior. Following the methodology of Engle and Colacito (2006), we develop a mean-variance portfolio approach to test the
Ht forecasts. Intuitively, if a conditional covariance process is misspecified, then the minimum
variance portfolio should emphasize such a shortcoming, compared to other models. Then, consider an investor who allocates a fixed amount between p stocks, according to a minimum-variancestrategy and independently at each time t. At each date t, he/she solves
min wt0 Ht wt , s.t. ι0 wt = 1,
wt

(18)

where wt is the p × 1 vector of portfolio weights chosen at (the end of) time t − 1, ι is a p × 1
vector of 1 and Ht is the estimated conditional covariance matrix of the asset returns at time t.
They are deduced from some dynamics that have been estimated on the sub-sample December
1998 - November 2015. Once the latter process is estimated in-sample, out-of-sample predictions
are plugged into the program (18) between December 2015 and March 2018. The solution of (18)
is given by the global minimum variance portfolio wt = Ht−1 ι/ι0 Ht−1 ι.
23Engle and Colacito (2006) show that the realized portfolio volatility is the smallest one whenthe variance-covariance matrices are correctly specified. As a consequence, if wealth is allocatedusing two different dynamic models i and j, whose predicted covariance matrices are (Hti ) and(Htj ), the strategy providing the smallest portfolio variance will be considered as the best one. To
do so, we consider a sequence of minimum variance portfolio weights (wi,t ) and (wj,t ), dependingon the model. Then, we consider a distance based on the difference of the squared returns of the
n
o2 n
o2
0 
0 
two portfolios, defined as uij,t = wi,t−w
. The portfolio variances are the same if thetj,t t
predicted covariance matrices are the same. Thus we test the null hypothesis H0 : E [uij,t ] = 0
by the Diebold and Mariano (1995) test. It consists of a least squares regression using HAC
standard errors, given by uij,t = α + u,t , E[u,t ] = 0, and we test H0 : α = 0. If the mean of uij,tis significantly positive (resp. negative), then the forecasts given by the covariance matrices ofmodel j (resp. i) are preferred.
We run the latter test to compare the scalar DCC (DCC), the Orthogonal GARCH (O-G), the
BEKK (BEKK), our OLS-MSV method (MSV) together with its penalised counterpart (denoted
as MSV-AL, MSV-BR, MSV-SCA, MSV-MCP for the adaptive Lasso, Bridge, SCAD, MCP respectively). The definition of the BEKK and O-GARCH processes are reported in Appendix D.
The matrix forecast comparisons are provided in Table 3. The results emphasises that the proposed penalised OLS-MSV method outperforms the MGARCH based competitors. Interestingly,
fostering sparsity yields to much better forecasting performances. Indeed, the non-penalised OLSMSV is outperformed by all alternative specifications, whereas the sparsity based specifications
outperform both MGARCH based models and the non-penalised OLS-MSV.

24

5

Conclusion

The focus of this paper was devoted to the estimation of high-dimensional MSV models. Our
main contribution consisted in proposing an estimation framework that does not rely on standard
MCMC/MCL methods but instead on a penalised OLS framework for state-space estimation. The
corresponding large sample properties of the two-step estimator are derived for a broad range ofpenalty functions. The performances of our proposed method compared to standard MGARCH
models are illustrated through simulated experiments together with an out-of-sample analysis forprediction accuracy, where our method clearly outperformed the competing MGARCH models.
These results also emphasized the gain of penalisation, which manages the overfitting problem.
Various issues and extensions can be further considered. Our proposed model could be extended to a dynamic correlation setting and include long memory and/or asymmetry, as discussed
in Asai, McAleer, and Yu (2006) and Chin, Omori, and Asai (2009). Another direction would
consists in modelling directly the variance-covariance matrix Ht without relying on the decomposition Dt ΓDt . To do so, a log-type dynamic on Ht could be considered and the estimation could
be managed through the development of a state-space based setting.

25

References
Abadir, K.M. and Magnus, J.R. (2005), Matrix algebra. Cambridge University Press.
Alexander, C. (2001). ”Orthogonal GARCH” in Alexander, C. (Ed.), Mastering Risk, Financial TimesPrentice Hall, London, pp. 21-28.
Asai, M., M. Caporin, M. McAleer (2015), “Forecasting Value-at-Risk Using Block Structure Multivariate
Stochastic Volatility”, International Review of Economics & Finance, 40, 40–50.
Asai, M. and M. McAleer (2009), “Multivariate Stochastic Volatility, Leverage and News Impact Surfaces”,
Econometrics Journal, 12, 292–309.
Asai, M., M. McAleer, and J. Yu (2006), “Multivariate Stochastic Volatility: A Review”, Econometric
Reviews, 25, 145–175.
Baba, Y., R. Engle, D. Kraft and K. Kroner (1985), “Multivariate Simultaneous Generalized ARCH”,
Unpublished Paper, University of California, San Diego. [Published as Engle and Kroner (1995)]
Bauwens L., S. Laurent, and J. K. V. Rombouts (2006), “Multivariate GARCH Models: A Survey”,
Journal of Applied Econometrics, 21, 79–109.
Billingsley, P. (1961), “The Lindeberg-Levy theorem for martingales”, Proceedings of the American Mathematical Society, 12, 788–792.
Billingsley (1995), Probability and measure, New York: John Wiley&Sons.
Boussama, F. (2006), ”Ergodicité des chaı̂nes de Markov à valeurs dans une variété algébrique: application
aux modèles GARCH multivariés”, Comptes Rendus de l’Académie des Sciences Paris, 343, 275–278.
Chernozhukov, V. and Hong, H. (2004), ”Likelihood estimation and inference in a class of nonregulareconometric models”, Econometrica, 72, No. 5, 1445-1480.
Chernozhukov, V., (2005), ”Extremal quantile regression”, The Annals of Statistics, 33, No. 2, 806-839.
Chib, S., F. Nardari, and N. Shephard (2006), “Analysis of High Dimensional Multivariate Stochastic
Volatility Models”, Journal of Econometrics 134, 341–371.
Chib, S., Y. Omori, and M. Asai (2009), “Multivariate Stochastic Volatility”, In: Andersen, T.G., R.A.
Davis, J.P. Kreiss, and T. Mikosch, T. (Eds.), Handbook of Financial Time Series, New York:
Springer-Verlag, pp.365–400.
Durbin, J., and S.J. Koopman (1997), “Monte Carlo Maximum Likelihood Estimation for Non-GaussianState Space Models”, Biometrika, 84, 669–684.
Durbin, J., and S.J. Koopman (2001), Time Series Analysis by State Space Methods, Oxford: Oxford
University Press.
Engle, R. F. (1982), “Autoregressive Conditional Heteroskedasticity with Estimates of the Variance ofUnited Kingdom Inflation”, Econornetrica, 50, 987–1007.
Engle, R. (2002). ”Dynamic Conditional Correlation”, Journal of Business & Economic Statistics 20:
339-350.
Engle, R., and R. Colacito. (2006). ”Testing and Valuing Dynamic Correlations for Asset Allocation”,
Journal of Business & Economic Statistics 24: 238-253.
Engle, R.F. and Kroner, K.F. (1995). ”Multivariate Simultaneous Generalized ARCH”, Econometric
Theory 11: 122-150.

26

Engle, R.F., O. Ledoit, and M. Wolf (2017), “Large Dynamic Covariance Matrices”, Journal of Business& Economic Statistics, 37, 363–375.
Fan, J. and Li, R. (2001), “Variable selection via nonconcave penalised likelihood and its oracle properties”,
Journal of the American Statistical Association, 96, No. 456, 1348-1360.
Fan, J., Y. Fan, and J. Lv. (2008). “Large dimensional covariance matrix estimation using a factor
model”, Journal of Econometrics, 147, 186–197.
Francq, C., and Zakoı̈an, J.-M. (2010), GARCH models structure, statistical inference and financial applications. John Wiley and Sons, Chichester, West Sussex, U.K.
Fu, W. J. (1998). “penalised regressions: The bridge versus the Lasso.” Journal of Computational andGraphical Statistics, 7, 397–416.
Gouriéroux, C. (1997), ”ARCH Models and financial Applications”. Springer.
Granger, C. W. J. and M. Morris (1976), “Time Series Modeling and Interpretation”, Journal of the Royal
Statistical Society, Series A, 139, 246–257.
Ghysels, E., A. C. Harvey, and E. Renault (1996), “Stochastic volatility”, In: Rao, C. R., Maddala,
G. S., eds. Statistical Models in Finance (Handbook of Statistics), Amsterdam: North-Holland, pp.
119–191.
Hannan, E. J. and L. Kavalieris (1984), “Multivariate Linear Time Series Models”, Advances in Applied
Probability,16, 492–561.
Hannan, E. J. and J. Rissanen (1982), “Recursive Estimation of Mixed Autoregressive-Moving Average
Order”, Biometrika, 69, 81–94.
Harvey, A. (1998), “Long Memory in Stochastic Volatility”, In: Knight, J. and S. Satchell (eds.), Forecasting Volatility in Financial Markets, Oxford: Butterworth-Haineman, 307–320.
Harvey, A. C., E. Ruiz, and N. Shephard (1994), “Multivariate Stochastic Variance Models”, Review of
Economic Studies, 61, 247–264.
Hastie,T., R. Tibshirani, and Wainwright, M. (2015), ”Statistical Learning with Sparsity: The Lasso andGeneralizations”, Monographs on Statistics and Applied Probability 143. Chapman and Hall.
Hjort, N. and Pollard, D. (1993), ”Asymptotics for minimisers of convex processes”, Statistical Research
Report, Institute of Mathematics, University of Oslo.
Kastner, G., S. Frü-Schnatter, and H. F. Lopes (2017), “Efficient Bayesian Inference for Multivariate
Factor Stochastic Volatility Models”, Journal of Computational and Graphical Statistics, 26, 905–
917.
Kim, J. and Pollard, D. (1990), ”Cube root asymptotics”, The Annals of Statistics, 18, No. 1, 191-219.
Knight, K. and Fu, W. (2000), “Asymptotics for Lasso-type estimators”, The Annals of Statistics, 28,
No. 5, 1356-1378.
Merlevède, F., Peligrad, M. and Rio, E. (2009). “Bernstein inequality and moderate deviations understrong mixing conditions”, Institute of Mathematical Statistics Collections, High Dimensional Probability, 5, 273–292.
Newey,W. K. and D. McFadden (1994), “Large Sample Estimation and Hypothesis Testing”, In R. F.
Engle and D. McFadden (Eds.), Handbook of Econometrics, Volume 4, 2111–2245. Amsterdam:
Elsevier.

27

Poignard, B. and J.D. Fermanian (2019), “High-dimensional penalised ARCH processes”, To appear in
Econometric Reviews.
Racine, J. (2000), “Consistent cross-validatory model-selection for dependent data: hv-block cross-validation”,
Journal of Econometrics. 99: 39-61.
Tibshirani, R. (1996). “Regression shrinkage and selection via the Lasso”, Journal of the Royal StatisticalSociety. Series B, 58, No. 1, pp. 267-288.
Tse, Y. K. and A. K. C. Tsui (2002), “A Multivariate Generalized Autoregressive Conditional Heteroscedasticity Model with Time-Varying Correlations”, Journal of Business & Economic Statistics,
20, 351–361.
Umezu, Y., Shimizu, Y., Masuda, H. and Ninomiya, Y., (2018), ”AIC for the non-concave penalised
likelihood method”, Annals of the Institute of Statistical Mathematics, 71, No. 2, 247-274.
White, H. (2001), “Asymptotic Theory for Econometricians”, 2nd edition, Emerald, UK.
Zhang, C.-H. (2010). “Nearly unbiased variable selection under minimax concave penalty.” The Annalsof Statistics, 38, 894-942.
Zou, H. (2006). “The adaptive Lasso and its oracle properties”, Journal of the American Statistical
Association, 101, No. 476, 1418-1429.

28

A

TablesTable 1: Average distance true/estimated covariance matrices - M-ARCH(q ∗ ) (100 replications)p = 15p = 50Ĥtdcc
9.12
13.27

Ĥtccc
9.89
14.25

Ĥtols
10.11
15.37

Ĥtols,br
8.41
12.13

Ĥtols,al
8.63
12.44

Ĥtols,scad
8.62
12.66

Ĥtmcp8.46
12.54

Table 2: Average distance true/estimated covariance matrices - BEKK (100 replications)p = 15p = 50

DCC
O-G
BEKK
MSV
MSV-ALMSV-BR
MSV-SCAMSV-MCPĤtdcc14.76
56.93

Ĥtccc14.91
61.98

Ĥtols
18.06
98.235Ĥtols,br
14.24
54.21

Ĥtols,al
14.07
54.86

Ĥtols,scad
14.38
54.15

Ĥtmcp14.12
54.08

Table 3: Diebold Mariano Test of Multivariate Variance-Covariance models
DCC
O-G
BEKK
MSV
MSV-AL MSV-BR MSV-SCA MSV-MCPcc−2.351−1.046 −3.415
7.671c7.587c7.673c7.577c
c
b
c
c
c
2.351
2.014
1.184
6.339
6.781
6.451
6.103c1.046−2.014b−2.092b8.234c8.175c8.241c8.164c
c
b
b
b
c
3.415
−1.184
2.092
1.895
1.982
2.445
2.241bcccb−7.671−6.339−8.234−1.895−0.497−0.196−0.484cccba−7.587−6.340−8.175−1.982
0.497
0.750
0.047
−7.673c −6.364c −8.241c −2.445c
0.196−0.750−0.721−7.577c −6.341c −8.164c −2.241c
0.484−0.047
0.721

Table 4: This table reports the out-of-sample t-statistics of the Diebold-Mariano test that checks
the equality between covariance matrix forecasts using the loss function uij,t over the period
December 2015 and March 2018. This loss function is defined as the difference between squared
realized returns of alternative Multivariate Variance-Covariance models. When the null hypothesisof equal predictive accuracy is rejected, a positive number is evidence in favour of the model inthe column. a, b, c: rejection of the null hypothesis at 10%, 5% and 1% respectively.

29

B

Technical results

We used the convexity argument to derive the asymptotic distribution in the Lasso case. Chernozhukov and Huong (2004), Chernozhukhov (2005) use this convexity argument to obtain theasymptotic distribution of quantile regression type estimators. This argument relies on the convexity Lemma, which is a key result to obtain an asymptotic distribution when the objectivefunction is not differentiable. It only requires the lower-semicontinuity and convexity of the empirical criterion. The convexity Lemma, as in Chernozhukov (2005), proof of Theorem 4.1, can bestated as follows.
Lemma 1. Convexity Lemma, Chernozhukov (2005)
Suppose(i) a sequence of convex lower-semicontinuous FT : Rd → R̄ marginally converges to F∞ : Rd → R̄over a dense subset of Rd ;(ii) F∞ is finite over a nonempty open set E ⊂ Rd ;(iii) F∞ is uniquely minimized at a random vector u∞ .
Then
d

d

arg min FT (z) −→ arg min F∞ (z), that is uT −→ u∞ .
z ∈Rd
z ∈Rd
As for the SCAD and MCP, due to the non-convexity of the penalty function, we used Lemma
3 of Umezu et al. (2018), which generalises Lemma 2 of Hjort and Pollard (1993) to the case of
convex non-penalised loss functions with non-convex penalties. Lemma 2 of Hjort and Pollard(1993) allows for deriving consistency and asymptotic normality of estimators that are defined byminimisation of convex criterion functions.
Lemma 2. Umezu, Shimizu, Masuda, Ninomiya (2018)30

Suppose that GT (u) is a strictly convex random function that is approximated by G̃T (u). Let ū
be a subvector of u, and let ζ(u) and η(ū) be continuous functions such that ζT (u) and ηT (ū)
converge to ζ(u) and η(ū) uniformly over u and ū in any compact set, respectively, and assume
that ζ(u) is convex and η(0) = 0. In addition, forνT (u) = GT (u) + ζT (u) + ηT (ū), and ν̃T (u) = G̃T (u) + ζ(u) + η(ū),
let uT and ũT be the arg min of νT (u) and ν̃T (u), respectively, and assume that ũT is unique and˜ T = 0. Then, for any  > 0, δ > 0, µ > δ, there exists γ > 0 such thatūP(kuT − ũT k ≥ δ) ≤ P(2∆T (δ) +  ≥ ΥT (δ)) + P(kuT − ũT k ≥ µ) + P(kūT k ≥ γ),
where∆T (δ) =sup|νT (u) − ν̃T (u)|, ΥT (δ) =inf|ν̃T (u) − ν̃T (ũT )|.
u:ku−ũT k=δ
u:ku−ũT k≤δFinally, for the large sample distribution of the Bridge penalised estimator, we relied on Theorem 2.7 of Kim and Pollard (1990).
Theorem 7. Kim and Pollard (1990)
Let {FT } be a random function into the space of all locally bounded real functions on Rd , and uT
random mapps into Rd such thatd(i) FT −→ Q for a Borel measure Q concentrated on Cmax (Rd )1 ;(ii) uT = Op (1);(iii) FT (uT ) ≥ sup {FT (u)} − αT for random variables (αT ) of order op (1).
udThen uT −→ arg max {F(u)} for a F(u) with distribution Q.
u
1

See their page 195 for a definition of this set

31

C

Proofs

Proof of Theorem 1. In a first step, we prove the uniform convergence of Gpen
T (y; .) to the limit
quantity Gpen
∞ (y; .) on any compact set B ⊂ Θ, idest
pen0.
sup |Gpen
T (y; x) − G∞ (y; x)| T−→
→∞
x∈B
P

(19)

We define C ⊂ Θ an open convex set and pick x ∈ C. Then, under assumptions 2 and 4, byLemma 2.4 of Newey and McFadden (1994), we havesup |GT (y; θ) − E[`(θ; ys , s ≤ t)]| = 0,x∈B
where the convergence of the sample criterion to the population level criterion is ensured by theergodic Theorem of Billingsley (1995), which thus justifies assumption 3. Thus, using λT /T → λ0and since the parameter is taken over a compact set, we obtainpen
0.
sup |Gpen
T (y; x) − G∞ (y; x)| T−→
→∞
x∈B
P

Now we would likearg min {GT (y; .) +

d
X
i=1

d

p(

X
λT
P
, |.|)} −→ arg min {G∞ (y; .) +p(λ0 , |.|)}.
T →∞
T
i=1

d
P

p( λTT , |θi |) ≥ GT (y; x), and arg min {GT (y; x)} = Op (1) by convexity of
i=1
x∈B
d
P
the criterion, it thus follows that arg min {GT (y; x) +p( λTT , |xi |)} = Op (1) with probability
i=1
x∈B
one.
First, GT (y; x) +Proof of Theorem 2. Let νT = T −1/2 +pcard(A)A1,T . We would like to prove that for any  > 0,

there exists C > 0 such thatP(1kθ̂ − θ0 k > C ) < .
νT
32

Following the reasoning of Fan and Li (2001), Theorem 1, and denoting Gpen
T (y; θ) = GT (y; θ) +dPp( λTT , |θi |), we have

i=1

P(

1
pen
kθ̂ − θ0 k > C ) ≤ P(∃u, kuk2 = C : Gpen
T (y; θ0 + νT u) ≤ GT (y; θ0 )),νTwhich implies that there is a local minimum in the ball {θ0 +νT u, kuk2 ≤ C } so that the minimumθ̂ satisfies kθ̂ − θ0 k = Op (νT ). Now by a Taylor expansion of the penalised loss function, we obtain
ν2

pen
0 2
0
T
Gpen
T (y; θ0 + νT u) − GT (y; θ0 ) = νT u ∇θ GT (y; θ0 ) + 2 u ∇θθ0 GT (y; θ0 )udP+{p( λTT , |θ0,i + νT ui |) − p( λTT , |θ0,i |)},
i=1

since the third derivative vanishes. We want to proveu0 ∇θ GT (y; θ0 ) + ν2T u0 Hu + ν2T RT (θ0 )dP+ νT−1 {p( λTT , |θ0,i + νT ui |) − p( λTT , |θ0,i |)} ≤ 0) < ,P(∃u, kuk2 = C :

(20)

i=1

where RT (θ0 ) = u0 {∇2θθ0 GT (y; θ0 ) − H}u. First, using the matrix derivatives formula of Abadirand Magnus (2005), we obtain for the score∇θ GT (y; θ0 ) ==T
1 P∇θ `(ys , s ≤ t; θ0 )
T
t=1
T

P
Zm,t−1 ⊗ {xt − Ψ0,1:m Zm,t−1 } .
− T1
t=1

As for the Hessian, we aim at extracting the form tr(L(dΛ)0 M (dΛ)) for L (resp. M ) any squarem × m matrix (resp. p × p). We thus obtain
1∇2θθ0 GT (y; θ0 ) =
T

T
X

0
(Zm,t−1 Zm,t−1⊗ Ip ).

t=1

Under assumption 3, by the Central Limit Theorem of Billingsley (1961)u0 ∇θ GT (y; θ0 ) = Op (n−1/2 u0 Mu),
with M = E[∇θ `(y; θ0 )∇θ0 `(y; θ0 )] assumed well defined by assumption 5. By the ergodic Theoremof Billingsley (1995), we haveP∇θθ0 GT (y; θ0 ) −→ H,
T →∞

33

which implies RT (θ0 ) = op (1). We now focus on the penalty terms. First, we show that thepenalty functions satisfy assumption 6. The Lasso satisfies∀i ∈ A, ∇θi p(

λT
λT
λT
, |θ0,i |) =sgn(θ0,i ), ∇2θi θi p( , |θ0,i |) = 0.
T
T
T

For 0 < q < 1, the Bridge satisfies∀i ∈ A,∇θi p( λTT , |θ0,i |) = λTT q|θ0,i |q−1 sgn(θ0,i ),∇2θi θi p( λTT , |θ0,i |) = λTT q(q − 1)|θ0,i |q−2 sgn(θ0,i ),thus the second order derivative of the Bridge converges to 0 when λT = o(T ). As for the scad,
we have(bscad λTT − |θ0,i |)+λT
λT
1{|θ |≤ λT } +1.
∀i ∈ A, ∇θi p( , |θ0,i |) =λT{|θ0,i |> T }
0,i
T
T
T
(bscad − 1) λTT
As a consequence, the scad penalty is twice continuously differentiable for λTT < |θ0,i |, which
implies that ∀i ∈ A, ∇2θi θi p( λTT , |θ0,i |) = o(1). In the MCP case, we have∀i ∈ A, ∇θi p(|θ0,i |λT
λT
, |θ0,i |) = (−)sgn(θ0,i )1{|θ |≤b λT } .
mcp T
0,i
T
T
bmcp

Under λT = o(T ), we straightforwardly obtain ∇2θi θi p( λTT , |θ0,i |) = o(1) for non-zero components.
Now for any i ∈ A ⊂ {1, · · · , d}, and since the penalties are coordinate-separable, we have
p(


ν2
λT
λT
λT
λT
, |θ0,i +νT ui |)−p( , |θ0,i |) = νT ui sgn(θ0,i )∇θi p( , |θ0,i |)+ T u2i ∇2θi θi p( , |θ0,i |) 1+o(1) .
T
T
T
2
T

Hence, using p( λTT , 0) = 0, we have
|

X

{p(

i∈A1

Using kuk1 ≤


ν2
λT
λT
, |θ0,i + νT ui |) − p( , |θ0,i |)}| ≤ νT kuk1 A1,T + T kuk22 A2,T 1 + o(1) .
TT2pcard(A)kuk2 , and under assumption 6, the third derivative being dominated, we

obtainXp
ν2λTλT| {p( , |θ0,i + νT ui |) − p( , |θ0,i |)}| ≤ νT card(A)kuk2 A1,T + T kuk22 A2,T .
T
T
2
i∈A

34

Then, denoting δT = λmin (H)C2 νT , and using ν2T E[u0 ∇2θθ0 `(y; θ0 )u] ≥ δT , we deduce that (20)can be bounded asP(∃u, kuk2 = C : u0 ∇θ GT (y; θ0 ) + ν2T u0 E[∇2θθ0 `(y; θ0 )]u + ν2T RT (θ0 )dP+νT−1 {p( λTT , |θ0,i + νT ui |) − p( λTT , |θ0,i |)} ≤ 0)
i=1

≤ P(∃u, kuk2 = C : u0 |∇θ GT (y; θ0 )| > δT /6) + P(∃u, kuk2 = C | ν2T RT (θ0 )| > δT /6)dPP(∃u, kuk2 = C : | {p( λTT , |θ0,i + νT ui |) − p( λTT , |θ0,i |)}| > νT δT /6).
i=1

We have for n and C sufficiently large enoughdP{p( λTT , |θ0,i + νT ui |) − p( λTT , |θ0,i |)}| > νT δT /6)
p
ν2≤ P(∃u, kuk2 = C : νT card(A)kuk2 A1,T + 2T kuk22 A2,T > νT δT /6) < /3.
P(∃u, kuk2 = C : |i=1

Moreover, if νT = n−1/2 +pcard(A)A1,T , for C large enoughP(∃u, kuk2 = C : u0 |∇θ GT (y; θ0 )| > δT /6) ≤C2 CstCst≤ 4 < /3,
C
T δT2

where Cst > 0 is a generic constant. Finally, using RT (θ0 ) = op (1), we obtainP(∃u, kuk2 = C : u0 |∇θ GT (y; θ0 )| > δT /6) + P(∃u, kuk2 = C | ν2T RT (θ0 )| > δT /6)
P(∃u, kuk2 = C : |p( λTT , |θ0 + νT u|) − p( λTT , |θ0 |)| > νT δT /6)≤ CCst4 + 2/3 ≤ ,
for T and C large enough. We deduce kθ̂ − θ0 k = Op (νT ).
√Proof of Theorem 3. Let u ∈ Rd such that θ = θ0 + u/ T , and the empirical criterion FT (u) =pen
1/2 (θ̂ − θ ) because θ̂
T {Gpen0
T (y; θ) − GT (y; θ0 )}. Note that FT (u) is minimized at ûT = T

minimizes Gpen
T (Z; θ). Thus ûT = arg min {FT (u)}.
u∈Rd
We first establish the finite distributional convergence of FT (.) to F∞ (.). Then we separate theproof depending on what penalty function we consider. We have the expansion
FT (u) =

d
X
√
√ 0
1
λT
λT
T u ∇θ GT (y; θ0 ) + u0 ∇2θθ0 GT (y; θ0 )u + T{p( , |θ0,i + ui / T |) − p( , |θ0,i |)}.
2
T
T
i=1

35

By assumption 3, by the Central Limit Theorem of Billingsley (1961) and the ergodic Theoremof Billingsley (1995)
√

d

P

T →∞

T →∞

T ∇θ GT (y; θ0 ) −→ NRd (0, M), ∇2θθ0 GT (y; θ0 ) −→ H.

As for the penalty terms, in the MCP and SCAD cases, we proceed as follows. We have
T

d
X

{p(

i=1

√
λT
λT
, |θ0,i + ui / T |) − p( , |θ0,i |)} = ζT (u) + ηT (u),TTwhere using the coordinate-separability property of the penalties
ζT (u) = T

X

{p(

i∈A

X
√
√
λT
λT
λT
, |θ0,i + ui / T |) − p( , |θ0,i |)}, ηT (u) = T{p( , ui / T )}.
T
T
T
c
i∈A

Then we have by a Taylor expansion for the indices i ∈ A
ζT (u) =X√

T ∇θi p(

i∈A

λT
λT
, |θ0,i |)ui sgn(θ0,i ) + ∇2θi θi p( , |θ0,i |)u2i /2(1 + o(1)).
T
T

Under assumption 6, we have A2,T = o(1). We then need to treat the first order term. For
both SCAD and MCP, since their derivatives respectively vanish outside [−bscad λTT , bscad λTT ],[−bmcp λTT , bmcp λTT ] we have∇θi p( λTT , |θ0,i |) = ∇θi p( λTT , |θ0,i |)1{|θ∇θi p( λTT , |θ0,i |)

=

λ

,

T
0,i |≤bscad T }λT∇θi p( T , |θ0,i |)1{|θ |≤b λT } ,
mcp
0,i
T

which implies for any  > 0 and i ∈ A that√> ) ≤ P(|θ0,i | ≤ bscad λTT ) → 0,
P( T ∇θi p( λTT , |θ0,i |)1{|θ |≤b
λT0,i
scad T }√P( T ∇θi p( λTT , |θ0,i |)1{|θ |≤b λT } > ) ≤ P(|θ0,i | ≤ bmcp λTT ) → 0.
0,iAs a consequence,

√

mcp T

T ∇θi p( λTT , |θ0,i |) = op (1). This is a direct consequence of the unbiasedness

property when regularising large coefficients. Thus ζT (u) → 0 as T → ∞. As for i ∈ Ac , we haveηT (u) =X√i∈AcT ∇θi p(



λT
1
λT
, |θi |) θ =0 |ui | + ∇2θi θi ( , |θi |) θ =0 u2i (1 + o(1)),
i
i
T
2
T
36

√
Based on the assumption that lim ∇x p( λTT , x) = λTT and λT = O( T ), we deducex→0+ηT (u) → λ0X|ui |.

i∈Ac

As a consequence, by Lemma 2, where in the latter we take GT (u) = G̃T (u) =√T ∇θ GT (y; θ0 )u+d
1 0 2
2 u ∇θθ0 GT (y; θ0 )u, we obtain arg min{FT } −→ arg min{F∞ }.

uuFor the Bridge estimator, we have for any index i ∈ A ∪ Ac and under the rate λT /T q/2 → λ0 ,
then
√
λT
λT
, |θ0,i + ui / T |) − p( , |θ0,i |)}
T
T
d
d
X
X
√
= λT
{|θ0,i + ui / T |q − |θ0,i |q } → λ0 |ui |q 1θ0,i =0 .

T {p(

i=1

k=1
d

Now we need to prove that arg min{FT } −→ arg min{F∞ } for any u and n sufficiently large. Touu
do so, we use Theorem 7 of Kim and Pollard (1990) (see Appendix B) and show arg min{FT } =uOp (1). We first have the expansiond√√P
FT (u) = GT (y; θ0 + u/ T ) − GT (y; θ0 ) + λT {|θ0,i + ui / T |q − |θ0,i |q }√≥ GT (y; θ0 + u/ T ) − GT (y; θ0 ) − λT

i=1
d
P

√
|ui / T |q

i=1

√

d
√
P
≥ GT (y; θ0 + u/ T ) − GT (y; θ0 ) − λ0 + |ui / T |q := F̃T (u),i=1where, following the proof of Theorem 3 of Knight and Fu (2000),  is such that λ/T q/2 ≤ λ0 + .
√
Then, expanding GT (y; θ0 + u/ T ) − GT (y; θ0 ) in F̃T (u), we have√F̃T (u) =

d
√
X
1 0 2
T u ∇θ GT (θ) + u ∇θθ0 GT (θ)u − λ0 + |ui / T |q .
2
0

i=1

The second term, which is quadratic in u, dominates the term with |ui |q . Hence arg min{F̃T } =uOp (1), which in turns implies arg min{FT } = Op (1). We then obtain for the Bridgeu√d
T (θ̂ − θ0 ) = arg min{FT } −→ arg min{F∞ }.
uu37Finally, for the Lasso estimator, we have√λT
λT
, |θ0,i + ui / T |) − p( , |θ0,i |)}
T
T
d
d
X
X
√

= λT
{|θ0,i + ui / T | − |θ0,i |} → λ0
ui sgn(θ0,i )1θ0,i 6=0 + |ui |1θ0,i =0 .

T {p(

k=1

k=1
d

We thus proved that FT (u) −→ F∞ (u), for a fixed u. Let us observe thatu∗T = arg min {FT (u)},u√and FT (.) admits as a minimizer u∗T = T (θ̂ − θ0 ). As FT is convex and F∞ is continuous, convex
and has a unique minimum, then by the convexity Lemma 1, we obtain√d

T (θ̂ − θ0 ) = arg min{FT } −→ arg min{F∞ }.
uu

0 , θ 0 )0 . To prove the support recovery consistency, we
Proof of Theorem 4. Let us define θ = (θAAcshow with probability tending to one when T → ∞, under kθA − θ0,A k = Op (T −1/2 ) and suitable
regularisation rates depending on the penalty, thatc
Gpen
T (y; θA , 0A ) =To prove (21), for any√minkθAc k≤CT −1/2c{Gpen
T (y; θA , θA )}.

(21)

T -consistent θA , we show that over the set {i ∈ Ac , θi : |θi | ≤ T −1/2 C}

for C > 0
−1/2 C,∇θi Gpen
T (y; θ) > 0 when 0 < θi < Tpen∇θi GT (y; θ) < 0 when −T −1/2 C < θi < 0,(22)

with probability converging to 1. For any index i ∈ Ac , by a Taylor expansion around the true
parameter, we have
λT
, |θi |)sgn(θi )TλT= ∇θi GT (y; θ0 ) + ∇2θi θi GT (y; θ0 ) θi − θ0,i + ∇θi p( , |θi |)sgn(θi ).
T∇θi Gpen
T (y; θ) = ∇θi GT (y; θ) + ∇θi p(

38

By the Central Limit Theorem of Billingsley (1961) and the ergodic Theorem of Billingsley (1995),
we have√T ∇θi GT (y; θ0 ) = Op (1), ∇2θθ0 GT (y; θ0 ) −→ E[∇2θθ0 `(ys , s ≤ t; θ0 )].
P

T →∞We thus obtain for the SCAD and MCP penalties∇θi GT (y; θ) = Op (T −1/2 ) + ∇θi p( λTT , |θi |)sgn(θi ) √= λTT { λTT ∇θi p( λTT , |θi |)sgn(θi ) + Op ( λTT )}.
As a consequence, under the condition lim lim inf λTT ∇x p( λTT , x) > 0 and if the regularisation
T →∞

x→0+

T
parameter satisfies Tλ1/2→ ∞, we deduce that the sign of the gradient entirely depends on thesign of θ̂. This this proves (22).
For the Bridge penalty, following the same reasoning as in the SCAD and MCP, the nonpenalised terms are of order Op (T −1/2 ). As for the penalty, we have∇θi p(

λT
λT
λT, |θi |) =q|θi |q−1 sgn(θi ) = (q+1)/2 q|T 1/2 θi |q−1 sgn(θi ).
T
TT

As a consequence, we obtainq/2

λT{q|T 1/2 θi |q−1 sgn(θi ) + Op ( TλT )}.
∇θi GT (y; θ0 ) = T (q+1)/2

Thus, under the assumption that λT /T q/2 → ∞, this this proves (22).
We now turn to the asymptotic distribution. We prove that θ̂Ac degenerates at 0Ac withprobability approaching one. Now by a Taylor expansion around θ0,i , for i ∈ A, we have
λT
, |θ̂i |)sgn(θ̂i )TXλT= ∇θi GT (y; θ0 ) +∇2θi θj GT (y; θ0 )(θ̂j − θ0,j ) + ∇θi p( , |θ0,i |)sgn(θ0,i )T∇θi GT (y; θ̂) + ∇θi p(j∈A+ ∇2θi θi p(

λT, |θ0,i |)(θ̂i − θ0,i )(1 + o(1)).
T

39

Then inverting this relationship and multiplying by√

T , we obtain in vector form with respect tothe elements in A
∇2AA GT (y; θ0 ) + ST,AA√

T {(θ̂ − θ0 )A + ∇2AA GT (y; θ0 ) + ST,AA

−1√
bT,A } = − T ∇A GT (y; θ0 ),which implies√

T (θ̂ − θ0 )A = − ∇2AA GT (y; θ0 ) + ST,AA√−1 √T ∇A GT (y; θ0 ) − T bT,A ,

where
0
bT,A
= ∇θ1 p( λTT , |θ0,1 |)sgn(θ0,1 ), · · · , ∇θk0 p( λTT , |θ0,k0 |)sgn(θ0,k0 ) ,ST,AA = diag(∇2θi θi p( λTT , |θ0,i |), i = 1, · · · , k0 ).
Now, since λT /T → 0, this implies that A1,T → 0 for the SCAD and MCP. As for the Bridge, the√√T -consistency requires λT = O( T ) by Theorem 2 and the oracle property requires λT T −q/2 →∞. As in the SCAD and MCP cases, A1,T → 0. Thus in all cases, bT,A and ST,AA vanish for T
large enough. We thus deduce that by the central limit theorem for U-statistics and the Slutsky
theorem√d−1
T θ̂ − θ0 A −→ NRk0 (0, H−1
AA MAA HAA ),
T →∞

with
M = E[∇θ `(ys , s ≤ t; θ0 )∇θ0 `(ys , s ≤ t; θ0 )], H = E[∇2θθ0 `(ys , s ≤ t; θ0 )].

Proof of Theorem 5. Under the Theorem’s assumption, θ̂ −→ θ0∗ . Now let us proveP

T →∞∀ > 0, lim P(kγ̂ − γ0 k > ) = 0.
T →∞0First, note that we work under the assumption that β0 = (θ0∗ , γ00 ) is a well-separated point.
Second, let us first establish that for every β̄ ∈ Θ0\2 with kγ̄ − γ0 k > 0, where Θ0\2 = {β : β =40(θ0∗ , γ) ∈ Θ1 × Θ2 } = {θ0∗ } × Θ2 , and every π > 0, there exists an open ball V(β̄, π) around β̄ inthe space Θ0\2 such thatEβ0 [f (ys , s ≤ t; β)] ≥ Eβ0 [f (ys , s ≤ t; β̄)] − π.

inf

(23)

β∈V(β̄,π)To prove this statement, for a given β̄ = (θ̄, γ̄) ∈ Θ0\1 , where Θ0\1 = {β : β = (θ, γ0 ) ∈ Θ}, andγ̄ 6= γ0 , consider a sequence of open balls with radius 1/k, k ∈ N defined by Vk (β0 ) = {β ∈ Θ0\1 :kβ − β0 k ≤ 1/k}. Since the sequence of random variablesinff (ys , s ≤ t; β) k is increasing,β∈Vk (β0 )then by the Beppo-Levi Theorem
lim Eβ0 [

k→∞

inf

f (ys , s ≤ t; β)] = Eβ0 [f (ys , s ≤ t; β0 )],β∈Vk (β0 )which thus proves (23). Now, under the Theorem’s assumption, θ̂ −→ θ0∗ . We thus need to proveP

T →∞∀ > 0, lim P(kγ̂ − γ0 k > ) = 0. Invoking (23), for any given π > 0 and β̄ ∈ Θ0\2 , β̄ 6= β0 , withT →∞kγ̄ − γ0 k ≥ /2, we can find an open ball U(β̄) ⊂ Θ0\2 such thatEβ0 [ inf f (ys , s ≤ t; β)] ≥ Eβ0 [f (ys , s ≤ t; β̄)] − π.
β∈U (β̄)

Since the function γ 7→ Eβ0 [f (ys , s ≤ t; θ0∗ , γ)] − Eβ0 [f (ys , s ≤ t; θ0∗ , γ0 )] defined on Θ0\2 , is strictly
positive because β0 is a well separated point and continuous on the compact subset C0 () = {β ∈Θ0\2 : kγ − γ0 k ≥ /2}, it reaches its minimum 2µ > 0. Hence, for any given β̄ ∈ C0 (), we have
π := π(β̄) = Eβ0 [f (ys , s ≤ t; β̄)] − Eβ0 [f (ys , s ≤ t; β0 θ0∗ , γ0 )] − µ > 0.
Moreover, define U(β0 ) = {β ∈ Θ0\2 : kβ − β0 k < }. Then[Θ0\2 ⊂ U(β0 ) ∪U(β).

β∈C0 ()Since Θ0\2 can be covered by a finite set of open ball by sequential compactness, there is a finite
set of points β1 , · · · , βT in C0 () such thatΘ0\2 ⊂ U(β0 ) ∪n[

U(βi ).

i=1

41

We thus deduceP(kγ̂ − γ0 k > ) ≤ P((θ0∗ , γ̂) ∈n[U(βi )) ≤

i=1

n
X

P((θ0∗ , γ̂) ∈ U(βi )).

i=1

By definition of β̂, for any i = 1, · · · , n,P((θ0∗ , γ̂) ∈ U(βi )) ≤ P( inf LT (y; β) ≤ LT (y; θ0∗ , γ̂))β∈U (βi )≤ P( inf LT (y; β) ≤ LT (y; θ̂, γ̂) + |LT (y; θ0∗ , γ̂) − LT (y; θ̂, γ̂)|)β∈U (βi )≤ P( inf LT (y; β) ≤ LT (y; θ̂, γ0 ) + |LT (y; θ0∗ , γ̂) − LT (y; θ̂, γ̂)|)β∈U (βi )≤ P(Eβ0 [ inf f (ys , s ≤ t; β)] ≤ Eβ0 [f (ys , s ≤ t; β0 )] + |LT (y; β0 ) − Eβ0 [f (ys , s ≤ t; β0 )]|β∈U (βi )+|LT (y; θ0∗ , γ̂) − LT (y; θ̂, γ̂)| + |RT (βi )|),
where RT (βi ) = T1

TP

inf f (ys , s ≤ t; β) − Eβ0 [ inf f (ys , s ≤ t; β)]. Using (23), we have

t=1β∈U (βi )β∈U (βi )Eβ0 [ inf f (ys , s ≤ t; β)] ≥ Eβ0 [f (ys , s ≤ t; β0 )] + µ.
β∈U (βi )Thus
P((θ0∗ , γ̂) ∈ U(βi ))≤ P(µ ≤ |RT (βi )| + |LT (y; β0 ) − Eβ0 [f (ys , s ≤ t; β0 )]| + |LT (y; θ0∗ , γ̂) − LT (y; θ̂, γ̂)|)≤ P(µ/3 ≤ |RT (βi )|) + P(µ/3 ≤ |LT (y; β0 ) − Eβ0 [f (ys , s ≤ t; β0 )]|) + P(µ/3 ≤ |LT (y; θ0∗ , γ̂) − LT (y; θ̂, γ̂)|).
Let us focus on P(µ/3 ≤ |RT (βi )|). Although the quantity f (ys , s ≤ t; β) is not necessarily
integrable, the Ergodic Theorem can still be applied to Eβ0 [f (ys , s ≤ t; β)]. Moreover, f (ys , s ≤t; β) is a measurable function of an ergodic process, by assumption 8, the Ergodic Theorem canbe applied to the process ( inf f (ys , s ≤ t; β)), that isβ∈U (βi )T

1X
inf f (ys , s ≤ t; β) = Eβ0 [ inf f (ys , s ≤ t; β)].
lim inf
T →∞ Tβ∈U (βi )β∈U (βi )t=142Thus, for T > T1 , we haveP(µ/3 ≤ |RT (βi )|) ≤ /3.
By the Ergodic Theorem, for T > T2 , thenP(µ/3 ≤ |LT (y; β0 ) − Eβ0 [f (ys , s ≤ t; β0 )]|) ≤ /3.
Finally, we need to control for the last probability. To do so, we need to prove
sup |LT (y; θ̂, γ) − LT (y; θ0∗ , γ)| = op (1).

γ∈ΘTo do so, by a Taylor expansion of LT (y; θ̂, γ) around θ0∗ , we have
LT (y; θ̂, γ) = LT (y; θ0∗ , γ) + (θ̂ − θ0∗ )0 ∇θ LT (y; θ̃, γ),
where kθ̃ − θ̂k ≤ kθ̂ − θ0∗ k. Using the consistency of θ̂, it is sufficient to prove that
T

1X
sup
k∇θ f (ys , s ≤ t; θ, γ)k = Op (1),Tθ∈Θ:kθ−θ0∗ k≤α
t=1

for a small α > 0. The score (applied here with respect to θ = vec(Ψ1:m )) is given by∇θ f (ys , s ≤ t; θ, γ)1[xt − c∗ − Φxt−1 − Ξ{xt−1 − Ψ1:m Zm,t−2 }]0 [xt − c∗ − Φxt−1 − Ξ{xt−1 − Ψ1:m Zm,t−2 }]2{Zm,t−2 } ⊗ {Ξ0 [xt − c∗ − Φxt−1 − Ξ{xt−1 − Ψ1:m Zm,t−2 }]} .

= ∇θ=As a consequence, there exists some positive constant C > 0 such that, for any α > 0,supθ∈Θ:kθ−θ0∗ k≤αk∇θ f (ys , s ≤ t; θ, γ)k ≤ Csupθ∈Θ:kθ−θ0∗ k≤αkZm,t−2 kkΞkkxt −c∗ −Φxt−1 −Ξ{xt−1 −Ψ1:m Zm,t−2 }k.

Based on xt = c∗ + Φxt−1 + Ξ{xt−1 − Ψ1:m Zm,t−2 } + vt , we obtain
Xt = c∗ + ΛXt−1 + v t ,43

where

 ∗
 

xt
c
vt
Φ + Ξ −ΞΨ1 −ΞΨ2
 xt−1 
0
0
 Ip
0
0


 
 

 xt−2 
0
0
 0
Ip
0
Xt = 
,v =  ,c =  ,Λ = 
 ..  t  ..  .. 
 ..
.
.
 . 
.
.
 .
.
···
xt−m
0
0
0
···
···



· · · −ΞΨm
···
0 

···
0 
.
.. 
···
. 
Ip
0

Assumption 9 provides a sufficient condition for stationarity, that is kΛks < 1, with k.ks the
spectral radius of Λ. Hence, under kΛks < 1,
T

1X
sup
k∇θ f (ys , s ≤ t; θ, γ)k = Op (1).
Tθ∈Θ:kθ−θ0∗ k≤α
t=1

As a consequence, for T > T3 , we deduceP(µ/3 ≤ |LT (y; θ0∗ , γ̂) − LT (y; θ̂, γ̂)|) < /3.
We can conclude that for T > T1 ∨ T2 ∨ T3 , we haveP((θ0∗ , γ̂) ∈ U(βi )) < .
This thus proves the desired consistency result.
Proof of Theorem 6. Through a Taylor expansion around β0 , we obtain for the γ̂ component
0 = ∇γ LT (y; θ̂), γ̂) = ∇γ LT (y; β0 ) + ∇θ LT (y; β̄)A (θ̂ − θ0 )A + ∇γ LT (y; β̄)(γ̂ − γ0 ),where kβ̄ − β0 k < kβ̂ − β0 k. Then inverting this relationship, multiplying by√

T and using theasymptotic expansion of the first step estimator, we obtain√

T (γ̂ − γ0 )√√= (−∇γ LT (y; β̄))−1 ∇θ LT (y; β̄)A T (θ̂ − θ0 )A + (−∇γ LT (y; β̄))−1 T ∇γ LT (y; β0 ) √= (−∇γ LT (y; β̄))−1 ∇θ LT (y; β̄)A {−bT,A + (− ∇2AA GT (y; θ0 ) + ST,AA )−1 T ∇A GT (y; θ0 )}√+(−∇γ LT (y; β̄))−1 T ∇γ LT (y; β0 ).
44Hence, by Slutsky’s Theorem and the Central Limit Theorem, we obtain the desired asymptotic
distribution since ST,AA = 0 and bT,A = 0 for T sufficiently large.

D

Some competing M-GARCH models

The BEKK model directly generates a variance-covariance process. Developed by Baba, Engle,
Kraft and Kroner, in a preliminary version of Engle and Kroner (1995), the BEKK is specified fora p-dimensional random vector t as
1/2
 t = Ht ηt , with Ht := E[t 0t |Ft−1 ]  0 so that
q P
K
r P
K
P
P
0 ,

Akj t−k 0t−k A0kj +Bij Ht−i Bij
 Ht = Ω +k=1j=1i=1i=1

where K is an integer, Ω, Akj and Bkj are square p × p matrices and Ω  0. One advantageof the BEKK model is there is no positive semi-definite constraint on the Akj and Bkj matrices.
However, it imposes highly artificial constraints on the volatilities and covariances of the components. As a consequence, the coefficients of a BEKK representation are difficult to interpret.
In our application, a scalar BEKK was considered, where Akj and Bkj are scalar with K = 1,
q = r = 1, together with a Gaussian QMLE estimation.
Beside BEKK type dynamics, factor models provide rather natural alternatives. The OGARCH assumes the decomposition Ht = P Λt P 0 , where Λt = diag(λ1,t , · · · , λK,t ), with K thenumber of factors. Here, we choose K = p factors and each λt is supposed to follow a univariate
GARCH(1,1) process that is estimated by maximum likelihood. The matrix P is nonsingular and
it is estimated by PCA on the empirical variance-covariance matrix of t : see Alexander (2001),
e.g.
Beside the latter direct specification of the covariance matrices (Ht ) dynamics, an alternativeroad is to split the task into two parts: individual volatility dynamics on one side, and correlation45

dynamics on the other side. The most commonly used correlation process is the Dynamic Conditional Correlation (DCC) of Engle (2002). In its BEKK form, the general DCC model is specified
as


1/2

t = Ht ηt , with Ht := E[t 0t |Ft−1 ]  0 so that



?−1/2
?−1/2
Ht = Dt Rt Dt , Rt = Qt
Qt Qt
,
q
r
P
P



Mk Qt−k Mk0 + Wl ut−l u0t−l Wl0 , Qt = Ω +

(24)

l=1

k=1


p
p
p
p
where Dt = diagh11,t , h22,t , . . . , hpp,t , ut = (u1,t , . . . , up,t ) with ui,t = i,t / hii,t , Qt =[qij,t ], Q?t = diag (q11,t , q22,t , . . . , qpp,t ). The model is parameterized by some deterministic matrices(Mk )k=1,··· ,q , (Wl )l=1,··· ,r and a positive definite p × p matrix Ω. Alternatively, Engle (2002)
considered a VEC-type specification too. Denoting by

the Hadamard matrix product, the(Qt )-dynamics becomeqX
Qt = Ω +Bk∗r
X
Qt−k +
Al

k=1

ut−l u0t−l ,

(25)

l=1

where the deterministic matrices (Bk )k=1,··· ,q and (Al )l=1,··· ,r must be positive semi-definite.
Since the number of parameters of the latter models is of order O(p2 ), the matrices Mk
and Wl (resp. Bk ’s and Al ) are often assumed to be scalar. This is typically a strong andquestionable constraint, particularly when the dimension p increases or when the variables in (t )are heterogeneous. Furthermore, their inference is usually carried out trough the QML method,
based on a Gaussian or Student quasi likelihood function. Under this methodology, applying a
regularisation method, even possible, is numerically arduous and no general asymptotic results
exist in this case (to the best of our knowledge), due to the non-convexity of the QML criterion.
If Rt = R a constant correlation matrix, then (24) becomes the Constant Conditional Correlation (CCC) model.

46