Approximating Correlation Matrices Using Stochastic Lie Group Methods

Abstract

Specifying time-dependent correlation matrices is a problem that occurs in several important areas of finance and risk management. The goal of this work is to tackle this problem by applying techniques of geometric integration in financial mathematics, i.e., to combine two fields of numerical mathematics that have not been studied yet jointly. Based on isospectral flows we create valid time-dependent correlation matrices, so called correlation flows, by solving a stochastic differential equation (SDE) that evolves in the special orthogonal group. Since the geometric structure of the special orthogonal group needs to be preserved we use stochastic Lie group integrators to solve this SDE. An application example is presented to illustrate this novel methodology.

1. Introduction

Correlation matrices play an important role e.g., in finance and risk management. A valid correlation matrix is a real matrix with the following properties:

1.

All diagonal elements of a correlation matrix are equal to one and absolute values of all non-diagonal elements are less than or equal to one.

2.

Correlation matrices are real symmetric and positive semi-definite, i.e. all eigenvalues are non-negative.

In this paper our goal is to construct time-dependent correlation matrices that fulfill the properties above and approximate the true correlation using real market data.

There are already methods available that were designed to tackle the same problem, see e.g., [1,2,3,4]. Newton-based methods for approximating covariance matrices can be found in [5,6,7]. Furthermore, there exist methods using hyperspherical decomposition [8] and unconstrained convex optimization [9]. But some methods show weaknesses in at least one of the desired properties of a correlation matrix mentioned above. Especially, the positive semi-definiteness is a criteria which is not well implemented. For example, the approach of [10] suffers from drawbacks in other portions of the matrix in order to maintain positive semi-definiteness.

Here we ensure the positive semi-definiteness of the correlation matrices constructed with our methodology by taking up the idea of [11]. The authors defined covariance flows based on isospectral flows by constructing matrices similar to an initial valid covariance matrix. This is a well-analyzed approach but still it lacks the stochastic nature of correlations.

In our methodology we include the stochastic behaviour of correlations by assuming that the orthogonal matrices needed for the covariance flows are driven by a stochastic differential equation (SDE) rather than an ordinary differential equation (ODE). Since the space of orthogonal matrices can be viewed as a Lie group, we use Lie group integrators [12] to solve this SDE. Lie group methods arose in the deterministic case for solving ODEs such that geometric properties of the Lie group are preserved, see e.g., [13,14] and ([15], p. 126). Their application in the stochastic setting has been studied in e.g., [16,17,18]. So far, stochastic Lie group methods were mostly applied to SDEs considered in information theory [19] and engineering [20]. For further concise literature on Lie group methods, we refer the interested reader e.g., to ([15], Chapter IV.8) or [13].

However, the application of Lie group methods on SDEs considered in finance has not been analyzed yet. Consequently, the contribution of our paper is twofold: we respect the stochastic behaviour of correlations by considering an SDE that evolves in a Lie group and we use a stochastic Lie group integrator in a financial mathematics context to get a numerical solution of the considered SDE.

Our methodology can be used e.g., if one wants to specify time-dependent correlation matrices in a given time interval and is only aware of the correlation matrix at the initial time point. In the following sections we aim to formulate an SDE that describes the correlation flows for the given time interval. Moreover, we approximate the solution of this SDE by applying the geometric Euler-Maruyama scheme presented in [18,21]. We demonstrate how techniques from geometric numerical integration which arose from a mechanical engineering background can be used in a financial mathematical setting.

The remainder of the paper is organized as follows. Covariance flows are introduced in Section 2. In Section 3 we present a numerical method to solve SDEs based on the relation between Lie group and Lie algebra. We then turn our attention to simulations and the application of our methodology in risk management by using real historical market data in Section 4. A discussion of our results is given in Section 5.

2. Covariance Flows

For creating valid time-dependent correlation matrices, we first introduce covariance flows $\{P_t:0\le t\le T\}$. A covariance flow is a set of similar, time-dependent matrices $P_t$ of the covariance space

