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.
57,438 characters · 17 sections · 87 citation commands
Sequential Estimation of Multivariate Factor Stochastic Volatility Models
\onehalfspacing \setcounter{page}{1}
Modeling and forecasting the multivariate volatility of financial returns is crucial for risk and portfolio management. Most investors hold large baskets of financial assets whose risks are time-varying and correlated in time. Dynamic approaches to capture, within a multivariate framework, the unobserved time-varying variation and correlation of the financial returns include the multivariate GARCH (MGARCH) models and the multivariate stochastic volatility (MSV) models. While the GARCH specifications are relatively restrictive as they treat the (co-) variance as conditionally deterministic, stochastic volatility models are more flexible as they allow for the variances and correlations to be stochastic. However, the MSV models are difficult to estimate since the (co-) variances and their noises are latent. Moreover, due to the curse of dimensionality, since the number of parameters increases at least quadratically in the dimension of the return vector, MSV models have found only limited application in practice.
A solution to the curse of dimensionality of MSV is to impose a factor structure on the vector of underlying returns, where the factors and idiosyncratic noises follow independent autoregressive stochastic volatility (ARSV) processes. This is known as the multivariate factor stochastic volatility (MFSV) model and has been introduced by harvey94 (further developed by pitt99). This factor representation substantially improves the feasibility of the MSV model in practice by significantly reducing the number of parameters that now increases only linearly with the number of series. However, MFSV still suffers from practical ineffectiveness, even under an exact factor structure, given that the computational burden is enhanced by the latency of the factors, additionally to the variances and idiosyncratic noises.
In this paper, we propose a simple method to estimate the parameters of MFSV models, which allows them to be easily applied in practice, even to very large dimensions and at little computational costs. It consists of two steps: in the first step, we estimate the parameters describing the factor structure, i.e., the loadings and the unconditional variances of the factors and of the idiosyncratic noises, by employing the maximum likelihood (ML) to a static factor representation and by using the convergence results of anderson03, bai12 and bai16 derived for exact factor models with arbitrary dynamics. In the second step, we apply the Efficient Method of Moments (EMM) of bansal94 and gallant96 to estimate the ARSV parameters by implementing simple univariate GARCH auxiliary models to each of the extracted static factors and noises from the first step.
This sequential estimation procedure, which has already been successfully applied by sentana08, calzolari13, halbleib18 and halbleib21 in other contexts, is computationally feasible and simple in terms of implementation and running time regardless of the dimension of the return vectors. To the best of our knowledge, it is also the first frequentist procedure to estimate MFSV models, as all existing ones are exclusively Bayesian. Unlike the Bayesian approaches, this paper focuses solely on estimating MFSV models with fixed parameters and not on filtration and forecasting of the latent variables of the model. We show that our method can be straightforwardly applied to larger dimensions (as well as many factors) compared to the existing literature without imposing any constraints on the parameters and at very low computational costs.
The existing Bayesian procedures, which simultaneously estimate the parameters and filter the latent variables of MFSV, build on the Markov Chain Monte Carlo (MCMC) method as first proposed by pitt99. Although theoretically appealing, this method becomes cumbersome when increasing the dimension of the return vector, as it samples from the entire posterior distribution of the model at each iteration. chib06 partially fix this problem and show that, conditionally on the factor loadings and the factors, the problem reduces to sampling from univariate ARSV processes, which significantly decreases the computational burden of the sampling algorithm and allows for efficient sampling from the posterior distribution also for higher dimensional vectors of returns. Thus, han05 and nardari07 apply the sampling algorithm of chib06 to vectors of around 30 returns and three factors and to vectors of 10 returns and five factors, respectively. kastner17 improves further the efficiency of the MCMC algorithm proposed by chib06, making it computationally faster, but still burdensome in terms of memory requirements, as it involves storing draws from a very high dimensional joint density function of all parameters and of all latent variables (factors, noises, and their respective variances) at each time observation. Thus, the procedure of kastner17 becomes computationally cumbersome for large vector dimensions, many factors, and many observations. For this reason, the algorithm is applied to only up to 26 series and 5 factors. kastner19 makes use of Bayesian shrinkage to constrain the dimension of the loading parameters in order to apply the procedure to larger dimensions, such as 300 returns. Nevertheless, as stated in the original paper, the possible parallelization makes the algorithm computationally more intensive in terms of memory requirements.
We show both in our simulation exercise and empirical application that the frequentist sequential method we propose in this paper can estimate the parameters of the MFSV model for very large dimensions of the return vector within minimal computational time and with no memory constraints, especially when parallelizing the equation-by-equation estimation of the second step of the procedure. In the simulation exercise we show that, even without parallelization, we can obtain accurate estimates for vectors of returns of dimension 10 and 2 factors in less than 13 seconds if feeding starting values close to the true parameter values and less than 82 seconds if the starting values are random. Moreover, we find the optimal number of simulations in the EMM estimation that increases the efficiency of the estimates with no additional computational costs as it is adaptive to the length of the series.
In the empirical application we provide further evidence on the speed of our procedure when estimating MFSV for a vector of 148 returns and up to three factors. Here, we show that, as expected, the computational advantage of our method becomes more pronounced once we use multiple CPU cores to run the second step of the estimation in parallel. With parallelization enabled, the computational advantage in the real data application increases by a factor larger than 200. While in simulations our method converges in seconds, in the empirical application it may take only up to 10 minutes. Thus, the method we propose in this paper is the first that manages to practically estimate large dimensional MFSV models without any (shrinkage) constraints on the parameters. While the filtration of the latent variables does not make the object of this paper, based on the estimates proposed here, it takes in the simulation up to approximately 160 minutes using the Bootstrap Particle Filter of gordon93, making it at least ten times faster than the most efficient Bayesian method, which runs roughly 24 hours. In the empirical application the Bootstrap Filter is not applicable due to the high number of latent series bengtsson08. We, therefore, apply the Auxiliary Multiple Particle Filter of muecher23 to the data. The run time of the filter with $10240$ particles per series is approximately 20.5 hours, making it roughly 2.9 times faster than the Bayesian method, which runs more than two days for the empirical exercise.
The remainder of this paper is structured as follows. Section 2 formally introduces the MFSV model. In Section 3 we present the estimation method proposed in this paper. Section 4 presents the results from a comprehensive simulation study, and in Section 5, we present the results from applying our procedure to real data. Section 6 concludes.
harvey94 and pitt99 introduce factor models in the context of MSV modeling to capture the common dynamics and conditional heteroskedasticity in financial returns and to reduce the curse of dimensionality of standard MSV models.
The MFSV is defined for $\bm{y}_{t}$, a vector of $N$ return series observed at time $t$ such that
where $\bm{f}_{t}$ is a vector of $k\times 1$ latent factors, independent of each other, $B$ is of dimension ${N\times k}$ and contains the factor loadings with $k \leq N$ and $\operatorname{rank}(B) = k$ and $\bm{\epsilon}_{t}$ is a $N\times 1$ vector of idiosyncratic components that are orthogonal to the factors. $\Sigma_{t}$ is the diagonal covariance matrix of $\bm{\epsilon}_{t}$ with the elements $\sigma_{t,1}^{2}, \hdots, \sigma_{t,N}^{2}$. $\Gamma_{t}$ denotes the diagonal covariance matrix of the $k$ factors with the elements $\gamma_{t,1}^{2},\hdots,\gamma_{t,k}^{2}$. For identification, we restrict $B$ such that $b_{jj} = 1 \text{ and } b_{ij} = 0 \, \forall j > i, \hspace{3mm} j = 1,\hdots,k$ and $i=j, \ldots,N$ pitt99,chib06,nardari07.\footnote{sentana01 show that assuming heteroskedasticity for the factors, one may relax one of these restrictions. However, given that our estimation strategy presented in the next section relies on estimating a static factor model, we need both restrictions on the factor loadings matrix, paired with the restrictions that $\Sigma_{t}$ and $\Gamma_{t}$ are diagonal (also implying that their unconditional counterparts, $\mathbb{V}[\bm{\epsilon}_{t}]$ and $\mathbb{V}[\bm{f}_{t}]$, are diagonal) in order to achieve full identification in the static factor model bai12,cox17,williams20.} Thus, the loading parameters left for estimation are $b_{ij}$ with $j=1,\hdots, k$ and $i=j+1,\hdots, N$.
Let $\bm{x}_{t} = (\epsilon_{t,1},\hdots,\epsilon_{t,N},f_{t,1},\hdots,f_{t,k})'$ such that each of the elements of $\bm{x}_{t}$ follows a univariate Autoregressive Stochastic Volatility (ARSV) process as it follows:
where
$|\varphi_{m}|<1$ that assures covariance-stationarity and $\sigma_{t,m}^{2}>0$. Thus, the conditional covariance matrix of the return series at time $t$ is given by:
where $\bm{h}_{t} = (h_{t,1},\hdots,h_{t,N+k})^{\prime}$ is the vector of log variances of the factors and idiosyncratic errors. The model has a total of $Nk - k(k+1)/2 + 3(N+k)$ parameters that increases linearly in the number of return series $N$ rather than quadratically, as is the case of the MGARCH and general MSV models. However, the estimation of the parameters of the MFSV is difficult since the factors $(f_{t,1},\hdots,f_{t,k})'$, the errors $(\epsilon_{t,1},\hdots,\epsilon_{t,N})'$ and their stochastic variances $(\gamma_{t,1}^{2},\hdots,\gamma_{t,k}^{2},\sigma_{t,1}^{2},\hdots\sigma_{t,N}^{2})'$ are latent.
In what follows, we present our simple frequentist approach to estimate the parameters of the MFSV model based on two steps. For this reason, we rewrite the vector of the model parameters as $\bm{\theta} = \left(\bm{\theta}_{1}, \bm{\theta}_{2}\right)^{\prime}$, where:
is estimated in the first step (to be described in Section (ref)) and
is estimated in the second step (to be described in Section (ref)).
Note that $\bm{\theta}_{2}$ does not contain the ARSV constants $(\mu_{1},\hdots,\mu_{N+k})'$, since we can identify them from the unconditional variances of the ARSV processes as it follows:
where $(\psi_1,\ldots,\psi_{N+k})'= \left(\sigma_{1}^{2}, \hdots,\sigma_{N}^{2}, \gamma_{1}^{2},\hdots,\gamma_{k}^{2}\right)'\equiv\bm\psi$.
This section provides the estimates of $\bm{\theta}_{1}$. For this we define a static factor model as it follows:
with the same assumptions and restrictions on $B^{*},\bm{g}_{t}, \Sigma^{*},\Gamma^{*}$ as in the MFSV model presented in Section 2. We achieve full identification of the static factor model from the restrictions on the factor loadings matrix and the unconditional covariance matrices of the static factors $\bm{g}_{t}$ and of the static factor model idiosyncratic errors $\bm{e}_{t}$. For details on the identification see williams20. Denote the parameter vector of this static factor model by
i.e., $\bm{\beta}_{1}$ gives the factor loadings, the unconditional variances of the factors, and the unconditional variances of the idiosyncratic errors of the static factor model defined in equations ((ref)) - ((ref)). bai12 and bai16 show that, in the exact factor model case, for $T \to \infty$,
i.e., the static factor model ML estimators converge to their true parameter counterparts in the dynamic factor model.\footnote{bai16 show that for approximate factor models, the convergence occurs for $N,T \to \infty$.} This result holds if the factors are independent of the error terms, if the factor variances are bounded away from zero (bai12; bai16) and if the distribution of $(\bm{\epsilon}_{t}^{\prime},\bm{f}_{t}^{\prime})^{\prime}$ has finite second-order moments anderson03, which are also the assumptions of our factor structure. This result simplifies our estimation procedure drastically, especially when applied to vectors of returns and factors of large dimensions.
However, as shown by bien11 and bai12, the log-likelihood of the static factor model is multi-modal. Thus, classical numerical optimization methods are infeasible for higher dimensional return vectors $\bm{y}_{t}$ as they are likely to get stuck in local minima. In this paper we use the Expectation Maximization (EM) algorithm proposed by bai12 to estimate the parameters of the static factor model combined with a modified version of the Gradient Descent algorithm of bien11, which considerably reduces the computational time (see also bai16b and daniele19 for a successful implementation of this combination to estimate other types of factor models). The detailed description of the estimation procedure at this step is presented in Appendix (ref).
After obtaining $\widehat{B^{*}}$, $\widehat{\Sigma^{*}}$ and $\widehat{\Gamma^{*}}$, we make use of the asymptotic results of bai12 and bai16, and take them as consistent estimates of ${B}$, $\Sigma$ and ${\Gamma}$, i.e. of $\bm{\theta}_{1}$, by simply setting: $\hat{\boldsymbol{\theta}}_{1} = \hat{\boldsymbol{\beta}}_{1}$.
Based on $\widehat{B^{*}}$, $\widehat{\Sigma^{*}}$ and $\widehat{\Gamma^{*}}$ obtained in the first step described above, we extract the static factors by the projection formula of bai12:
where $\bar{\bm{y}}$ is the mean of the return series over the whole sample. Define:
and let $\hat{\bm{x}}_{t} = (\hat{e}_{t,1},\hdots,\hat{e}_{t,N}, \hat{g}_{t,1},\hdots,\hat{g}_{t,k})'$.
At this second step, we estimate the elements of $\bm{\theta}_{2}$ by applying the EMM approach of bansal94 and gallant96. Given the independence among and between the factors and the noises of the factor structure defined in Section (ref), the EMM approach can be simplified by estimating the parameters of the ARSV models equation-by-equation by applying the appropriate univariate auxiliary model to each of the series composing $\hat{\bm{x}}_{t}$. We choose here to apply GARCH(1,1) as the auxiliary model, as it has already been successfully implemented for the EMM estimation of univariate ARSV models by calzolari04 and monfardini98, among others.\footnote{As an alternative to GARCH, we also implement the ARMA(1,1) model fitted to the log squared transformation of the elements of $\hat{\bm{x}}_{t}$ monfardini98. However, due to the "flatness" in the ARMA(1,1) likelihood for at least one of the elements of $\hat{\bm{x}}_{t}$, we were not able to get reliable estimates and standard errors of the MFSV parameters based on this choice of the auxiliary model. The results are, however, available from the authors upon request.}
In Appendix (ref) we provide the detailed description of the steps taken to undergo the EMM estimation of $\bm{\theta}_{2}$. For this we define here below the structure of the GARCH(1,1) that we apply to each of the series composing $\hat{\bm{x}}_{t}$ as it follows:
where $m=1,\ldots,N, N+1, \ldots, N+k$, $\xi_{t,1},\ldots,\xi_{t,m}$ are independent white noises, $\alpha_{1,m}>0$, $\alpha_{2,m}>0$ and $\alpha_{1,m}+\alpha_{2,m}<1$ $\forall m$ and $\hat\psi_{m}$ is the estimate of $m$-th element of $\bm\psi$ computed from the estimated elements of $\bm{\theta}_{1}$ obtained in the first step. Here, we follow sentana01 and sentana08 and identify GARCH constant through a consistent estimator of the unconditional variance of the respective process. We define $\bm{\beta}_{2}=(\alpha_{1,1},\ldots,\alpha_{1,N+k},\alpha_{2,1},\ldots,,\alpha_{2,N+k})'$ containing all parameters for the $N+k$ univariate GARCH(1,1) models described above that are estimated equation-by-equation by means of Pseudo Maximum Likelihood (PML). To speed up the estimation at this step, we can parallelize the estimation of $\bm{\theta}_{2} = \left(\varphi_{1},\hdots,\varphi_{N+k},\sigma_{\eta,1}, \hdots,\sigma_{\eta,N+k}\right)'$ in $N+k$ EMM estimation algorithms for each of the parameter pairs $\bm{\theta}_{2,m} = (\varphi_{m},\sigma_{\eta,m})'$ with $m=1,\ldots, N+k$ to which corresponds the auxiliary parameter vector $\bm{\beta}_{2,m}=(\alpha_{1,m},\alpha_{2,m})'$. The EMM estimation provides consistent estimators of $\bm{\theta}_{2}$ under very general conditions that are fulfilled by our model specification.
Both parts of the estimation procedure we present above involve the estimation of an auxiliary model: i.e., in the first part, the static factor model, in the second part, univariate GARCH models on the extracted static factors from the first part. While the estimates of the auxiliary model of the first part $\hat{\bm{\beta}}_{1}$ are taken to be the final corresponding estimates (of loadings and unconditional variances) of the true dynamic model, the auxiliary estimates of the second part $\hat{\bm{\beta}}_{2}$ help at computing the EMM estimates of $\bm{\theta}_{2}$. However, in order to compute the standard errors of the estimators of the MFSV parameters we propose here, we implement the asymptotic variance-covariance matrix of the EMM estimator given by gourieroux93:
where $\widehat{\mathcal{I}}(\hat{\bm{\beta}})$ denotes the Fisher information matrix of the auxiliary model and $ \mathcal{Q}\left(\tilde{\bm{y}}_{t}(\hat{\bm{\theta}}),\hat{\bm{\beta}}\right)$ is obtained by stacking vertically together the auxiliary score vectors $\mathcal{Q}(\bm{y}_{t};\hat{\bm{\beta}}_{1})$ and $\mathcal{Q}(\hat{\bm{x}_t};\hat{\bm{\beta}}_{2})$, i.e., $\mathcal{Q}(\bm{y}_{t};\hat{\bm{\beta}}_{1},\hat{\bm{\beta}}_{2}) =(\mathcal{Q}(\bm{y}_{t};\hat{\bm{\beta}}_{1})^{\prime},\mathcal{Q}(\hat{\bm{x}_t};\hat{\bm{\beta}}_{2})^{\prime})^{\prime}$, where $\mathcal{Q}(\bm{y}_{t};\hat{\bm{\beta}}_{1})$ are obtained by computing the first derivative of the log-likelihood of the static factor model with respect to ${\bm{\beta}}_{1}$ at $\hat{\bm{\beta}}_{1}$. The closed form expression of this score vector is presented in Appendix (ref). $\mathcal{Q}(\hat{\bm{x}_t};\hat{\bm{\beta}}_{2})$ stacks together the $N+k$ GARCH(1,1) scores vectors: $\mathcal{Q}(\hat{\bm{x}_t};\hat{\bm{\beta}}_{2})=(\mathcal{Q}_{1}(\hat{\bm{x}}_{1};\hat{\bm{\beta}}_{2,m}), \ldots,\mathcal{Q}_{N+K}(\hat{\bm{x}}_{N+K};\hat{\bm{\beta}}_{2,N+K}) )'$ defined in Appendix (ref). \footnote{Note that $\mathcal{Q}(\bm{y}_{t};\hat{\bm{\beta}}_{1},\hat{\bm{\beta}}_{2})$ does not contain $\hat{\bm{x}_t}$, since it is a function of $\bm{y}_{t}$ and $\hat{\bm{\beta}}_{1}$} We define $\hat{\bm{\beta}} = (\hat{\bm{\beta}}_{1}^{\prime},\hat{\bm{\beta}}_{2}^{\prime})^{\prime}$.
In Equation ((ref)), the Fisher information matrix of the auxiliary model, $\widehat{\mathcal{I}}(\hat{\bm{\beta}})$ is not available in complete form, as the parameters $\hat{\bm{\beta}}$ are estimated in two steps. Here we follow halbleib21 and replace it with a consistent simulation-based estimator: the sample variance-covariance matrix of 1000 independently simulated score vectors of the overall auxiliary model. These score vectors are computed after the last iteration upon convergence. As all components of the auxiliary model score vector are available in closed form, the simulation-based estimation of the Fisher information matrix and the whole variance-covariance matrix $W(H)$ in Equation ((ref)) are computationally very fast.
As one may see from Equation ((ref)), $W(H)$ decreases with $H$, the number of the simulated series in the auxiliary EMM estimation. As there is no rule on how to choose $H$, in our empirical application we follow the main stream of the literature, and choose $H$ in such a way that the precision of the estimates improves, but not at very high computational costs (usually $H$ is set to 10). However, in our simulation exercise, we find also an optimal choice of $H$ adapted to the number of observations available for the empirical study.
In order to get the standard errors of the constant terms of our ARSV processes $\hat{\mu}_{1},\hdots,\hat{\mu}_{N+k}$ are computed from the estimates of first step estimation, we use the Delta Method as described below: Remember that
and denote $\bm{\xi}_{m} = (\psi^{*}_{m},\varphi_{m},\sigma_{\eta,m})^{\prime}$. The derivative of $\mu_{m}$ with respect to the parameter vector $\bm{\xi}_{m}$ evaluated at $\hat{\psi}_{m},\hat{\varphi}_{m},\hat{\sigma}^{2}_{\eta,m}$ is given by:
The variance of $\hat{\mu}_{m}$ is then given by
where $\mathbb{V}[\hat{\bm{\xi}}_{m}]$ is extracted correspondingly from $W(H)$ defined above.
This section provides results on the statistical properties of the estimates obtained based on the procedure described in the previous section for various choices of $T$, $N$ and $k$ that are empirically relevant. We use $ N = \{10,20,30,100\}$ simulated return series, $k = \{1,2,3\}$ simulated factors to generate $ T = \{1000,4000,10000\}$ observations. We choose the parameter values similarly to kastner17:
\\[1em] The dots indicate an evenly spaced grid between the first and the last value. Therefore, we choose the elements of the first column of $B$ such that they are evenly spaced between the values $0.9$ and $0.1$; the elements of the second column are chosen such that they are evenly spaced between $0.2$ and $0.8$; the elements of the third column are chosen to be evenly spaced between $0.1$ and $0.7$, such that $\varepsilon=(0.7-0.1)/(N-4)$. For the simulation cases with $k<3$, we choose the first $k$ columns of $B$ and the first $N+k$-th values from each of the parameter vectors ${\bm\varphi}$, ${\bm\mu}$ and ${\bm\sigma}_{\eta}$. For all simulations, we choose $R = 1000$ replications.
For the step size of the Gradient Descent Method of bien11, we use the adaptive moments algorithm with $d = 0.005$.\footnote{The estimation is conducted on the bwHPC Cluster with multiple jobs, each using 40 CPU cores (two Octa-core Intel Xeon E5-2670 (Sandy Bridge) Processors with 2.6GHz) in parallel.}
We compare two strategies for selecting the starting values for the parameters of the EMM procedure of the second step. The first strategy is to provide user-specified starting values where we choose them to be 20% smaller than the true ARSV autoregressive parameters and 20% larger than the ARSV standard deviation parameters. The second strategy is to obtain starting values from the QML approach of ruiz94 applied to the extracted static factors and residuals, i.e., we apply the Kalman filter to their log squared transformations. This second strategy is especially appealing for empirical applications where the choice of good starting values is not straightforward.
Next, we analyze the effect that $H$, the number of simulation paths in the EMM, has on the estimation results. We first set $H = 10$ as in calzolari14, halbleib18 and halbleib21. However, we find that, for small $T$, increasing $H$ to $100$ improves the performance of our procedure, whereas, for large $T$, the improvement in the performance is negligible while the computational burden hugely increases.\footnote{We do not report the results for $H = 100$. They are available from the authors upon request.} To balance this behavior, we follow the idea of monfardini98 and choose, besides fixing it to the value of 10, $H$ as a function of the sample size, $T$. Specifically, we set $H$ such that $H\cdot T = 10^5$, i.e. $H = 10^{5}/T$. This is mainly motivated by the reduction in the number of outliers of the EMM estimation procedure, as described below.
We summarize the results of the simulation study by means of the Mean Squared Error (MSE) of the estimated parameters $\hat{\bm{\theta}}^{r}$ over the $r = 1,\hdots,R$ MC replications. We compute the MSE of $\hat{\bm{\theta}}$ by
i.e., the average squared deviation of the estimated values from the true values. We also report results on the ratio between the empirical and asymptotic standard deviations of the estimated parameters.
Some of the replications yield outliers in the second step of the estimation. Such outliers are often characterized by very large estimated ARSV standard deviation parameters or by negative ARSV autoregressive parameters. Not all cases of outliers are accompanied by large values of the EMM distance, and therefore their identification based on it is very sensitive to the distance threshold we choose. Therefore, we choose a less sensitive way to identify outliers based on the values of the estimated parameters in each MC replication. In particular, we discard the replication with at least one of the following outliers. First, the estimated autoregressive parameter is at least ten times smaller than the true value, or it is negative. Second, the absolute value of the constant or standard deviation parameter is at least ten times larger than its true value. \footnote{In fact, we discard replications for which the absolute value of the estimated constant of the factors is greater than 9, since the true value is zero.}Although these types of outliers only affect the estimated values of the parameters of the corresponding series in the second step procedure (and not the estimates of the first step), they also affect all standard errors of the entire vector of parameter estimates due to the way we compute the variance-covariance matrix of the estimates as described in Section (ref).\footnote{Alternatively, we could choose not to discard the outliers and present summary results by means of quantile measures, such as median and the interquartile ranges. These results are available from the authors upon request.}
Tables (ref) to (ref) in Appendix (ref) report the percentage of outliers for different numbers of factors $k$, number of return series $N$, sample size $T$, number of simulations in the EMM step $H$ and starting values. From analyzing the entries of the tables, we see that the choice of $T$ is crucial in generating outliers: while large values of $T$ are accompanied by only a few outliers (for $T=10000$, only up to 0.5% of all simulated series are discarded), working with series of smaller sample sizes is much more affected by the outliers: e.g., for $T=1000$, large $N$ ($N=100$), small $H$ (e.g., $H=10$) and poor starting values, we need to discard up to 37% of the replications. However, increasing $H$ (to 100) reduces the number of outliers drastically, e.g., only 10% are discarded. A further significant reduction is obtained by choosing the QML starting values: for $H=10$ the percentage of discarded replications decreases to 5%, while for $H=100$, it further reduces to 3.5%. Also, for $T=4000$, we find that increasing $H$ and choosing the QML starting values significantly reduces the outliers, especially for large $N$. Moreover, by setting $H=10^5/T$ and using QML starting values, we need to discard only up 3.6% of the replications, regardless of the choice of $T$ and $N$, which is much less than by keeping $H$ fixed (e.g., to $H=10$). Therefore, in our empirical application presented in Section (ref) we choose $H=10^5/T$ and use the QML starting values.
Appendix (ref) and (ref) report simulation results for the whole vector of the MFSV parameters, i.e., $\hat{\bm{\theta}}$ as well as for sub-vectors of it, i.e., for the factor loadings, the unconditional variances and the ARSV parameters, respectively. In Appendix (ref) we focus on presenting the results on the ratio between the MC standard deviation of the estimated parameters and their asymptotic standard errors, computed as described in the Section (ref), averaged over the components of the corresponding (sub)-vectors. Tables (ref) and (ref) present results for fixed $H=10$ and tables (ref) and (ref) report results for $H = 10^{5}/T$. As expected, increasing $H$ in the EMM estimation procedure reduces the asymptotic variance of the estimated parameters. For $H = 10$ and large $T$, the average ratio between the empirical standard deviation of the estimated parameter values over the Monte Carlo replications and the estimated asymptotic standard errors is close to 1, regardless of the choice of the starting values. For lower values of $T$, however, the ratio is smaller than 1 in most cases, indicating that the asymptotic standard errors are too large compared to the Monte Carlo standard deviation. In these cases, the QML-based starting values improve the ratio towards 1. However, by setting $H = 10^{5}/T$, also for smaller values of $T$, the ratios become close to 1.
Appendix (ref) provides detailed numerical results for the MSE of the estimated parameters for the two choices of $H$ and the starting values. From the tables one may see that, as expected, the MSEs of the parameters estimated in the first part of our procedure are unaffected by the choice of the EMM starting values of $H$. However, for the quality of ARSV parameters estimated in the second part of our procedure, both the choice of $H$ and of the starting values is important and exhibits complementary effects for larger $T$s. Thus, choosing QML starting values of $T$ larger than $1000$ provides ARSV estimates with small MSE also in the case of fixed $H=10$.
Generally, the MSE of the parameter vector $\hat{\bm{\theta}}$ decreases when $T$ increases and increases when $k$ increases. These effects are particularly pronounced for smaller values of $N$ and reduce when $N$ increases. Similar results are obtained by halbleib21 that show that the burden of adding a new factor to the estimation reduces when the underlying number of returns increases, as it provides more information on their commonalities than when increasing $T$.
By looking in detail at each set of parameters, one may see that increasing $N$ has a stronger positive effect on the quality of the ARSV estimates than on the quality of the loading and the unconditional variance estimates. While the quality of the loading estimates is mainly improved by increasing $T$, the quality of the unconditional factor variances also improves by increasing the number of factors. As the constant parameters are identified from the unconditional variance of the error term and the AR and the standard deviation parameters, it inherits their MSE pattern, especially in what regards the dependence to $k$ and $N$.
Figure (ref) plots the MSE values for $H = 10^{5}/T$ and QML starting values. The figure plots the MSE against the number of observations $T$ for different numbers of returns series $N$. The numerical MSE results for this case are detailed in table (ref) in Appendix (ref). The figure summarises the MSE results averaged over all parameters. In a nutshell: increasing $N$ and $T$ significantly reduces the MSE, while increasing $k$ has a detrimental effect on the MSE, especially if $N$ and $T$ are small.
Due to the identifying assumptions we make for our model, the ordering of the returns in the vector $\bm{y}_{t}$ matters since the first series only loads on the first factor, the second series only loads on the first two factors, and so on. While the factor loadings matrix, the unconditional variance-covariance matrix of the factors and the constant parameter of the factor ARSV models change by changing the order, the other set of parameters should not change and this is confirmed in a small Monte Carlo simulation study where we invert the original ordering of the vector of returns. The results, which are computed with $H=10^5/T$ and QML starting values, can be obtained from the authors upon request.
One very important advantage of our method is its computational efficiency. Figures (ref) and (ref) plot the computational time per replication for $H=10$ and $H = 10^{5}/T$, respectively against $T$ for different numbers of factors and return series, however by implementing the QML starting values for the EMM estimation. The upper panel shows the average computational time over the Monte Carlo replications for $N = 100$ return series for one, two, and three factors, while the lower panel depicts smaller dimensions, i.e., $N=10, 20, 30$.
From the two graphs one may see that the computational time increases linearly in $T$. Moreover, surprisingly the computational time for $H = 10^{5}/T\geq 10$ is only maximally twice the time needed for $H=10$. For instance, for $T=4000$ (i.e., $H = 25$), there are almost no differences in the computational time between the two choices of $H$. Therefore, one can conclude that adapting the choice of $H$ to $T$ is optimal for any empirical application as it increases the efficiency of the estimates, reduces the number of outliers and is computationally feasible.
In order to present results on the computational time of our procedure we focus here on the simulation design of kastner17 and choose $N = 10$, $k = 2$ and $T = 1000$ as well as their choice of parameters. Our procedure estimates the parameters of the MFSV model under the settings of kastner17 on average in 13 seconds when we use the QML starting values and set $H=10$. Setting $H = 10^{5}/T$- and with poorer starting values, our procedure takes, on average, 82 seconds (roughly 1.37 minutes). Therefore, we can conclude that our procedure is very effective and computationally feasible in estimating the parameters of MFSV even for very large vectors of returns as also illustrated in the empirical application presented below.
Although it does not make the object of this paper, we also report here some computational time from filtering the unobserved series of the model, namely: $15.7$ minutes using $5000$ particles, $28.43$ minutes using $50000$ particles and $158.2$ minutes using 500000 particles. We obtain these numbers replicating the simulation study of kastner17, i.e., the case for $n = 10$, $k = 2$, and $T = 1000$ for $1000$ Monte Carlo replications. We perform the filtering using a straightforward implementation of the Bootstrap Particle Filter of gordon93. In comparison, the Bayesian procedure of kastner17, which estimates the parameters and filtered values simultaneously, requires 1446.208 minutes (roughly 24 hours) using $550000$ sampling periods, of which we discard the first $50000$ as burn-in periods.\footnote{We replicate their MC study in $R$ using the factorstochvol package provided by the authors. We run the code on the bwHPC Cluster. Each job uses 1 core of a Octa-core Intel Xeon E5-2670 (Sandy Bridge) Processors with 2.6GHz), since the package does not allow for parallel execution.}
In this section, we apply our proposed algorithm to a vector of $N = 148$ log-returns of large-cap stocks from the S&{P}500 from January 1981 till December 2018 (a total of $T=9584$ observations), ordered by market cap.\footnote{We provide a list of the ticker symbols and company names ordered by the market capitalization of the stocks in Table (ref) in Appendix (ref).} We demean and standardize the returns such that each series has zero sample mean and unit sample variance to facilitate the estimation. Since we have close to $T=10000$ observations, we use the adaptive choice of $H=10^5/T$ and set, thus, it to $H =10$.
As our procedure only regards the estimation of the MFSV model for a prespecified number of factors, we apply the test of onatski09 to find the number of dynamic factors in our dataset. The test is implemented iteratively with the hypotheses
where $k_1$ and $k_2$ are pre-specified. We set $k_2=3$ and $k_1=0,1,2$. We run the test using the Matlab code provided by onatski09 and present the results in Table (ref). The table clearly shows that we do not reject $H_0$ that $k=1$ at 10% significance level. Figure (ref) in Appendix (ref) shows the number of principle components of the data set we use in this empirical application. The cutoff between the first and the rest of the principal components is pretty clear, which confirms the results of the test of onatski09. That is why we present below first the empirical results for estimating the MFSV model with $k=1$ factors and, afterwards, we discuss the results from estimating the model for $k=2$ and $k=3$.
Figures (ref) to (ref) in Appendix (ref) display the estimated parameters of the ARSV dynamics of the idiosyncratic error term series. While the estimates of the constant parameter seem to vary less among the 148 idiosyncratic noises, the estimates of the autoregressive parameter and of the noise standard deviation display much more variation. Most of the idiosyncratic noises have a very persistent stochastic volatility, with the AR parameter being close to 1, while others display less persistency in their dynamic variation, with the AR parameter being around 0.8. Also the standard deviation of their SVs vary largely among the series from values close to zero to values larger than 0.3. This mirrors the large heterogeneity in the dynamics of the 148 returns considered in the series study
The AR parameter in the SV structure of the factor is estimated to be around 0.9752, while the constant and the standard deviation of the noise are estimated to be -1.5866 and 0.1863, respectively. All parameters of the factor SV are statistically significant at 5% significance level.
In Table (ref) we show the computational time needed for our estimation procedure to estimate the parameters of the MFSV model with one factor on the vector of 148 daily returns of this empirical application. When unable that the estimation of the second step runs in parallel, the run time is roughly 3,5 hours. Parallelizing the second step on 40 CPU cores, the run time reduces to approximately 10 minutes. Because we have no comparison of our estimation time to any other purely estimation algorithm in the literature, we can only present here the time of convergence of the Bayesian procedure of kastner17, which is roughly 59 hours. In order to make the comparison complete, we also present some preliminary results from filtering the latent variables based on our estimates. The classical Bootstrap Particle Filter by gordon93 we use for the runtime comparison for the $n=10$ and $k=2$ case in the simulations, suffers from the curse of dimensionality and therefore does not yield reliable filtering results. However, in muecher23, we introduce the Auxiliary Multiple Particle Filter for the MFSV, to deal with this problem. Applying this filter here with 10240 particles per series (i.e., a total of 1525760 particles) takes 20.5 hours, making the whole procedure roughly 2.9 times faster than the Bayesian procedure. However, the filter uses seven times more particles than the Bayesian procedure uses sampling iterations. \footnote{ We use the factorstochvol package provided by the authors for the estimation and, as suggested by them, we use the deep interweaving code. After 20000 burn-in periods we use 200000 draws from the posterior distribution. The authors use 50000 burn-in and 500000 sampling periods in their paper. Our choice of the burn-in and sampling periods is due to RAM and runtime limitations on our computational device.}
Although as suggested by the test of onatski09 and PCA, we only need one factor to capture the commonality in the stochastic volatility dynamics of the series, we present here, for completeness, the results from estimating the MFSV model for more factors. We choose $k=3$ to be the maximum, as it has also been a common choice in the existing Bayesian estimation literaturehan05,nardari07. We report the estimation results of our procedure in Appendix (ref), namely in Table (ref) for estimating the SV parameters of the factors and in figures (ref) to (ref), for estimating the SV parameters of the idiosyncratic error series.
The first finding of applying our estimation procedure to MFSV models with more factors that the "true" ones is that it detects this mismatch by delivering not significant estimates or anomalies in the estimates, as reported in Table (ref): while the pattern of the estimates of the SV model for the idiosyncratic noises is similar to the ones of $k=1$ and the AR parameters of the SV model of the factors are all very close to 1, which indicate high persistency in the volatility dynamics of the series, similar to the case $k=1$, most of the SV parameters of the second and third factor are estimated to be not significant and some of them have unusual high standard errors. This might be due to the fact that the accurate ML estimation of a static factor model in the first step requires some degree of variation in the factors (bai12; bai16). Suppose this is not the case, as it happens when estimating more factors than necessary, which in fact, should have variation zero. In that case, the rank condition on the factor loadings matrix is violated anderson56, and the distribution of the static factors is not identified cox17. In other words, estimating factors with almost zero variance may indicate that the model is not correctly specified and that there are fewer factors describing the true structure. Therefore, the empirical results for the case of $k=2$ and $k=3$ should be handled with care as they refer to misspecified MFSV models.
To bring further evidence in this direction, we run a small simulation study where we simulate $N=10$ and $N=100$ return series with $T = 4000$ observations and $k = 2$ factors. Figures (ref) and (ref) in Appendix (ref) show the first 10 principal components of the $R=1000$ simulated data sets for $N = 10$ and $N = 100$ respectively. The plots clearly indicate the number of factors to be equal to 2, which is precisely the number of factor we simulate. On this simulated series, we apply our estimation procedure by choosing $k = 1,2,3$ factors. We find that the Fisher information matrix of the auxiliary model for $k=3$ is not of full rank with serious effects on the standard errors and the inference based on them. In fact for $k=1$ and $k = 2$ the estimation works smoothly, while for $k=3$, the first step estimation reveals zero variance of the third factor, as depicted in Figure (ref) in Appendix (ref) for $N=10$ and in Figure (ref) for $N=100$.
The fact that for $k=2$ and $k=3$ we deal with misspecified MFSV models for our real dataset is also detectable when running the Bayesian procedure of kastner17, which in fact, sets the constant in the factor ARSV model to 0 and instead specifies a scaling parameter for the variance of the factor. The variance of the factor is multiplied by this scaling parameter broto04. In the specific context of the MFSV model, the scaling parameters are multiplied then with the respective column of the loadings matrix, allowing for diagonal elements that are not restricted to 1. Figures (ref) and (ref) in Appendix (ref) show that, for $k=3$ the scaling of the second factor is ill-behaved. Specifically, we see in Figure (ref) that the Markov Chain does not converge for the scaling factor and Figure (ref) shows that the estimated marginal density of the scaling factor is bi-modal, with a lot of probability mass around zero. Thus, the Bayesian procedure also faces "problems" when the number of factors becomes larger than the true one.
This paper proposes a two step procedure to estimate multivariate factor stochastic volatility models with exact factor structure and stochastic volatility dynamics for both the factors and idiosyncratic noises. In the first step we use use the convergence results of bai12 and bai16 to estimate the factor loadings and the variances of the idiosyncratic errors and of the factors by means of ML applied to a static factor model. In the second step we estimate the parameters of the stochastic volatility of the factors and of the idiosyncratic noises equation-by-equation by using the EMM procedure of bansal94 and gallant96 and univariate GARCHs as auxiliary models.
Our procedure is the first in the literature to implement "frequentist" techniques to estimate such a model. Due to the two-step procedure and equation-by-equation estimation in the second step, the procedure we propose is computationally very fast and can easily be applied to very large dimensions of vectors of returns and number of factors. The computational efficiency can be further increased by choosing appropriate starting values and parallelizing the EMM estimation in the second step.
In a comprehensive simulation exercise we show that our estimation provides accurate estimates of the parameters as well as how to choose the optimal number of simulations in the EMM procedure in order to increase the efficiency of the estimation at minimal computational costs. In the empirical application to a large vector of 148 returns we provide further evidence of the efficacy of our method. While in the simulation exercise we can estimate the model to a vector of 100 returns in a matter of seconds, in the empirical application the estimation converges in about 10 minutes.
Cristian M\"ucher acknowledges financial support from the Graduate School of Decision Sciences (GSDS), University of Konstanz, Germany and from the German federal state of Baden-W{\"u}rttemberg through a Landesgraduiertenstipendium. Roxana Halbleib acknowledges financial support from the German Science Foundation through the project HA 8672/1. We want to thank Robin Braun, Ralf Br\"uggemann, Maurizio Daniele, Giampiero Gallo, Lyudmila Grigoryeva, Julie Schnaitmann, Winfried Pohlmeier, the participants of the SoFiE Conference 2019, QFFE International Conference 2019, CFE 2018 and Econometrics Colloquium at the University of Konstanz, for helpful comments. All remaining errors are ours. We acknowledge computational support by the state of Baden-W\"urttemberg through bwHPC.