Estimation of the Cholesky Multivariate Stochastic Volatility Model Using Iterated Filtering

Aim: The paper aims to propose a new estimation method for the Cholesky Multivariate Stochastic Volatility Model based on the iterated filtering algorithm (Ionides et al., 2006, 2015). Methodology: The iterated filtering method is a frequentist-based technique that through multiple repetitions of the filtering process, provides a sequence of iteratively updated parameter estimates that converge towards the maximum likelihood estimate. Results: The effectiveness of the proposed estimation method was shown in an empirical example in which the Cholesky Multivariate Stochastic Volatility Model was used in a study on safe-haven assets of one market index: Standard and Poor’s 500 and three safe-haven candidates: gold, Bitcoin and Ethereum. Implications and recommendations: In further research, the iterating filtering method may be used for more advanced multivariate stochastic volatility models that take into account, for example, the leverage effect (as in Ishihara et al., 2016) and heavy-tailed errors (as in Ishihara and Omori, 2012). Originality/Value: The main contribution of the paper is the proposition of a new estimation method for the Cholesky Multivariate Stochastic Volatility Model based on iterated filtering algorithm This is one of the few frequentist-based statistical inference methods for multivariate stochastic volatility models.


Introduction
Understanding the correlation structure in asset allocations is crucial for effective portfolio construction and risk management as it allows investors to make decisions that align with their financial goals and risk tolerance.Firstly, by combining assets with low or negative correlations, investors may smooth out the portfolio's performance and reduce the impact of individual asset price fluctuations on the overall portfolio.Secondly, correlation affects also the volatility of a portfolio: assets with high positive correlations lead to higher portfolio volatility, and conversely, assets with low or negative correlations reduce overall portfolio volatility.Thirdly, knowing the correlation structure allows investors to perform better stress testing and scenario analysis to assess portfolio performance under different market conditions.Finally, certain asset allocation strategies and complex financial products, such as derivatives, rely heavily on understanding the correlation structure.
The volatility of asset prices can be highly correlated, and individual asset volatilities can exhibit complex and time-varying patterns.Over the last several decades, significant progress has been made in modelling the volatilities of multivariate stock market returns.Nowadays, there are four main approaches: multivariate generalised autoregressive conditional heteroskedasticity (MGARCH) models (Engle, 2002;Engle and Kroner, 1995), multivariate stochastic volatility (MSV) models (Chib et al., 2009) realised covariance models (Bollerslev et al., 2018;Jin and Maheu, 2013), and machine learning (ML) algorithms (Bejger and Fiszeder, 2021;Fiszeder and Orzeszko, 2021).All of these approaches have some strengths and weaknesses, and the choice between them depends on the specific characteristics of the data, the objectives of the analysis, and the trade-offs between modelling flexibility and computational complexity.
Multivariate stochastic volatility models are a class of models developed to describe the joint dynamics of volatility for multiple assets or financial instruments.The crucial feature of MSV models is the assumption that for each asset, the volatility is assumed to follow a stochastic process.Consequently, the conditional covariance matrix is always time-varying.However, the corresponding correlation matrix may be constant or dynamic.Both approaches are used in finance to capture the relations between asset returns, but differ in their assumptions and flexibility.In a constant correlation model, it is assumed that the correlation between the asset returns remains constant over time: the correlation matrix of asset returns is fixed and does not change with market conditions.This simplifies the modelling process and reduces computational complexity since the correlation matrix is estimated once and applied uniformly across all time periods.In contrast, dynamic correlation models allow for the correlation between asset returns to vary over time.These models try to capture changing relations and fluctuations during different market conditions.Dynamic correlation models are more complex to estimate and require more parameters to capture the time-varying correlations accurately and are more computationally demanding.
There are several approaches to incorporating dynamic correlation into multivariate stochastic volatility models.Yu and Meyer (2006) proposed the time-varying correlation bivariate model based on the Fisher transformation.However, their proposition has an obvious drawback, as it cannot be generalised for higher dimensions.Philipov andGlickman (2006a, 2006b) obtained a time-varying correlation matrix assuming that the conditional covariance matrix follows an inverted Wishart distribution.Asai et al. (2006) proposed a dynamic correlation model based on matrix exponential transformation which was further developed by Ishihara et al. (2016).In this paper, the author used the dynamic correlation multivariate stochastic volatility model based on the Cholesky decomposition (Cholesky Multivariate Stochastic Volatility Model, ChMSV) introduced by Tsay (2005) and detailed by Asai et al. (2006).This model is a versatile tool that can be used in various areas of financial research, especially when the main task is to examine correlations over time.The ChMSV model has been used in a number of empirical studies on economic time series, including: investigation of market return predictability on the basis of three predictor variables: dividend yield, consumption-wealth ratio and bond-yield (Lopes et al., 2011), examination of time variation in US monetary policy (Hartwig, 2019), two analyses of stock market co-movements: between three stock indices (Warsaw Stock Exchange Index WIG, the Standard and Poor's S&P500, The Financial Times Stock Exchange FTSE 100 Index) in Pajor (2010), and between five stock indices (US Dow Jones Industrial, the German DAX 30 Performance, the European EuroStoxx50, the Japanese Nikkei 225, and the Chinese Shanghai Shenzen CSI 300) in Zaharieva et al. (2020).
The main contribution of the paper is the proposed estimation technique of ChMSV using iterated filtering algorithm (Ionides et al., 2006(Ionides et al., , 2015)).Iterated filtering allows the estimation of parameters for general non-linear or non-Gaussian State-Space models (SSM) based on the frequentist approach.This makes the method particularly attractive in the context of stochastic volatility models for which the Bayesian approach dominates.Frequentist-based statistical inference for multivariate stochastic volatility is still limited.Harvey et al. (1994) proposed to use the quasi-maximum likelihood method, but their proposal was restricted only to a simple basic stochastic model with a constant correlation matrix.Jungbacker and Koopman (2006) introduced importance sampling Monte Carlo techniques for maximum likelihood estimation in three specific multivariate extensions of the basic stochastic volatility (SV) model.
The paper is organized as follows: Section 2 presents the Cholesky Multivariate Stochastic Volatility Model, Section 3 outlines the estimation methodology, Section 4 provides an empirical example, Section 5 conducts a simulation study, and finally, Section 6 concludes the paper.

The Cholesky Multivariate Stochastic Volatility Model
The main challenge in constructing multivariate stochastic volatility models with a dynamic correlation matrix is to ensure positive definiteness of the conditional covariance matrix.For this purpose, Tsay (2005) proposed to use the Cholesky decomposition.In fact, Tsay used the LDL variant of this method, which allows to decompose positive definite  to the form  = ′, where  is lower triangular matrix with unit diagonal elements, and  is diagonal matrix with positive elements.This transformation allows to separate elements related to variance from those related to correlation.Asai et al. (2006) suggested to use the autoregressive process of order 1 to specified elements of both matrix  and .This study followed Pajor (2010, p. 82) and Chib et al. (2009), and used instead the mean-reverted process of order 1.Mean-reversion of volatility is a commonly observed phenomenon, and it seems that this assumption can also be applied to correlations.The Cholesky Multivariate Stochastic Volatility Model employed in this work takes the form (1) Moreover, it was assumed (Pajor, 2010, p. 82) that: and   ,   are statistically independent for ,  ϵ ℤ.
In the ChMSV the conditional distribution of vector of returns given the past information ℱ −1 is p-dimensional normal distribution where Σ  =        is a time-variant covariance matrix.The elements of this matrix can be determined according to the following formulas (Asai et al., 2006): a) diagonal elements: b) non-diagonal elements: Consequently, both conditional variances and covariances are time-varying, but their dynamics are not determined separately, because both are dependent on ℎ , and  , processes.The standard deviation of i-th asset return and correlation between i-th and j-th assets return may be calculated correspondingly as and The number of parameters in the ChMSV model equals 3( 2 + )/2.For a 2, 3 and 4-dimensional model, respectively, the study obtained 9, 18, and 30 parameters.A limitation of this specification is that it is not straightforward to interpret the parameters of ℎ , and  , processes, because they affect both variances and covariances.