$$P\left(n\right)=\{P\in {\mathbb{R}}^{n\times n}:P=P^{\top },x^{\top }Px\ge 0\mathrm{for}\mathrm{all}x\in {\mathbb{R}}^n,n\in \mathbb{N}\}.$$

(1)

These covariance matrices can be easily transformed into corresponding correlation matrices $R_t$. The computation of covariance flows is based on isospectral flows, see [11].

Due to the symmetry of covariance matrices the principal axis theorem can be applied, i.e., there exists an orthogonal matrix Q and a diagonal matrix D consisting of the eigenvalues of P such that

$$P=Q^{\top }DQ.$$

(2)

Without loss of generality, we can assume that Q is a rotation matrix whose determinant is always equal to 1. Thus, we assume Q to be an element of the special orthogonal group

$$\mathrm{SO}\left(n\right)=\{X\in \mathrm{GL}\left(n\right):X^{\top }X=I,det\left(X\right)=1\},$$

(3)

where GL(n) is the general linear group, i.e., the set of $n\times n$ invertible matrices. Since SO(n) defines a differentiable manifold and the matrix multiplication is a differentiable mapping, SO(n) is a matrix Lie group. The corresponding Lie algebra is denoted by $\mathfrak{so}\left(n\right)$ and consists of skew-symmetric $n\times n$-matrices, for details we refer to e.g., ([22], p. 58).

Now let an initial covariance matrix $P_0$ be given. We consider the covariance flow

$$P_t=Q_t^{\top }{P}_0{Q}_t.$$

(4)

Creating time-dependent covariance matrices $P_t$ that are similar to $P_0$ implies the generation of time-dependent, orthogonal matrices $Q_t$, $t\ge 0$. Therefore, we consider the following stochastic differential equation (SDE):

$$\mathrm{d}{Q}_t=Q_t{K}_t\mathrm{d}t+Q_t\left(V_{t,1}\mathrm{d}{W}_{t,1}+V_{t,2}\mathrm{d}{W}_{t,2}\right),Q_0=I,$$

(5)

where $K_t,V_{t,1},V_{t,2}\in {\mathbb{R}}^{n\times n}$ and $W_{t,1}$ and $W_{t,2}$ are independent, standard Brownian motions, i.e. it holds $\mathrm{d}{W}_{t,1},\mathrm{d}{W}_{t,2}\sim \mathcal{N}(0,\mathrm{d}t)$. In order to ensure that the resulting matrices $Q_t$ have the desired properties, we use the characterization of $K_t$, $V_{t,1}$ and $V_{t,2}$ defined in the following theorem.

Theorem 1.

If $Q_t$ obeys the matrix SDE (5) then $Q_t\in$ SO(n) if and only if $V_{t,1},V_{t,2}\in \mathfrak{so}\left(n\right)$ and $K_t+K_t^{\top }=V_{t,1}^2+V_{t,2}^2$.

A more general version and the proof of this theorem can be found in [18]. Theorem 1 suggests that the easiest way to determine the unknown matrices is to first find symmetric matrices $V_{t,1}^2$ and $V_{t,2}^2$ such that skew-symmetric square roots $V_{t,1}$ and $V_{t,2}$ exist. The matrix $K_t$ can then be identified as the lower triangular matrix of $Y_t^2:=V_{t,1}^2+V_{t,2}^2$ where the entries on the diagonal of $K_t$ are equal to 0.5 times the diagonal elements of $Y_t^2$. The matrices $V_{t,1},V_{t,2}$ and $K_t$ can be fixed according to a given problem.

Regarding our goal, which is the approximation of correlation matrices using real market data, we estimate the matrices $V_{t,1},V_{t,2}$ and $K_t$ from the considered market data by using a least-squares approach. We describe our initialization of these matrices more detailed in Section 4.

3. Stochastic Lie Group Method

