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.
59,726 characters · 11 sections · 107 citation commands
Multivariate Fractional Components Analysis
\thispagestyle{empty} \setcounter{page}{0}
\paragraph{\bf Abstract.}
We propose a setup for fractionally cointegrated time series which is formulated in terms of latent integrated and short-memory components. It accommodates nonstationary processes with different fractional orders and cointegration of different strengths and is applicable in high-dimensional settings. In an application to realized covariance matrices, we find that orthogonal short- and long-memory components provide a reasonable fit and competitive out-of-sample performance compared to several competing methods.
\paragraph{\bf Keywords.}
Long memory, fractional cointegration, state space, unobserved components, factor model, realized covariance matrix.
\paragraph{\bf JEL-Classification.}
C32, C51, C53, C58.
Multivariate fractional integration and cointegration models have proven valuable in a wide range of empirical applications from macroeconomics and finance. They generalize the standard concept of cointegration by allowing for non-integer orders of integration both for the observations and for equilibrium errors; see GilHua2008 for a literature review. In the field of macroeconomics, such models have turned out to be relevant in analyses of purchasing power parity beginning with CheLai1993, of the relation between unemployment and input prices CapGil2002 and of broader models for economic fluctuations Mor2006. The empirical finance literature has considered fractional cointegration, e.g., for analysing international bond returns DueSta1998, for modeling co-movements of stock return volatilities BelMor2006, for assessing the link between realized and implied volatility Nie2007 and for quantifying risk in strategic asset allocation problems SchTscBu2008. From a methodological point of view, semiparametric techniques for inference on the cointegration rank, the cointegration space and memory parameters have been very popular among empirical researchers, although the development of optimal parametric inferential methods for models with triangular or fractional vector error correction representations has recently made considerable progress RobHua2003,AvaVel2009,Las2010,JohNie2012.
Despite their flexibility and their computationally simple treatment, semiparametric models are limited in scope since they aim to describe low-frequency properties only and are hence not appropriate for impulse response analysis and forecasting. While semiparametric techniques have been developed to cope with multivariate processes of different integration orders and multiple fractional cointegration relations of different strenghts CheHur2006,HuaRob2010,Hua2009, there seems to be a lack of parametric models of such generality. Furthermore, the usual error correction and triangular models with their abundant parametrization are not deemed appropriate for time series of dimension, say, larger than five.
In this paper, we propose new models for multivariate fractionally integrated and cointegrated time series which are formulated in terms of latent purely fractional and additive short-memory components. With a “type II” definition of fractional integration Rob2005, this approach allows for a flexible modeling of possibly nonstationary time series of different fractional integration orders. It permits cointegration relations of different strengths as well as polynomial cointegration GraLee1989, i.e., cointegration between the levels of some time series and their (fractional) differences, and guarantees a clear representation of the long-run characteristics. Consequently, our model is among the most general setups regarding its integration and cointegration properties, compared to popular existing models for cointegrated processes. The unobserved components formulation benefits the modeling of relatively high-dimensional time series. For this situation we propose a parsimonious parametrization based on dimension reduction and dynamic orthogonal components in the spirit of PanYao2008 and MatTsa2011.
In contrast to our parametric approach, latent fractional components have mostly been studied in semiparametric frameworks. RayTsa2000 use semiparametric memory estimators and canonical correlations to infer the existence of common fractional components, Mor2004 proposes a frequency domain principal components estimator, Mor2007 estimates components of a single fractional integration order by univariate permanent-transitory (or persistent-transitory) decompositions followed by a principal component analysis of the permanent (or persistent) components and LucVer2015 estimate their fractional factor model by fitting long-memory models to the principal components of a large panel of time series. In a setup closest to ours, CheHur2006 suggest a semiparametric frequency domain methodology to identify and estimate cointegration subspaces which annihilate fractional components of different memory.
Recent parametric frameworks competing to ours have either been much more restrictive, or have a different focus, e.g., on numerical simulation-based estimation methods, or on panel data analysis. Using a Bayesian approach, HsuRayBr1998 discuss a bivariate process sharing one stationary long-memory component, while MesKooOo2016 consider simulated maximum likelihood estimation of models with one or more latent stationary ARFIMA components. On the other hand, ErgVel2017 as well as Erg2017 focus on the elimination of common fractional components and are motivated as alternatives to prior unit root testing. Contrary to our approach, they eliminate the common factor structure, which they treat as nuisance.
As the second main contribution of this paper, our model is applied to forecasting daily realized covariance matrices. In this setup, the strengths of our approach become apparent. In realized covariance modelling, typically high-dimensional processes with strong persistence and a pronounced co-movement in the low-frequency dynamics are considered. In time series of log variances and z-transformed correlations for six US stocks, we find that common orthogonal short- and long-memory components with two different fractional integration orders provide a reasonable fit. Since the dimension of the dataset is reduced to a smaller number of latent processes, our model becomes a factor model. A pseudo out-of-sample study shows that the fractional components model provides a superior forecasting accuracy compared to several competitor methods. In addition to the favorable forecast properties, our methods can be applied to study the cointegration properties of stock market volatilities. These are of particular importance for longer-term portfolio hedging and the analysis of systematic risk.
The paper is organized as follows. Section (ref) introduces the general setup and clarifies its integration and cointegration properties. In section (ref), its relation to existing models for multivariate integrated time series is discussed. In section (ref), a specific model appropriate for relatively high-dimensional processes is considered. The empirical application to realized covariance matrices and a pseudo out-of-sample assessment are contained in section (ref) before section (ref) concludes.
We consider a linear model for a $p$-dimensional observed time series $y_t$, which we label a fractional components (FC) setup,
The model is formulated in terms of the latent processes $x_t$ and $u_t$ where $\varLambda$ will always be assumed to have full column rank and the components of the $s$-dimensional $x_t$ are fractionally integrated noise according to
In principle, $s>p$ is possible, but we only consider cases where $s \leq p$ here. For a generic scalar $d$, the fractional difference operator is defined by
where $L$ denotes the lag or backshift operator, $Lx_t = x_{t-1}$. We adapt a nonstationary type II solution of these processes Rob2005 and hence treat $d_j\geq0.5$ alongside the asymptotically stationary case $d_j<0.5$ in a continuous setup, by setting starting values to zero, $x_{jt}=0$ for $t\leq 0$. Nonzero initial values have been considered for observed fractional processes by JohNie2012, but are not straightforwardly handled for our unobserved processes. The solution is based on the truncated operator $\Delta_+^{-d_j}$ Joh2008 and given by \[ x_{jt} = \Delta_+^{-d_j} \xi_{jt} = \sum_{i=0}^{t-1} \pi_i(-d_j) \xi_{j,t-i}, \qquad j=1,\ldots,s. \] Without loss of generality let the components be arranged such that $d_1 \geq \ldots \geq d_s$.
We assume $d_j>0$ for all $j$ in what follows, so that $x_t$ governs the long-term characteristics of the observations $y_t$. These are complemented by additive short-run dynamics which we describe by stationary vector ARMA specifications for $u_t$ in the general case. This ARMA process is given by
where $\varPhi(L)$ and $\varTheta(L)$ are a stable vector autoregressive polynomial and an invertible moving average polynomial, respectively. The disturbances $\xi_t$ and $e_t$ jointly follow a Gaussian white noise ($N\!I\!D$) sequence such that
where at this stage, before turning to identified and empirically relevant model specifications below, we do not consider restrictions on the joint covariance matrix, but only require $\varSigma_{\xi}$ to have strictly positive entries on the main diagonal.
Some remarks regarding the general FC setup are in order. The model as given in (ref) is not identified without further restrictions on the loading matrix $\varLambda$, on the vector ARMA coefficients and on the noise covariance matrix. While restrictions on $\varSigma_{\xi}$ and $\varLambda$ may be based on results in dynamic factor analysis as will be seen below, choosing specific parametrizations for $u_t$ will depend on characteristics of the data and on the purpose of the empirical analysis. Identified vector ARMA structures like the echelon form Luet2005 can be used for a rich parametrization, while a multivariate structural time series approach as described in Har1991 integrates nicely with the unobserved components framework considered in this paper and allows for more restricted parameterization, e.g., by individual or common stochastic cycle components. Below, we introduce a parsimonious model well-suited to relatively high dimensions which is conceptually based on dimension reduction and orthogonal components.
For a characterization of the integration and cointegration properties of our model, we adapt the definitions of these concepts from HuaRob2010, which prove useful here. Hence, a generic scalar process $\rho_t$ is called integrated of order $\delta$ or $I(\delta)$ if it can be written as $\rho_t = \sum_{i=1}^l \Delta_+^{-\delta_i}\nu_{it}$, where $\delta = \max_{i=1,\ldots,l} \{ \delta_i\}$ and $\nu_t=(\nu_{1t},\ldots,\nu_{lt})'$ is a finite-dimensional covariance stationary process with spectral density matrix which is continuous and nonsingular at all frequencies. A vector process $\tau_t$ is called $I(\delta)$ if $\delta$ is the maximum integration order of its components. We call the process $\tau_t$ cointegrated if there exists a nonzero vector $\beta$ such that $\beta'\tau_t$ is $I(\gamma)$ where $\delta-\gamma>0$ will be referred to as the strength of the cointegration relation. The number of linearly independent cointegration relations with possibly differing $\gamma$ is called cointegration rank of $\tau_t$.
By these definitions, $x_{jt}$ is clearly $I(d_j)$ while both $x_t$ and $y_t$ are integrated of order $d_1$. We observe at least two different integration orders in the individual series of $y_t$ whenever $\varLambda_{i1}=0$ for some $i$ and $d_1>d_2$. More generally, $y_{it} \sim I(d_j)$, if $\varLambda_{i1} = \ldots = \varLambda_{i,j-1}=0$ but $\varLambda_{ij}\neq 0$.
To state the cointegration properties of the FC setup (ref), we assume that $s \leq p$, so that all fractional components are reflected by the integration and cointegration structure of $y_t$ and that $\varSigma_{\xi}$ is nonsingular. It is useful to identify all $q$ groups of $x_{jt}$ with identical integration orders and denote their respective sizes by $s_1$, \ldots, $s_q$, such that $d_{s_1+\ldots+s_{j-1}+1}= \ldots = d_{s_1+\ldots+s_j}$ and $s=\sum_{j=1}^q s_j$. Of course, if $q=s$, then $s_1 = \ldots = s_q=1$ and all components of $x_t$ have mutually different integration orders, while for $q=1$ it holds that $s=s_1$ and we observe $d_1=\ldots=d_s$.
To keep notation simple, for a generic matrix $A$ for which a specific grouping of rows and columns is clear from the context, we denote by $A^{(i,j)}$ the block from intersecting the $i$-th group of rows with the $j$-th group of columns. A stacking of several groups of rows $i,\ldots,j$ and columns $k,\ldots,l$ is indicated by $A^{(i:j,k:l)}$. For a grouping in only one dimension we write $A^{(i)}$ or $A^{(i:j)}$, where it shall be clear from the context whether a grouping of rows or columns is considered. Furthermore, we denote the column space of a generic $k \times l$ matrix $A$ by $sp(A)\subseteq \mathds{R}^k$ and its orthogonal complement by $sp^{\bot}(A)$. Further, for $k > l$, the $k \times (k-l)$ orthogonal complement of $A$ will be denoted by $A_{\bot}$, which spans the $(k-l)$-dimensional space $sp^{\bot}(A)$.
According to the grouping of equal individual integration orders in $x_t$, we may therefore rewrite the FC process (ref) as \[ y_t = \varLambda^{(1)} x_t^{(1)} + \ldots + \varLambda^{(q)} x_t^{(q)} + u_t. \] Here, $\varLambda^{(j)}$ is a $p \times s_{j}$ submatrix of $\varLambda$ consisting of columns $\varLambda_{\cdot i}$ for which $s_1+\ldots+s_{j-1}<i\leq s_1+\ldots+s_{j}$, and $x_t^{(j)}$ is a $s_j$-dimensional subprocess of $x_t$ corresponding to components with memory parameter $d^{(j)}:=d_{s_1+\ldots+s_{j-1}+1}=\ldots = d_{s_1+\ldots+s_j}$. Whenever $s_1<p$, there exist $p-s_1$ linearly independent linear combinations $\beta_i'y_t \sim I(\gamma_i)$ and $\gamma_i<d_1$, so that fractional cointegration occurs. Due to our definition of cointegration, this may be a trivial case where a single component $y_{it}$ with integration order smaller than $d_1$ is selected. Since \[ \varLambda^{(1)'}_{\bot} y_t = \varLambda^{(1)'}_{\bot}\varLambda^{(2)} x_t^{(2)} + \ldots + \varLambda^{(1)'}_{\bot}\varLambda^{(q)} x_t^{(q)} +\varLambda^{(1)'}_{\bot} u_t \] is integrated of order $d^{(2)}$, the columns of $\varLambda^{(1)}_{\bot}$ qualify as cointegration vectors and ${\cal S}^{(1)}:=sp^{\bot}(\varLambda^{(1)})$ is the $(p-s_1)$-dimensional cointegration space of $y_t$.
Whenever $s_1+s_2<p$, there are subspaces of ${\cal S}^{(1)}$ forcing a stronger reduction in integration orders. More generally, it holds that $\varLambda_{\bot}^{(1:j)'}y_t \sim I(d^{(j+1)})$ whenever $\sum_{i=1}^js_i<p$ and where we set $d^{(j+1)}=0$ for $j>s$. Analogously to HuaRob2010, for $s=p$ and $j=1,\ldots,q-1$, we call ${\cal S}^{(j)}:=sp^{\bot}(\varLambda^{(1:j)})$ the $j$-th cointegration subspace of $y_t$, for which ${\cal S}^{(q-1)}\subset \ldots \subset {\cal S}^{(1)}$. For $p>s$, ${\cal S}^{(q)}\subset{\cal S}^{(q-1)}$ is a further such subspace. Cointegration vectors in ${\cal S}^{(q)}$ cancel all fractional components and hence reduce the integration order from $d_1$ to zero, the strongest reduction possible in our setup.
Besides this general pattern of cointegration relations, our model features an interesting special case with so-called polynomial cointegration, that is, cointegration relations where lagged observations nontrivially enter a cointegration relation. To see this possibility, consider a bivariate example similar to GraLee1989, where $q=p=2$ and $\xi_{1t}=\xi_{2t}$, so that $\varSigma_{\xi}$ is singular and $x_{2t} = \Delta^{d_1-d_2} x_{1t}$. Augmenting the variables by a fractional difference as $\tilde{y}_t:=(y_{1t}, y_{2t}, \Delta^{d_1-d_2}y_{2t})'$, we obtain a three-dimensional system where levels of $y_t$ enter a nontrivial cointegration relation with a fractional difference to achieve a reduction in integration order from $d_1$ to $\max\{2d_2-d_1,0\}<d_2$. Hence, our setup complements the model of Joh2008, which was the first to handle polynomial cointegration in a fractional setup, and the results of CarSan2018, who derive a Granger representation for the fractional VECM of Gra1986 under polynomial cointegration.
In this section, we clarify the relation of the fractional components model (ref) to popular existing representations for cointegrated processes and show how our model can be represented in alternative ways brought forward in the literature. While our model is among the most general setups with respect to its integration and cointegration properties, the additive modeling of short-run dynamics is new to the literature and gives rise to distinct parametrizations not possible within other representations in a similarly convenient way.
\paragraph{Error correction models.} The most popular representation of cointegrated systems in the $I(1)$ setting is the vector error correction form. Since an early mention by Gra1986, in the fractionally integrated case, e.g., AvaVel2009, Las2010 and JohNie2012 have recently considered such models. In terms of the integration and cointegration properties, the fractional error correction setups are typically restricted to the special case with $q=2$ and $s=p$, such that the observed variables are integrated of order $d^{(1)}$ and there exist $p-s_1$ cointegration relations with errors of order $d^{(2)}$.
Defining the fractional lag operator $L_b:=1-\Delta^b$ Joh2008, we are able to derive the error correction representation for this special case of our model; see appendix (ref). It is given by
where we find $\alpha \beta'=-\varLambda^{(2)}(\varLambda^{(1)'}_{\bot} \varLambda^{(2)})^{-1}\varLambda^{(1)'}_{\bot}$ to precede the error correction term, while \[ \kappa_t:=M(\varLambda^{(1)}\xi_t^{(1)} + \Delta^{d^{(1)}} u_t )-\alpha \beta'(\varLambda^{(2)}\xi_t^{(2)} + \Delta^{d^{(2)}} u_t) \] is integrated of order zero and $M$ is defined in (ref).
The model differs both from the models of AvaVel2009 and from the representation of Joh2008 in the way short-run dynamics are modeled. The literature has considered (fractional) lags of differenced variables and possibly of error correction terms in the VECM representation. Our setup, in contrast, generates autocorrelated $\kappa_t$ by filtering the latent $u_t$ with fractional difference operators. Hence, adding lags of $\Delta^{d^{(1)}} y_t$ in the model (ref) is only an approximate solution and achieving a desired approximation quality may require estimating a large number of parameters.
As we have discussed above, Joh2008 proposes a polynomially cointegrated generalization of his $\mathrm{VAR}_{d,b}$ model which allows terms integrated of orders $d$, $d-b$ and $d-2b$ in the Granger representation Joh2008. Even compared to that specification, our model allows for more general patterns of integration orders and cointegration strengths, since we only assume $d_j > 0$ for all $j$. More in line with the generality envisaged in this paper, TscWebWe2013a present a model with error correction term and different integration orders, while LaVel2014 sequentially fit error correction models to test for cointegration relations of possibly different strengths.
\paragraph{Vector ARFIMA.} An interesting special case of (ref) occurs for $s=p$ and $\varLambda=I$, where each series in $y_{it}$ is driven by a single fractional component and $y_{it} \sim I(d_i)$. This resembles standard vector ARFIMA models with possibly different integration orders; see, e.g., Lob1997 who labels the popularly termed vector ARFIMA class considered here as “model A”. A frequently used submodel is the fractionally integrated vector autoregressive model discussed by Nie2004b. The main difference to these approaches is our additive modeling of short-run dynamics, whereas in the vector ARFIMA setup weakly dependent vector ARMA instead of white noise processes are passed through the fractional integration filters.
Our model belongs to the class of vector ARFIMA processes for integer $d_j\in \{1,2,\ldots\}$, but not for general fractional integration orders. For the case of integer $d_j$, note that $(x_t',u_t')'$ is a finite-order vector ARMA process, and hence $y_t$ as a linear combination is itself in the ARMA class; see Luet1984. For general vector ARFIMA processes, a similar conclusion does not hold. To see this, consider a stylized univariate case of our model with $p=s=1$, where $\Delta^d x_t =\xi_t$ and $(1- \phi L) u_t = e_t$. First note that $(\Delta^d x_t, u_t)'$ has an ARMA structure, and hence $(x_t, u_t)'$ is a vector ARFIMA process. Expanding $(1-\phi L) \Delta^d x_t = (1-\phi L)\xi_t$ and $(1-\phi L) \Delta^d u_t = \Delta^de_t$, we can write the sum, belonging to the fractional components model class, as
The right hand side of this expression is not a finite-order MA process in general, as it has nonzero autocorrelations for all lags, and hence, the process does not belong to the ARFIMA class for non-integer $d$.
\paragraph{Triangular representations.} The models discussed so far have restricted integration or cointegration properties as compared to our model. Even in the most general setup of Joh2008, the integration orders are restricted to be $d$, $d-b$, $d-2b$ for polynomial cointegration. In contrast, Hua2009 and HuaRob2010 have proposed a very flexible model which adapts the triangular form of Phi1991 and its generalization to processes with multiple unit roots StoWat1993 to the fractional cointegration setup.
To derive the triangular representation for our model, we assume that the variables in $y_t$ are ordered in a way that $\varLambda^{(1:j,1:j)}$ is nonsingular for $j=1,\ldots,q$ and restrict attention to the case $s=p$ for notational convenience. The variables are partitioned according to the groups of different integration orders in $x_t$ as $y_t^{(j)}:=(y_{s_1+\ldots+s_{j-1}+1}, \ldots, y_{s_1+\ldots+s_j})'$, $j=1,\ldots,q$. The first block in the triangular system is
where $\omega_t^{(1)}$ is integrated of order zero. The general expression for the $j$-th block of the triangular system is derived in appendix (ref) for $j=2,\ldots,q$, and given by
where also $\omega_t^{(j)}$ is integrated of order zero for $j=2,\ldots,q$. By inverting the fractional difference operators we obtain
where $B$ has a block triangular structure such that $B^{(i,i)}=I$ and $B^{(i,j)}=0$ for $i<j$. A re-ordering of the variables in $y_t$ yields the representation of HuaRob2010.
This representation allows for a semiparametric cointegration analysis of our model using the methods of Hua2009 and HuaRob2010. However, our model differs significantly from straightforward parametrizations of the triangular system, e.g., from assuming a vector ARMA process for $\omega_t$, since in our setup $\omega_t$ as stated in (ref) generally contains fractional differences that cannot be represented within the ARMA framework.
\paragraph{State space approaches.}
BauWag2012 have presented a state space canonical form for multiple frequency unit root processes of different (integer-valued) integration orders. Their discussion is based on unit root vector ARMA models which are separated in pure unit root structures and short-term dynamics. Although the analogy to our model is striking, there are notable differences between their unit root and our fractional setup. Firstly, as discussed in the paragraph on vector ARFIMA models (see (ref)), the fractional components setup (ref) is not nested within a general class comparable to the vector ARMA models, which form the basis of the discussion in BauWag2012. Secondly, in their setting, the introduction of different integration orders is achieved by repeated summation of lower order integrated processes which themselves enter the observations to achieve polynomial cointegration. This is in contrast to the continuous treatment of integration orders in our (type II) fractional setup.
However, fractional components models could be constructed to straightforwardly extend the setup of BauWag2012. Using the fractional lag operator $L_b=1-\Delta^b$ instead of $L$ in the short-run dynamic specification (ref), a stable vector ARMA$_b$ process can be defined by $\tilde{\varPhi}(L_b)\tilde{u}_t = \tilde{\varTheta}(L_b) e_t$ under suitable stability conditions Joh2008. Then, replacing $u_t$ by $\tilde{u}_t$ in the model setup (ref) with $d_j$ restricted to some multiple of $b$ ($d_j=i_jb$, $i_j \in \{1,2,\ldots\}$), the process $y_t$ is in the class of vector ARMA$_b$ models itself, while unit roots in the vector autoregressive polynomial generate the fractional $I(d_j)$ processes. Such a framework could be treated analogously to BauWag2012, but the restriction that all integration orders are multiples of $b$ makes such a framework somewhat less flexible than ours.
So far, we have considered a general modeling setup and discussed its integration and cointegration properties as well as its relation to existing approaches in the literature. We now turn to the discussion of a specific model from this class which bears potential for parsimonious modeling of long- and short-run dynamics in relatively high-dimensional applications. Besides its general interest, this will be the workhorse specification for the empirical application to realized covariance modeling in section (ref).
To introduce the model and emphasize its restrictions as compared to (ref), we decompose the short-term dependent process $u_t$ into an autocorrelated component, $\varGamma z_t$, where $z_t$ is a vector of $s_0$ mutually uncorrelated components with $s+s_0 \leq p$, and a Gaussian white noise component $\varepsilon_t$, respectively. We label the result the dynamic orthogonal fractional components (DOFC) model,
where $x_t$ is generated by purely fractional processes (ref) as above, while \[ (1-\phi_{j1}L -\ldots - \phi_{jk}L^k) z_{jt} = \zeta_{jt}, \qquad j=1,\ldots,s_0, \] are $s_0$ univariate stationary autoregressive processes of order $k$. Regarding the noise processes $\xi_t$, $\zeta_t$ and $\varepsilon_t$, we assume mutual independence over leads and lags, \[ \xi_t \sim N\!I\!D (0, I),\qquad \zeta_t \sim N\!I\!D (0, I)\qquad \text{and} \quad \varepsilon_t \sim N\!I\!D (0, H), \] where $H$ is diagonal with entries $h_{i}>0$, $i=1,\ldots,p$. Note that for $s + s_0 < p$ the DOFC model is a factor model as it allows for dimension reduction.
The model as specified in (ref) and below is not identifiable without further information. Considering $\tilde{y}_t:= \Delta^{d^{(1)}} y_t$ instead of $y_t$ to meet the assumptions of HeaSol2004, their theorem 4 suggests that groups of common components $\Delta^{d^{(1)}}x_t^{(1)}$, \ldots, $\Delta^{d^{(1)}}x_t^{(q)}$, $\Delta^{d^{(1)}}z_t$ can be disentangled (up to rotations within these groups) through their different shapes in spectral densities whenever $d^{(1)}>\ldots>d^{(q)}>0$. Still, there exist observationally equivalent structures with $\tilde{\varLambda}^{(j)} = \varLambda^{(j)} M^{-1}$ and $\tilde{x}_t^{(j)}=M x_t^{(j)}$ which satisfy the model restrictions for orthonormal $M$. Hence, we impose further restrictions on the loading matrices. As is standard practice in dynamic factor analysis, we set the upper triangular elements to zero such that $\varLambda_{rl}^{(j)}=0$ for $r<l$, $j=1,\ldots,q$, and $\varGamma_{rl}=0$ for $r<l$. Certain observables are thus assumed not to be influenced by certain factors.
The model (ref) is very parsimonious considering that it includes both a rich fractional structure as well as short-run dynamics with co-dependence. This is possible by comprising three components of parsimony which have been brought forward in the statistical time series literature. Firstly, there are $p-s-s_0\geq 0$ white noise linear combinations of $y_t$. A strict inequality implies a reduced dimension in the dynamics of $y_t$ which is characteristic for so-called statistical factor models; see PanYao2008, LamYaoBa2011 and LamYao2012. In contrast, the model (ref) does not belong to this class in general, since it allows for $s \geq p$ and general forms of autocorrelation in $u_t$. Secondly, all cross-sectional correlation stems from the common components which is a familiar feature from classical factor analysis AndRub1956. Thirdly, both the fractional and the nonfractional components are mutually orthogonal for all leads and lags.
Combined with semiparametric techniques of fractional integration and cointegration analysis, existing methods for statistical factor and dynamic orthogonal components analysis MatTsa2011 can be used to justify the model assumptions and may be useful in the course of model specification. For final model inference, maximum likelihood estimation based on a state space representation is the preferred method. Both steps will be illustrated in the empirical application of the next section.
We apply the fractional components approach to the modeling and forecasting of multivariate realized stock market volatility which has recently received considerable interest in the financial econometrics literature.
We use the dataset of ChiVoe2011 which comprises realized variances and covariances from six US stocks, namely (1) American Express Inc., (2) Citigroup, (3) General Electric, (4) Home Depot Inc., (5) International Business Machines and (6) JPMorgan Chase & Co for the period from 2000-01-01 to 2008-07-30 ($n=2156$). The data are available from \url{http://qed.econ.queensu.ca/jae/2011-v26.6/chiriac-voev}.
Different transformations of the realized covariance matrices have been applied to fit dynamic models to data of this kind. Wei2014 discusses these transforms and considers a general framework nesting several previously applied approaches. His results suggest that applying linear models to a multivariate time series of log realized variances along with z-transformed realized correlations is a reasonable choice in practice. We follow this approach and base our empirical study on the 21-dimensional time series
where $X_t$ is the $6 \times 6$ realized covariance matrix at period $t$, and the z-transforms are \[ Z_{ij,t} = 0.5[\log(1+R_{ij,t})-\log(1-R_{ij,t})], \quad R_{ij,t}=\frac{X_{ij,t}}{\sqrt{X_{ii,t}X_{jj,t}}}. \] All time series (grey) of log variances and their maxima and minima for a given day $t$ (black) are depicted in figure (ref), while z-transformed correlations are shown in figure (ref).
Recent approaches to modeling realized covariance matrices have successfully used long-memory specifications ChiVoe2011, or found co-movements between the processes well-re\-pre\-sen\-ted by dynamic factor structures; see BauVor2011 and Gri2013. In the related problem of forecasting univariate realized variances, factor models with long-memory dynamics have already been proposed. While BelMor2006 use frequency-domain principal components techniques to assess the low-frequency co-movements, LucVer2015 apply time-domain principal components to their high-dimensional series and apply fractional integration techniques to both estimated factors and idiosyncratic components. Recently, AsaMcA2014 have considered long-memory factor dynamics also for the modeling of realized covariance matrices, where again a semiparametric factor approach precedes a long-memory analysis in their two-step approach.
Our fractional components model DOFC (ref), applied to the time series (ref), offers various advantages to researchers and practitioners in the field. (a) Our methods offer new insights in the integration and cointegration properties of stock market volatilities, for which fractional components structures of different integration orders have not been investigated so far. (b) Fractional cointegration between variances and correlations is of particular interest for the understanding of longer-term portfolio hedging and systemic risk assessment, but has not found attention in the existing literature. (c) Our state space approach for variances and correlations also features other relevant aspects of volatility modeling. It offers a separation into short-term and long-term components in the spirit of EngLee1999, directly accounts for measurement noise, and is applicable in datasets of higher dimensions. The parameter-driven state space approach our specification enables yields (d) practicability in case of missing values, while it (e) straightforwardly carries over to stochastic volatility frameworks for daily return data in the spirit of HarRuiSh1994.
We investigate whether the constraints imposed in the DOFC model (ref) are reasonable for the dataset under investigation. Semiparametric methods are used to assess these restrictions and to obtain reasonable starting values for the parametric estimation of our model.
The model (ref) implies that there are $s+s_0$ components which govern the dynamics of $y_t$, and hence, for $p>s+s_0$, there is a dimension reduction in terms of the autocorrelation characteristics. PanYao2008 study time series with such properties and propose a sequential test to infer the dynamic dimension of the process, allowing for nonstationarity of the autocorrelated components. The algorithm sequentially finds the least serially correlated linear combinations of $y_t$, subsequently testing the null of no autocorrelation of these linear combinations. We apply 3 lags when detecting autocorrelations in what follows.
Applying this approach to our dataset, we do not reject the null for eight linear combinations which can hence be treated as white noise. For the ninth such combination, the p-value for the multivariate Ljung-Box test drops from 0.1935 to 0.0002, so that the white noise hypothesis is rejected for reasonable significance levels. We conclude that there are $s+s_0=21-8=13$ components which account for the dynamic properties of the process. PanYao2008 also propose an estimator for the space of dynamic components $(x_t', z_t')'$. We call these estimates (rotated by principal components) the factors in what follows.
Our model implies that $(x_t', z_t')'$ and hence a suitable rotation of the factors can be modelled as $s+s_0$ univariate time series which are mutually orthogonal at all leads and lags. This corresponds to the notion of dynamic orthogonal components as introduced by MatTsa2011 who provide methods to test for the presence of such a structure and to estimate the appropriate rotation. Using first differences of the factors to achieve stationarity as required by MatTsa2011 for suitable values of $d_j$, we find highly significant cross-correlations of the raw factors (the test statistic takes the value 4198.94 for a level 0.01 critical value of 625.80) while a dynamic orthogonal structure is not rejected for the rotated series, with a test statistic of 445.55 and a corresponding p-value close to one. The test result also holds if the test is conducted in levels. In what follows, the dynamic orthogonal components are computed from the factors in levels which slightly outperforms the difference-approach in simulations with fractional processes.\footnote{Results are available from the authors upon request.}
Due to their dynamic orthogonality, the rotation of MatTsa2011 identifies the single processes in $(x_t', z_t')'$ up to scale, sign and order. A preliminary analysis of the integration orders of $x_t$ can hence be undergone by a univariate treatment of these series. We investigate these integration orders by the exact local Whittle estimator allowing for an unknown mean Shi2010a.
A possible grouping of components with equal integration orders is assessed by the methods proposed by RobYaj2002, with the modifications for possibly nonstationary integration orders by NieShi2007. The specific-to-general approach of RobYaj2002 sequentially tests for existence of $j=1,2,\ldots$ groups of equal integration orders. The sequence is terminated if for some $j^*$ there is a grouping for which within-group equality is not rejected, and for $j^*>1$ the grouping with highest p-value is selected. In our application, we restrict attention to possible groupings where, for $\hat{d}_{i_1}>\hat{d}_{i_2}>\hat{d}_{i_3}$, there is no group including both $i_1$ and $i_3$ but not $i_2$. For the tests of equal integration orders within the sequential approach, we consider the Wald test proposed by NieShi2007, jointly testing all hypothesized equalities for a given grouping. We choose $m=\lfloor n^{0.5} \rfloor = 46$ as bandwidth and set the trimming parameter $h$ to zero, since the dynamic orthogonal components structure does not permit fractional cointegration.
The estimated integration orders for the 13 dynamic orthogonal components range from 0.0087 to 0.7328 and indicate that some of the components may have short memory while others behave like stationary or nonstationary fractionally integrated processes. We clearly reject equality of all integration orders, while also each of the groupings in two groups can be rejected on a 0.01 significance level. For three groups, we do not reject the hypothesis of equal integration orders within groups. The sequential test for groups with equal memory yields $j^*=3$ with a p-value of 0.3181, where groups of three ($\hat{d}^{(1)} = 0.6717$), seven ($\hat{d}^{(2)} = 0.3448$) and three ($\hat{d}^{(3)} = 0.0523$) components are identified, respectively. The hypothesis that $d^{(3)}=0$ is not rejected. We may therefore treat the members of the third group as short-range dependent and belonging to $z_t$. Thus, $s_1=3$, $s_2=7$ and $s_0=3$ appear as a reasonable specification for model (ref) due to the preliminary analysis.
We obtain starting values for the parametric estimator from this procedure. Firstly, $d$ and $\phi$ are estimated from the dynamic orthogonal components. Secondly, from regressing observed data on standardized estimated orthogonal components with unit innovation variance, we obtain starting values for $h$, $\varLambda$ and $\varGamma$, while certain columns of the latter matrices are rotated to satisfy the zero restrictions.
In very high-dimensional cases, the approach of PanYao2008 is not applicable, but LamYaoBa2011 and LamYao2012 provide feasible methods for stationary settings and comment on possible extensions to nonstationarity. In cases where the dynamic orthogonal components specification (ref) is not appropriate, but the general setup (ref) is, a specification search and preliminary estimates for the integration and cointegration parameters of the more general model could be based on the algorithm of Hua2009 which is capable of identifying and estimating cointegration subspaces by semiparametric methods.
We proceed with maximum likelihood estimation of the fractional components model using the EM algorithm of the state space representation. Although the exact state space respresentation is easily obtained using the current type II definition of fractional integration, the state dimension grows linearly with $n$ and becomes computationally infeasible. Instead, the latent fractionally integrated components are mapped to approximating ARMA(3,3) dynamics as described and justified by HarWei2018. There, we show by simulation that low-order ARMA approximations (with parameters depending both on $d_j$ and on $n$) provide an excellent approximation performance and outperform truncated moving average and autoregressive representations by large amounts.
We note that an asymptotic theory for maximum likelihood estimation in the fractionally cointegrated state space setup is not available. Certain functions of the parameter estimates are expected to exert nonstandard asymptotic behavior, especially in the nonstationary case $d_j>0.5$ for some $j$. However, normal and mixed normal asymptotics have been established and conventional tests and confidence intervals have been justified in different parametric fractional cointegration settings as well as in state space models with common unit root components ChaMilPa2009,ChaJiaPa2012. We thus use standard parameter tests in what follows, bearing the preceding caveats in mind.
Constant terms are included by a further column $c$ in the observation matrix and estimated along with the free elements of $\varLambda$ and $\varGamma$. Setting the autoregressive order of $z_t$ to one and using starting values as described above, we estimate models with $q\in \{1,2,3\}$ groups of equal integration orders $d^{(j)}>0$ and additional autoregressive components. The Bayesian information criterion (BIC) is used to select sizes $s_0,\ldots,s_q$ and the value of $q$ with appropriate in-sample fit.\footnote{Instead of estimating all reasonable combinations of $s_0$, \ldots, $s_q$ for each $q$, we begin by the optimal grouping for a given $q$ obtained from the semiparametric methods of the previous section. From this specification, denoted as $s^{\{0\}}_j$, $j=0,\ldots,q$, we estimate all models characterized by $s_j \in \{s^{\{0\}}_j-1,s^{\{0\}}_j,s^{\{0\}}_j+1\}$, $j=0,\ldots,q$, given that they satisfy $s+s_0-1 \geq s_j^{\{0\}}\geq1$. The model with the least value of the BIC is selected and its indices denoted as $s^{\{1\}}_j$, and again models with indices close to $s^{\{1\}}_j$ are estimated and compared. This process is iterated until $s^{\{i\}}_j=s^{\{i-1\}}_j$ holds for all $j=0,\ldots,q$. As a result, also the number of white noise combinations may differ from 8, the result of the semiparametric analysis in the previous section.} We apply the BIC even if consistency is not established in this fractional setting. We expect that existing results hold for specification choices not involving the fractional components, while it is not clear to what extent the results of ChaJiaPa2012 carry over to the fractional setup. There, consistency of the BIC is shown for the number of stochastic trends in a unit root state space model.
We complement the semiparametric results of the previous section by a parametric specification search. After diagnostic checking of the selected model, we will take a closer look at its parameter estimates and implied long-run characteristics. The best models for each $q$ are shown in table (ref), where estimated integration orders are given along with the log-likelihood (log-lik) and the BIC. Regarding the integration orders, we find that for $q>1$ estimates of $d^{(1)}$ are always above 0.5 suggesting nonstationarity of at least $s_1$ series in $y_t$. Overall, the models with $q=2$ are superior, in particular the grouping in $s_1=2$ and $s_2=9$ fractional and $s_0=2$ nonfractional components. This specification is similar to the one selected by the semiparametric approach and also suggests a dynamic dimension of $s+s_0=13$. Interestingly, the same specification with full noise covariance matrix $H$ is inferior ($BIC=-16.626$) as is the model with a full vector autoregressive matrix $\varPhi$ ($BIC=-17.150$). Furthermore, considering more lags in $z_t$ does not sufficiently improve the fit ($BIC=-17.155$ for $k=2$, $BIC=-17.046$ for $k=3$ and $BIC=-17.139$ for $k=4$).
We conduct several diagnostic tests on standardized model residuals $e_{it}= v_{it} / \sqrt{F_{ii,t}}$, where $v_t$ and $F_t$ are filtered residuals and forecast error covariance matrices, respectively. The residuals corresponding to log variances and z-transformed correlations for the first three assets are plotted in figure (ref), while residual autocorrelations are depicted in figure (ref), autocorrelations of squared residuals in figure (ref) and histograms of the residuals along with the normal density in figure (ref). The visual inspection shows some but no overwhelming evidence against the model assumptions. Autocorrelation both of residuals and squared residuals are generally below 0.1 in absolute value and mostly within the $\pm$ 2 standard error bands which are shown as horizontal lines. Some deviations from normality are visible, but not the sort of skewness and fat tails observed for models of untransformed residual variances and covariances.
Table (ref) presents the diagnostic tests on standardized residuals. The p-values are shown for the Ljung-Box test (LM) and the ARCH-LM test for conditional heteroscedasticity (CH) for different lag length 5, 10 and 22. Additionally, the Jarque-Bera test result (JB) is shown in the last column. The null of no autocorrelation is not rejected at the 0.01 level for all but two or three residuals, depending on lag length. Clear evidence of conditional heteroskedasticity is found for the residuals of the log variance series, that is $e_{2t}$,$e_{3t}$, $e_{5t}$, and $e_{6t}$, where also the normality assumption is clearly rejected, but also for a few correlation series such as $e_{15,t}$ or $e_{19,t}$. A more flexible data transformation like the matrix Box-Cox approach of Wei2014 would typically ameliorate these findings, but we do not follow this approach further here.
Estimates of several of the model parameters are shown in table (ref). Along with the maximum likelihood estimates, we also show the mean of the estimators from a model-based bootstrap resampling exercise with 1000 iterations and generally find a low bias for the corresponding estimates. We also show standard errors, obtained in three ways, namely by the bootstrap (SE.boot), using the information matrix Har1991, denoted by SE.info, and by the sandwich form Whi1982, labelled SE.sand in the table. The different methods of computing standard errors give similar results, except for the variance parameters $h_i$, where the sandwich estimates are large compared to the others. Overall, including the parameters not shown in the table, the median ratio between bootstrap and sandwich standard errors is 1.31, while a typical sandwich estimate is 1.20 times larger than the corresponding estimate from the information matrix. We hence use the bootstrap methods in order to avoid a possible underestimation of the variances and spurious inference.
The estimated memory parameters $d_1$ and $d_2$ exert a marked difference in the integration orders of fractional components. The two series in the first group are the cause of significant nonstationarity in our dataset. The second group of nine series introduces stationary long-memory persistence. In contrast, the nonfractional components in $z_t$ are only mildly autocorrelated, with small but significant autoregression parameters. Figure (ref) gives a visual impression of the factor dynamics, showing full sample (smoothed) estimates of the two nonstationary components (above), of the first two stationary long-memory components (middle) and of the short-memory components (below). The $\pm$ 2 standard error confidence intervals suggest a relatively precise estimation of the components. The different persistence of the three groups is clearly visible.
We turn to a discussion of the cointegration properties of the estimated system. In our preferred specification with a cointegration rank of $p-s_1=19$, and an 11-dimensional cointegration subspace, the loadings of fractional components provide an easier interpretation than the corresponding cointegration vectors, although the latter can be easily obtained and suitably normalized.
With the abovementioned caveat that asymptotic properties are not available for this fractional cointegration setting, we show $t$-ratios for constants, for fractional loadings and for nonfractional loadings in table (ref), where the bootstrap standard errors are used. The $t$-ratios for $\varLambda^{(1)}$ suggest that each of the series in $y_t$ is influenced by the nonstationary components, and hence all components of $y_t$ are nonstationary themselves. The first component loads very significantly on all variances with the same sign and can hence be interpreted as the main common risk factor. The second component represents joint common nonstationarity of the correlations, which is negatively associated with the IBM return variances. Except those corresponding to the first, the second and the forth stationary components with their equal signs, the columns of $\varLambda^{(2)}$ have a rather mixed pattern. Like the nonstationary factors, also the $I(d^{(2)})$ components affect variance and correlation dynamics at the same time and therefore induce fractional cointegration between log variances and z-transformed correlations.
The finding of nonstationary fractional components affecting variances and correlations at the same time is new to the literature and may have remarkable consequences on portfolio selection and hedging opportunities, even at longer horizons. These effects should also be relevant to systemic risk measures as considered by central banks and regulators worldwide. To shed further light on the practical value of our approach, we turn to an evaluation of the forecasting precision in a real-world scenario in the following section.
We assess the forecasting performance of our model by means of an out-of-sample comparison. To avoid reference of the forecasts on the out-of-sample periods, we conduct a semiparametric specification search along the lines of section (ref) for the first estimation sample only, i.e.\ for $y_t$, $t=1,\ldots,1508$, while $t=1509,\ldots,2156$ is reserved for prediction and therefore not used for selecting the specification. In this way, the model for the forecasting comparison includes $s_1=2$, $s_2=7$ and $s_0=3$ components of different integration orders. Rather than conducting comprehensive comparisons of a wide range of available methods which is beyond the scope of this paper, we select straightforward and simple benchmark models which have performed well in previous studies.
We choose the same out-of-sample setup as in Wei2014. Thus, for each $T'\in [1508;2156-h]$, various competing models are estimated for a rolling sample with $n=1508$ observations, $y_{T'-1507}, \ldots, y_{T'}$. From these estimates, forecasts of $y_{T'+h}$, $h=1,5,10,20$, are computed. Also in line with Wei2014, we compute bias-corrected forecasts of the realized covariance matrices $\hat{X}_{T'+h|T'}$ by the simulation-based technique discussed there. We evaluate the forecasting accuracy using the ex-post available data of the respective period.
The forecasting precision is assessed using different loss functions defined in appendix (ref). We consider the Frobenius norm $LF_{T',h}$ (ref), the Stein norm $LS_{T',h}$ (ref) and the asymmetric loss $L3_{T',h}$ (ref); see LauRomVi2011 and LauRomVi2013. Additionally, the ex-ante minimum variance portfolio is computed from the forecast and its realized variance $LMV_{T',h}$ (ref) used as a loss with obvious economic relevance. Furthermore, we assess density forecasts $f_r$ of the daily returns using covariance matrices, which are evaluated at the daily returns $r_{T'+h}$ in a logarithmic scoring rule $LD_{T',h}$ (ref).
As benchmarks, we consider two linear models for the log variance and z-transformed correlation series $y_t$, namely a diagonal vector ARMA(2,1) and a diagonal vector AR\-FI\-MA(1,$d$,1) model, which have been found to perform well by Wei2014. Additionally, the diagonal vector AR\-FI\-MA(1,$d$,1) model is applied to the Cholesky factors of the covariance matrices ChiVoe2011. Furthermore, we consider models with a conditional Wishart distribution, namely the conditional autoregressive Wishart (CAW) model of GolGriLi2012, a dynamic correlation specification (CAW-DCC) of BauStoVi2012, and additive and multiplicative components Wishart models as proposed by JinMah2013. For further details on the comparison models consult appendix (ref).
For each loss function and horizon $h$, we compute the average losses (risks) for all models and obtain model confidence sets of HanLunNa2011, bootstrapping the max-$t$ statistic with a block lengths of $\max\{5,h\}$. In tables (ref), (ref), (ref) and (ref), we present the risks for $h=1,5,10,20$. The best performing model ($^{***}$) as well as members of the 80% model confidence set ($^{**}$) and models contained in the 90% but not in the 80% set ($^{*}$) are indicated.
The fractional components model is among the best competitors for all horizons and loss functions. It has lowest risks for almost all setups. Exceptions occur for $h\geq10$ where the ARFIMA model for log variances and z-correlations performs best in some cases. Overall, the ARFIMA model on $y_t$ appears as second best in terms of forecasting precision.
The DOFC model is always contained in the 80% model confidence set whereas all other models are rejected at least in some cases. For the Stein loss and the minimum-variance loss, the DOFC model is significantly superior than most competitors for small horizons, while with the Frobenius and asymmetric loss, rejections of other models are achieved for $h=10$ and $h=20$.
The performance of the fractional components model in terms of density forecasting is noteworthy. In each case there, our model is either the single member or one of two models in the confidence set and hence significantly outperforms most of the competitors. Since the behaviour of future daily returns is usually more important than the realized measures themselves, this finding is particularly strong from a practitioner's perspective.
Overall, we find a very good forecast performance of the model proposed in this paper. Although for some criteria and horizons statistical significance is lacking, the model yields very precise forecasts in relation to different competitors for all considered horizons and for several ways to measure this precision.
We have suggested a general setup and a parsimonious model with very general fractional integration and cointegration properties. We discussed the usefulness of our approach for multivariate realized volatility modeling. In our application it was shown to provide a reasonable in-sample fit and competitive out-of-sample forecasting accuracy.
Several questions remain for further research. From an empirical point of view, we have shown the relevance of a very restricted specification in financial econometrics, but the general setup we introduced has a broader scope. Fractional components models with rich short-run dynamics may be considered for models of smaller dimension. In several empirical setups, fractional integration and cointegration has been found relevant, so that dynamic modeling, forecasting, identification of structural shocks and impulse response analyses in an according framework is a fruitful direction of ongoing research.
The research of Roland Weigand has mostly been done at the Institute of Economics and Econometrics of the University of Regensburg and at the Institute for Labour Market Research (IAB) in Nuremberg. Very valuable comments by Rolf Tschernig, by Enzo Weber and by participants of the Interdisciplinary Workshop on Multivariate Time Series Modeling 2011 in Louvain La Neuve, at the Statistische Woche 2011 in Leipzig, and of research seminars at the Universities of Regensburg, Augsburg and Bielefeld are gratefully acknowledged. The authors are also thankful to Niels Aka for providing R codes to estimate model confidence sets. Tobias Hartl gratefully acknowledges support through the projects TS283/1-1 and WE4847/4-1 financed by the German Research Foundation (DFG).