State Space Model Representation of the ChMSV Model
The research goal was to estimate vector of parameters of the ChMSV model using iterated filtering.Firstly, it had to be shown that the ChMSV model is in fact element of a broader class of statistical models called State Space Models (SSM).This class is a statistical framework used to describe and analyse complex systems that evolve over time, and allow to estimate the hidden states based on the observed data (filtering), predict future states of latent process and perform parameter estimation.They consist of two components: a hidden or unobserved state process   = (  ) ≥0 and an observed process   = (  ) ≥0 .The latent and observed processes are specified respectively by the set of densities �  (  | −1 ; )� ≥1 (with the initial density  0 ( 0 ; )) and �  (  |  ; )� ≥0 .Latent process   is Markovian and measurement process   is assumed to be conditionally independent given   (for details of this and two alternative definitions see Chapter 2 in Chopin and Papaspiliopoulos (2020)).
As a result, the ChMSV model may be considered as a example of a State Space Model.This makes it possible to use the filtration, prediction and estimation of parameters as in other SSM models.However, the ChMSV model is nonlinear, precluding the use of the Kalman filter.Nowadays, particle filters are a well-established technique for estimating the latent state in nonlinear or non-Gaussian State Space Models with fixed parameters.However, a direct application of particle filters for likelihood estimation was not possible due to a limitation: the likelihood estimator was non-continuous as a function of the parameters (Malik and Pitt, 2011).A number of estimation methods indirectly using particle filters have been proposed: approximation of gradient ascent (Liu and West, 2001), expectation-maximisation algorithms, Particle Markov Chain Monte Carlo methods (Andrieu et al., 2010), SMC2 (Chopin et al., 2013).A detailed discussion of this issue can be found in Kantas et al. (2015) and in Chapter 14 in Chopin and Papaspiliopoulos (2020).In this study the author used the iterated filtering algorithm proposed by Ionides et al. (2015).