In the following we are taking a closer look on the SDE (5) and how it can be solved numerically. In general, there is no closed form solution. However, a solution can be defined via a Magnus expansion $Q_t=Q_0{e}^{{\mathsf{\Omega }}_t}$, ${\mathsf{\Omega }}_0=0$, where ${\mathsf{\Omega }}_t\in {\mathbb{R}}^{n\times n}$ obeys a matrix SDE. This auxiliary SDE is given in the following theorem and can also be found in [18].

Theorem 2.

If $Q_t=Q_0{e}^{{\mathsf{\Omega }}_t}$, ${\mathsf{\Omega }}_0=0$, obeys the SDE (5) then ${\mathsf{\Omega }}_t$ obeys the SDE

$$\mathrm{d}{\mathsf{\Omega }}_t=A_t\mathrm{d}t+{\Gamma }_{t,1}\mathrm{d}{W}_{t,1}+{\Gamma }_{t,2}\mathrm{d}{W}_{t,2},$$

(6)

where

$$\begin{array}{cc}\hfill A_t& =\mathrm{d}{exp}_{-{\mathsf{\Omega }}_t}^{-1}\left(A_{t,0}\right)-\frac{1}{2}\mathrm{d}{exp}_{-{\mathsf{\Omega }}_t}^{-1}\left(C_{t,1}\right)-\frac{1}{2}\mathrm{d}{exp}_{-{\mathsf{\Omega }}_t}^{-1}\left(C_{t,2}\right),\hfill \\ \hfill A_{t,0}& =K_t-\frac{1}{2}{V}_{t,1}^2-\frac{1}{2}{V}_{t,2}^2,\hfill \\ \hfill C_{t,r}& =\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }\frac{1}{(p+q+2)}\frac{{(-1)}^p}{p!(q+1)!}{\mathrm{ad}}_{{\mathsf{\Omega }}_t}^p\left({\mathrm{ad}}_{{\Gamma }_{t,r}}({\mathrm{ad}}_{{\mathsf{\Omega }}_t}^q\left({\Gamma }_{t,r}\right))\right),r=1,2,\hfill \\ \hfill {\Gamma }_{t,r}& =\mathrm{d}{exp}_{-{\mathsf{\Omega }}_t}^{-1}\left(V_{t,r}\right),r=1,2.\hfill \end{array}$$

The expression $\mathrm{d}{exp}_{-\mathsf{\Omega }}^{-1}\left(X\right)$ is given by

$$\mathrm{d}{exp}_{-\mathsf{\Omega }}^{-1}\left(X\right)=\sum _{k=0}^{\infty }\frac{B_k}{k!}{\mathrm{ad}}_{-\mathsf{\Omega }}^k\left(X\right),$$

where $B_k$ are the Bernoulli numbers and ${\mathrm{ad}}_X\left(Y\right)=[X,Y]$ is the adjoint operator which is used iteratively

$$\begin{array}{cc}\hfill {\mathrm{ad}}_X^0\left(Y\right)& =Y,\hfill \\ \hfill {\mathrm{ad}}_X^i\left(Y\right)& =\left[X,{\mathrm{ad}}_X^{i-1}\left(Y\right)\right]={\mathrm{ad}}_X\left({\mathrm{ad}}_X^{i-1}\left(Y\right)\right),i\ge 1.\hfill \end{array}$$

(7)

In the case where the solution $Q_t$ of the SDE (5) is in the Lie group SO(n), it holds that the matrix SDE (6) evolves in the corresponding Lie algebra $\mathfrak{so}\left(n\right)$ [18]. The Lie algebra $\mathfrak{so}\left(n\right)$ is the tangent space of the differentiable manifold SO(n) at the identity, see e.g., ([22] p. 71). This simple Euclidean-like geometry of the Lie algebra can be used to solve SDEs that evolve in the associated Lie group.