Iterated Filtering
Iterated filtering, initially proposed by Ionides et al. (2006) was later supported by theoretical details in Ionides et al. (2011), and underwent further development with the introduction of the second generation, IF2, initiated by Ionides et al. (2015).Although both generations of iterated filtering utilise a recursive filtering approach with the augmented model, their theoretical foundations differ.The first generation, IF1, relies on approximating the Fisher score function, whereas the second generation, IF2, integrates the idea of data cloning (proposed by Lele et al., 2007) with the convergence of an iterated Bayes map, as presented by Nguyen (2016).Empirical investigations demonstrated that IF2 surpasses IF1 in performance (Ionides et al., 2015).Consequently, the calculations presented in this article used the second generation of the algorithm.
Iterated filtering uses a basic bootstrap particle filter (Gordon et al., 1993) and, therefore, avoids the need to directly evaluate the transition density   (  | −1 ; ).Instead, it only requires the capability to simulate from this density (simulation-based approach).However, iterated filtering employs particle filters not directly in a target model, but in a similar model with parameters that take a random walk in time.This added variability smooths the likelihood surface, which addresses the main challenge faced by particle filters in parameter estimation and mitigates the issue of particle depletion.Through multiple repetitions of the filtering process, each using a particle filter, the variance of the random walk tends to zero, and the augmented model progressively converges to the original one.Consequently, iterated filtering provides a sequence of iteratively updated parameter estimates that converge towards the maximum likelihood estimate (likelihood approach), as detailed in Ionides et al. (2015).Unlike other simulation-based methods such as Approximate Bayesian Computation (Toni et al., 2009), Particle Markov Chain Monte Carlo (Andrieu et al., 2010), and SMC2 (Chopin et al., 2013) based on the Bayesian approach, iterated filtering relies on the frequency interpretation of probability (frequentist-based approach).As indicated by King et al. (2016), iterated filtering is one of the few (and possibly the only) method that combines these three approaches and applies to most SSM models.Iterated filtering has demonstrated successful applications to various State Space Models, predominantly within the context of epidemiology (Bhadra et al., 2011;He et al., 2009;King et al., 2008;Stocks et al., 2020;You et al., 2020).Additionally, it has been employed in economic modelling, particularly for univariate (Bretó, 2014;Szczepocki, 2020) and bivariate (Szczepocki, 2022) stochastic volatility models.