Therefore, we compute a numerical solution of the SDE (5) by solving the SDE (6) e.g., with the Euler-Maruyama scheme in the Lie algebra $\mathfrak{so}\left(n\right)$ and projecting this solution ${\mathsf{\Omega }}_t$ back onto the Lie group via the exponential map, $exp:$$\mathfrak{so}\left(n\right)\to \mathrm{SO}\left(n\right)$, to get a solution for $Q_t$. Since this scheme preserves the geometry of the manifold SO(n), it is called the geometric Euler-Maruyama scheme [18]. One can easily check that the geometry of SO(n), namely the condition $Q_t^{\top }{Q}_t=I$, is not preserved if a numerical integration scheme is applied directly to (5) instead of (6). A simple version of the geometric Euler-Maruyama scheme applied to (5) is given in the following algorithm. It can be interpreted as a special case of a stochastic version of the Runge-Kutta–Munthe-Kaas schemes [14] and a more general version of the algorithm can also be found in [18].

There are only limited convergence results available for stochastic Lie group methods. However, it was proven in [21] that the geometric Euler-Maruyama scheme converges with rate $\mathcal{O}\left({\Delta }_j\right)$ with respect to mean uniform squared error over the whole interval $[0,T]$.

In the Euler-Maruyama step of Algorithm 1 we truncated each of the infinite sums involved in (6) to the first summand only. The geometry is preserved under these truncations because for skew-symmetric matrices X and Y, the adjoint operator ${\mathrm{ad}}_X^k\left(Y\right)$ is also skew-symmetric for any $k\ge 0$, which can be easily proved by induction. It follows that any truncation of ${\mathsf{\Omega }}_{j+1}$ is in $\mathfrak{so}\left(n\right)$ and thus, any projection of ${\mathsf{\Omega }}_{j+1}$ is in SO(n).

In the Projection step of Algorithm 1 and in Theorem 2 the exponential map is used as a parametrization for the Lie group. However, the basic concepts of Lie group methods are not limited to this specific parametrization. One could also use other mappings, e.g., the Cayley transform $\mathrm{cay}\left(\mathsf{\Omega }\right)={(I-\mathsf{\Omega })}^{-1}(I+\mathsf{\Omega })$. Since a truncation of infinite sums induced by the definition of the matrix exponential can be avoided, considering the Cayley map instead might be beneficial in cases of higher dimensions where no closed form expressions for the exponential map are available. A comparison of the usage of these two maps, $exp\left(\mathsf{\Omega }\right)$ and $\mathrm{cay}\left(\mathsf{\Omega }\right)$, in the context of Lie group methods can be viewed in [23].

Algorithm 1: Geometric Euler-Maruyama
	We divide the time interval $[0,T]$ into subintervals $[t_j,t_{j+1}]$, $j\ge 0$. Starting with $t_0=0$ and $Q_0=I$ the following steps are repeated over successive intervals $[t_j,t_{j+1}]$ until $t_{j+1}=T$.
	Initialization step: Let $Q_j$ be the approximation of $Q_t$ at $t=t_j$. Analogously: $K_j$, $V_{j,1}$ and $V_{j,2}$. Euler-Maruyama step: Compute $$\begin{array}{cc}\hfill {\Delta }_j& =t_{j+1}-t_j,\hfill \\ \hfill \Delta W_{j,r}& \sim \mathcal{N}(0,{\Delta }_j)\mathrm{for}r=1,2,\hfill \\ \hfill {\mathsf{\Omega }}_{j+1}& =\left(K_j-\frac{1}{2}{V}_{j,1}^2-\frac{1}{2}{V}_{j,2}^2\right){\Delta }_j+V_{j,1}\Delta W_{j,1}+V_{j,2}\Delta W_{j,2}.\hfill \end{array}$$ Projection step: Set $Q_{j+1}=Q_{j}exp\left({\mathsf{\Omega }}_{j+1}\right)$.

4. Simulation

In this section we want to apply the method described in the previous sections to approximate correlations that can be observed in a real market. For this purpose we consider historical prices of the S&P 500 index and the Euro/US-Dollar exchange rate on a daily basis. We compute moving correlations with a window size of 30 days and obtain correlations from 3 January 2005 to 6 January 2006 (see Figure 1).

Assume the following scenario: A risk manager retrieves from the middle office’s reporting system the initial correlation matrix at $t=0$ (which corresponds to 3 January 2005) of the regarded historical data

$$R_0^{\mathrm{hist}}=\left(\begin{array}{cc}1& -0.0159\\ -0.0159& 1\end{array}\right),$$

(8)

i.e., the Pearson correlation coefficient between the S&P 500 index and the Euro/US-Dollar exchange rate was computed as $\rho =-0.0159$. Due to the stochastic behaviour of correlations and the path shown in Figure 1 being only one of many possible realizations a density plot of the historical correlations is needed. Therefore, we suppose that the risk manager is aware of the density function of the considered correlation and estimate a density function from the historical data using kernel smoothing functions, which is also plotted in Figure 2. For more details on the density estimation see [24]. Then, the goal is to create valid time-dependent correlation matrices that reflect the stochastic nature of correlations while trying to match the density function of the historical data.

4.1. Construction of Covariance and Correlation Flows

Our methodology, which is designed to solve the risk manager’s problem described in the scenario above, for approximating correlation matrices can be summarized by the following steps:

• Find matrices $V_{t,1},V_{t,2}$ and $K_t$ such that the conditions in Theorem 1, namely $V_{t,1},V_{t,2}\in \mathfrak{so}\left(n\right)$ and $K_t+K_t^{\top }=V_{t,1}^2+V_{t,2}^2$, are fulfilled.
• Insert the matrices computed in the previous step into (5) and solve this SDE, i.e.

  $$\mathrm{d}{Q}_t=Q_t{K}_t\mathrm{d}t+Q_t\left(V_{t,1}\mathrm{d}{W}_{t,1}+V_{t,2}\mathrm{d}{W}_{t,2}\right),Q_0=I,$$

  by using Algorithm 1, the geometric Euler-Maruyama scheme.
• Compute for a given initial covariance matrix $P_0$ the covariance flow $P_t=Q_t^{\top }{P}_0{Q}_t$.
• Transform the so computed covariance matrices $P_t$ to corresponding correlation matrices $R_t={\Sigma }_t^{-1}{P}_t{\Sigma }_t^{-1}$ with ${\Sigma }_t={\left(diag\left(P_t\right)\right)}^{\frac{1}{2}}$.

4.1.1. Setting the Coefficient Matrices

For the first step, we construct a symmetric matrix $Y_t^2:=V_{t,1}^2+V_{t,2}^2$ such that skew-symmetric matrices $V_{t,1}$ and $V_{t,2}$ can be derived. According to Rinehart [25], every real square root of a symmetric matrix is similar to a real skew square root of this matrix if it is negative semi-definite with nonzero eigenvalues that have even multiplicities. Regarding this required matrix structure, this step might seem like an obstacle. However, this initialization step is actually the part that gives the degrees of freedom needed to incorporate historical market data. To explain this further, we consider the following generator of $\mathfrak{so}\left(2\right)$:

$$G=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right).$$

(9)

An element of $\mathfrak{so}\left(2\right)$ is then any scalar multiple of the generator:

$$\omega \in \mathbb{R}\Rightarrow \omega G\in \mathfrak{so}\left(2\right).$$

(10)

Together with a function $g\left(t\right):{\mathbb{R}}_{\ge 0}\to \mathbb{R}$, an arbitrary time-dependent skew-symmetric matrix $Y_t$ can be defined as

$$Y_t=g\left(t\right)G.$$

(11)

The function $g\left(t\right)$ can be chosen arbitrarily, e.g., as a cubic polynomial. Experimenting with different functions, we found that the following function performed best given the historical data

$$g\left(t\right)=x_1+x_{2}t+x_3{t}^2+x_4{t}^3+x_{5}tsin\left(x_{6}t\right)+x_{7}tcos\left(x_{8}t\right),$$

(12)

where $x_1,\cdots ,x_8\in \mathbb{R}$ represent parameters that can be associated with possible degrees of freedom. The symmetric matrices $V_{t,1}^2$ and $V_{t,2}^2$ can now be arranged, e.g., in a convex combination