Iterated Filtering Setup of the ChMSV Model
The main advantage of iterated filtering is that the iteration algorithm, to apply to specific SSMs, requires only three conditions: 1) ability to simulate from initial density  0 ( 0 ; ), 2) ability to simulate from transition density   (  | −1 ; ), 3) evaluate measurement density   (  |  ; ).
The study showed that in the case of the ChMSV model, all the conditions were fulfilled.The initial assumed conditions are with probability 1, which in practical implementations are treated as fixed.Simulating from the transition density requires drawing from a p-dimensional normal distribution.However, a diagonal form of the covariance matrix allows drawing from one-dimensional normal distributions according to formulas: and Moreover, evaluating the measurement density   (  |  ; ) is straightforward, because it requires the determination of the density value of the -dimensional normal distribution with zero mean vector and Σ  covariance matrix with elements defined by equations ( 4) and ( 5).
The calculations for the empirical part and the simulation study in the article were carried out using the POMP package (Partially Observed Markov Processes, King et al., 2016) designed specifically for the R statistical computing environment (R Development Core Team, 2010), hence the author used the random-walk perturbations of the parameters incorporated in the POMP package.As a sequence of the decay of perturbations, a geometric sequence was chosen.In order to increase computational efficiency, the code responsible for handling the initial function, transition, and measurement densities was implemented in the C programming language.