$$V_{t,1}^2=x_9·Y_t^2,V_{t,2}^2=(1-x_9)·Y_t^2$$

(13)

where $x_9\in [0,1]$. The so constructed $V_{t,1}^2$ (resp. $V_{t,2}^2$) fulfills the aforementioned requirements such that skew square roots $V_{t,1}$ (resp. $V_{t,2}$) can be computed. As mentioned before we set $K_t$ equal to the lower triangular matrix of $Y_t^2$ where the diagonal elements of $K_t$ are equal to the diagonal entries of $Y_t^2$ multiplied with 0.5.

4.1.2. Preparation for the Geometric Euler-Maruyama Scheme

For the second step, solving (5) with Algorithm 1, we chose to simulate 100 independent realizations of Brownian motions. As a result, we obtained $M=100$ paths of skew-symmetric matrices ${\mathsf{\Omega }}_t$ at each time step $t=t_{j+1}$, $j\ge 0$. Next, we computed the mean of these different paths as an estimator for the expectation value which we then projected onto SO(n) via the exponential map, i.e.

$$\widehat{E}\left({\mathsf{\Omega }}_t\right)=\frac{1}{M}\sum _{k=1}^M{\mathsf{\Omega }}_t^{\left(k\right)},Q_t=Q_{0}exp\left(\widehat{E}\left({\mathsf{\Omega }}_t\right)\right).$$

(14)

In the projection step of Algorithm 1 we used

$$exp\left(\mathsf{\Omega }\right)=exp\left(\begin{array}{cc}0& -\omega \\ \omega & 0\end{array}\right)=\left(\begin{array}{cc}cos\omega & -sin\omega \\ sin\omega & cos\omega \end{array}\right)$$

(15)

for fast computation which results from Taylor series expansion of $sin\omega$ and $cos\omega$. If more than two correlations are considered, e.g., $\mathsf{\Omega }\in \mathfrak{so}\left(3\right)$, then the Rodrigues formula ([26], p. 261) can be used to avoid dealing with the infinite sum expression of the matrix exponential.

4.1.3. Computation of Covariance Flows

For the third step of computing covariance flows, we first need an initial $P_0$ to begin with. As in [11], we estimate the covariance matrix of the whole historical data by using the MATLAB function cov and as a result we got

$$\widehat{C}=\left(\begin{array}{cc}0.2327\times {10}^{-4}& 0.0056\times {10}^{-4}\\ 0.0056\times {10}^{-4}& 0.4269\times {10}^{-4}\end{array}\right).$$

(16)

We keep the eigenvalues $\{0.2325\times {10}^{-4}$, $0.4271\times {10}^{-4}\}$ of this historical covariance matrix $\widehat{C}$ and construct a $P_0$ with the same eigenvalues such that the corresponding correlation matrix $R_0$, given by $R_t={\Sigma }_t^{-1}{P}_t{\Sigma }_t^{-1}$ with ${\Sigma }_t={\left(diag\left(P_t\right)\right)}^{\frac{1}{2}}$, is a good approximation for the given historical correlation matrix at $t=0$. More precisely, we are looking for an orthogonal matrix $H\in {\mathbb{R}}^{3\times 3}$ such that

$$P_0=H^{\top }DH\mathrm{and}{\parallel R_0-R_0^{\mathrm{hist}}\parallel }_F\to min,$$

(17)

where D is a diagonal matrix containing the eigenvalues of $\widehat{C}$ multiplied with 1000. Since the eigenvalues of $\widehat{C}$ are very small, multiplying with the factor 1000 simplifies the optimization procedure of finding a compatible initial matrix $P_0$. The choice of this factor can be adapted according to the given historical data. In our experiments, the initial covariance matrix is found as

$$P_0=\left(\begin{array}{cc}0.0233& -0.0005\\ -0.0005& 0.0427\end{array}\right),$$

(18)

whereas the orthogonal matrix used to get $P_0$ is given by