Empirical Example
This section applied the author's estimation method of the Cholesky Multivariate Stochastic Volatility Model to the returns of four assets: Standard and Poor's 500 index (SPX), the price of gold in US dollars, and the prices in US dollars of two cryptocurrencies, Bitcoin and Ethereum.The sample period from 4 June 2018 to 4 June 2023, yielded a total of 1258 observations of daily logarithmic returns multiplied by 100, only keeping the days when observations were present for all time series.The data came from stooq.pl database.The same set of variables was part of the empirical study of safe haven assets in Będowska-Sójka and Kliber (2021), which, however, used the data to estimate three bivariate Yu and Meyer stochastic volatility models, each time pairing the SPX index with a different safe haven candidate: gold, Bitcoin and Ethereum.In this example, these four assets are analysed together as a four-dimensional model of stochastic volatility which allows to estimate not only the dynamic correlation between the index and the candidate for a safe haven asset, but also to determine the correlation between the candidates.
The literature on safe assets is extensive.One of the key works is was Baur and Lucey (2009), in which a precise conceptual differentiation between two financial terms 'safe haven' and 'hedge' was introduced: "a safe haven is defined as an asset that is uncorrelated or negatively correlated with another asset or portfolio in times of market stress or turmoil, whereas a hedge is a security that is uncorrelated with the stock market on average".Baur and McDermott (2010) further expanded on this concept by introducing the distinction between strong and weak safe haven effects: "a strong (weak) safe haven is defined as an asset that is negatively correlated (uncorrelated) with another asset or portfolio in certain periods only, e.g. in times of falling stock markets".Therefore, to verify whether a given asset is a safe haven, it is not enough to calculate the sample Pearson correlation coefficient over a certain period.On the contrary, to determine the correlation in moments of turbulence in the stock markets, it is necessary to estimate the time-varying correlation.Many different methods have been used to estimate dynamic correlation, including regression models with a heteroscedastic random error (Baur and McDermott, 2010), bivariate DCC GARCH models (Aftab et al., 2019;Bouri et al., 2017) and bivariate stochastic volatility models (Będowska-Sójka and Kliber, 2021;Kliber et al., 2019).This empirical example only illustrates that the ChMSV model also may be used to estimate the time-varying correlation in the case of safe haven asset studies, and is not an attempt to determine whether these assets actually play such a role.
Figure 1 presents time series of asset prices (left column) and returns (right column).Table 1 shows summary statistics of returns.Among the analysed assets, Ethereum was the most volatile, followed by Bitcoin.Standard and Poor's 500 index, and gold were less volatile.Cryptocurrencies also had the lowest minimum values and the highest maximum values, however the kurtosis of Standard and Poor's 500 index was higher than the kurtosis of the cryptocurrencies.Such a high kurtosis value combined with a more moderate standard deviation may indicate that Standard and Poor's 500 index is less volatile than cryptocurrencies most of the time yet with exceptionally strong extreme returns compared to typical values.Figure 1 shows for the index that the most extreme values are related to the outbreak of the pandemic, while cryptocurrencies volatility is more even over time.It is also significant that gold is the least volatile, has the lowest kurtosis and has the least extreme values.
The estimation results are presented in Table 2.As mentioned in Section 2, interpretation of the obtained estimates is not possible.Additionally, the results cannot be compared with one-dimensional stochastic volatility models.Algorithm settings: 200 iterations, 1,000 particles, random-walk perturbations with the initial 0.01 perturbations, and a geometric decay of perturbations of α = 0.5.Standard errors of the parameters were calculated via numerical approximation to the Hessian (see supporting text by Ionides et al. (2006) for details of the procedure).The log-likelihood estimation with the standard errors were obtained by averaging twelve likelihood evaluations of a bootstrap particle filter with 1,000 particles.
Source: own work.
Having the model parameter estimates, one can apply a bootstrap particle filter to estimate the standard deviation of the returns and the dynamic correlation.This is possible by presenting the model in the form of a state space model and using formulas ( 6) and (7). Figure 2 presents estimates of the standard deviation of returns for the analysed time series, which shows that the estimates reflect well the volatility of returns observed in the right column of Figure 1. Figure 2 also shows a huge rise in volatility at the beginning of the pandemic for Standard and Poor's 500 index and for gold, increased volatility of gold in 2016, Bitcoin at the turn of 2017-2018 and Ethereum at the beginning of 2016.
One of the most important advantages of the ChMSV model is the ability to estimate time-varying (dynamic) correlations.Figure 3 presents the correlation between Standard and Poor's 500 index and three candidates for a safe haven asset.One can see that all assets are slightly positively correlated with the index for most of the sample.The correlations are stable and during the index's greatest decline at the beginning of the pandemic, the correlations became reached their peak.It is clear from Figure 3 that none of the three assets can be treated as a strong safe haven asset because they are non-negatively correlated throughout the period.A subject of further analysis is the examination of whether they can be treated as a weak safe haven asset.
Furthermore, the four-dimensional model allowed to study the correlations between safe haven candidates.This can be useful in finding more safe haven assets so that they would be negatively correlated or at least uncorrelated with each other.Figure 4 shows that gold is moderately correlated with both cryptocurrencies for most of the sample, but at selected moments the direction of the correlation turned negative for a short time.Bitcoin's correlation with Ethereum is very strong, but also very volatile.The results were obtained by using a particle filter with 1000 particles Source: own work.

Simulation Study
Due to the inability to analytically demonstrate the convergence of the iterated filtering algorithm to the maximum likelihood estimates, the author evaluated its performance using a simple simulation study.The objective was to assess the precision of parameter estimation, conducting 100 simulations of the ChMSV model.To ensure consistency, each simulation had the same dimensions as the empirical example (four time series x 1258 observations) and the same parameter values obtained from the empirical study outlined in Section 4 (see Table 2).The same estimation procedure was applied as in the empirical study for each simulation.The findings, summarised in Table 3, present the mean errors (ME) and the root mean square errors (RMSE) of the parameter estimates.It can be seen from the ME that for any parameter there is no evidence that the errors tend to be overestimated in one direction.
The biggest errors are related to the long-term mean parameter, perhaps because this one parameter has unlimited support.The RMSE values are comparable to the estimates of standard errors obtained in an empirical study.The results indicate that the proposed method is sufficiently reliable.Algorithms settings: 200 iterations, 1000 particles, the random-walk perturbations with initial 0.01 perturbations, and geometric decay of perturbations of α = 0.5 for all parameters.
Source: own work.

Conclusion
The main novelty of this paper is the proposition of iterated filtering algorithm for estimating the Cholesky Multivariate Stochastic Volatility Model.The author proved that this model is an element of a broader class of statistical models called State Space Models and therefore the techniques used for this class (among others, iterated filtering) can also be applied in this case.Iterated filtering provides a sequence of iteratively updated parameter estimates that converge towards the maximum likelihood estimate.This is one of the few frequentist-based statistical inference methods for multivariate stochastic volatility models.Additionally, a bootstrap particle filter may be used for estimating the filtering distribution of log-volatilities and dynamic correlation in the frequentist-based approach.The effectiveness of the proposed estimation method was shown in an empirical example in which the ChMSV model was used in a study on safe-haven assets of one market index: Standard and Poor's 500 and three safe-haven candidates: gold, Bitcoin and Ethereum.The author presented parameter estimates and estimates of standard deviation and dynamic correlations, and the parameter estimates obtained from the empirical study were further used for the simulation experiment.It was confirmed that if the form of the model is correct, the proposed estimation method reproduces well the parameter values.
In further research, it would be interesting to try to apply the iterating filtering to more advanced multivariate stochastic volatility models that take into account, for example, the leverage effect (as in Ishihara et al. 2016) and heavy-tailed errors (as in Ishihara and Omori, 2012).The leverage effect in stochastic volatility modelling refers to the phenomenon where the volatility of an asset is negatively correlated with its return.This results in the tendency for the volatility of an asset to increase when the asset's returns are negative, and decrease when the returns are positive.In contrast, the use of heavy-tailed errors in multivariate stochastic volatility models is motivated by the desire to capture the presence of extreme events in financial time series.The heavy-tailed distribution allows for a more realistic representation of the fat tails observed in financial returns, meaning that extreme events are more probable than would be implied by a normal distribution.

Fig. 1 .
Fig. 1.Time series of Standard and Poor's 500 index (SPX) and the price of gold, Bitcoin, Ethereum (left column) and corresponding logarithmic returns multiplied by 100 (right column) in the period from 4 June 2018 to 4 June 2023 Source: own work.

Fig. 2 .
Fig. 2. Estimates of the standard deviation of returns for four time series: Standard and Poor's 500 index (SPX) and the price of gold, Bitcoin, and Ethereum.The sample period: from 4 June 2018 to 4 June 2023.The results were obtained by using a particle filter with 1000 particles Source: own work.

Fig. 3 .
Fig. 3. Estimates of dynamic correlation calculated for returns of Standard and Poor's 500 index (SPX) and returns of gold, Bitcoin, and Ethereum in the period from 4 June 2018 to 4 June 2023.The results were obtained by using a particle filter with 1000 particles Source: own work.

Fig. 4 .
Fig. 4. Estimates of dynamic correlation calculated for three pairs of returns: gold and Bitcoin (top), gold and Ethereum (middle), and Ethereum and Bitcoin (bottom), in the period from 4 June 2018 to 4 June 2023.The results were obtained by using a particle filter with 1000 particles

Table 1 .
Summary statistics of logarithmic returns multiplied by 100 calculated for four time series of Standard and Poor's 500 index (SPX) and the price of gold, Bitcoin, and Ethereum

Table 2 .
Parameter and log-likelihood estimates and standard errors (in brackets) of the Cholesky Stochastic Volatility model were obtained using iterated filtering for the analyzed data

Table 3 .
The results of the simulation study based on 100 time series simulation of the ChMSV model with the same parameter values as those obtained in the empirical study and the same length of time series as the analysed of S&P500 index and gold