$$H=\left(\begin{array}{cc}-0.9997& -0.0258\\ -0.0258& 0.9997\end{array}\right).$$

(19)

Due to the construction of our covariance flow it follows that every computed covariance matrix $P_t=Q_t^{\top }{P}_0{Q}_t$ will be positive semi-definite and contain information from the whole historical data.

4.1.4. Computation of Correlation Flows

Finally, we convert the covariance flows to correlation flows by

$$R_t={\Sigma }_t^{-1}{P}_t{\Sigma }_t^{-1},$$

(20)

where ${\Sigma }_t={\left(diag\left(P_t\right)\right)}^{\frac{1}{2}}$. The covariance flows or rather the correlation flows are computed such that the mean squared error between the empirical density function of the historical data $f^{\mathrm{hist}}\left(z\right)$ and the empirical density function of the correlation flow $f^{\mathrm{flow}}\left(z\right)$ is minimized

$$\frac{1}{N}\sum _{j=1}^N{\left(f^{\mathrm{hist}}\left(z_j\right)-f^{\mathrm{flow}}\left(z_j\right)\right)}^2\to min,$$

(21)

where N is the number of points where the kernel smoothing function estimate is evaluated at. For the computation of both empirical density functions, $f^{\mathrm{hist}}$ and $f^{\mathrm{flow}}$, we used the MATLAB function ksdensity.

4.2. Results

Taking up the risk manager’s problem we assume that he followed the steps described above and implemented them in a software package of his choice. Implementing the steps of our methodology in the software package MATLAB we found that the error defined in (21) is minimized by the following choice of parameters

$$\begin{array}{cc}\hfill x& =(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9)\hfill \\ \hfill & =(8.7492,-29.2169,21.4603,-14.7219,18.0489,-6.6036,-9.4284,-9.3455,0.0975),\hfill \end{array}$$

where $x_1,\cdots ,x_9$ are the parameters in (12) and (13). The mean squared error using these parameters is $3.527\times {10}^{-3}$. The corresponding density function of our correlation flow compared to the historical density function is plotted in Figure 3. The plot shows that the fitting by least-squares works well. The optimal bandwidth of the kernel smoothing function estimate of the historical data was computed to be 0.0559. We chose the same bandwidth for the density estimate of our correlation flow by inserting the bandwidth 0.0559 as an input argument in the computation of $f^{\mathrm{flow}}$ with the MATLAB function ksdensity. The density functions were evaluated at $N=100$ equally spaced points in the interval [−0.5473, 0.4788].

5. Discussion

We have presented a method to produce feasible covariance and correlation matrices. Based on isospectral flows we produced matrices similar to an initial valid covariance matrix, which we determined beforehand using historical data. In these covariance flows we assumed the required rotation matrices to be driven by an SDE in order to mimic the stochastic behaviour of correlations. These rotation matrices can be used to control the tendency of the corresponding correlation flows. For instance, one can require that the correlation flows match a desired density function.

Our generated covariance and correlation matrices are not only real symmetric and positive semi-definite but also exhibit stochastic behaviour. Hence, the correlation flows computed with our methodology have proven to fulfill the requirements for a valid correlation matrix while reflecting the randomness of correlations. In other words, we were able to expand the approach for constructing correlation flows presented in [11] by modelling an SDE for the underlying problem which we then solved by applying stochastic Lie group methods, i.e. we have shown that a financial mathematical problem can be solved by using methods of geometric numerical integration.

There are multiple possibilities to extend our presented methodology. For example, the function $f\left(t\right)$ used in the linear combination (11) can be chosen such that even more degrees of freedom are involved. In the case where $n>2$ correlations are considered and thus, a basis for $\mathfrak{so}\left(n\right)$ with more elements (similar to (9)) is needed, one could incorporate as many functions for the linear combination as there are basis matrices. For the numerical integration of (6) one could also use a method of higher order. The Euler-Maruyama scheme in Algorithm 1 could for example be replaced by the Milstein scheme. By analyzing modified equations one could construct even higher order methods, see e.g., [27].