A Short Note on Gaussian Distribution with Non-Constant Correlation

1. Background

Gaussian copula is widely used in quantitative finance modelling. The Gaussian distribution is closely related to an underlying Brownian motion: the standard multi-variate normal distribution is the terminal distribution of an underlying multi-variate Brownian motion where the correlations are constant over time. However the correlation being constant is a limitation of this model which might not fit the actual market. On the other hand, if the correlations are not constant, the result terminal distribution has no closed-form representation in general. Without the closed-form solution or analytic tractability, it becomes less attractive for practical usage. There are research in alternative directions which bypass this tractability issue, for example in [1,2], the respective authors created different terminal distributions which can admit shape with the desired correlation skew effect. In this paper, we still focus on the terminal distribution result from the Brownian motion itself. We study the PDE for the density function and show that with some assumption on the correlation function, the PDE can be decomposed to lower dimensional ones and therefore make the calculation fast and practical. With this technique, the result distribution can be a useful variation to the standard multi-variate normal distribution and it can be used for purpose like modelling correlation skew effect in quant finance.

We also think the terminal distribution with non-constant correlation might be an interesting mathmatical object on itself.

2. Methodology

We study this math problem below. This is a 2-dimensional case however we show later that similar techniques can be applied to higher dimensions.

$$\begin{array}{ccc}\hfill x\left(0\right),y\left(0\right)& =& 0,0\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill dx& =& d{w}_1\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill dy& =& \rho (x,y,t)d{w}_1+\sqrt{1-{\rho }^2(x,y,t)}d{w}_2\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill <d{w}_1,d{w}_2>& =& 0\hfill \end{array}$$

(4)

The Fokker-Planck equation [3,4,5] describes the joint probability density function $p(x,y,t)$ by:

$${\displaystyle \frac{\partial p}{\partial t}}={\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial p}{\partial x^2}}+2{\displaystyle \frac{\partial \left(\rho p\right)}{\partial x\partial y}}+{\displaystyle \frac{\partial p}{\partial y^2}})$$

(5)

This is a 2d-PDE in the convention of quant finance industry (2d refers to 2-dimension in space variables $(x,y)$ while in fact it is a 3-d PDE if counting t, given the common presence of t in this type of PDE we refer the dimensions to only the space variables) and the general numerical method is slow. However, we can decompose the 2d-PDE into two 1d-PDEs if we make a reasonable assumtion on the correlation function as below:

$$\rho (x,y,t)=\rho (x+y,t)$$

(6)

This means the correlation depends on the $(x,y)$ in terms of the total $(x+y)$, which can be interpreted as: correlation depends on a market factor which is the average of the underlyers. With this extra assumption, we can simplify the problem as below:

Lets make change of variables below

$$\begin{array}{ccc}\hfill u& =& {\displaystyle \frac{1}{2}}(x+y)\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill v& =& {\displaystyle \frac{1}{2}}(x-y)\hfill \end{array}$$

(8)

Then we have

$$<du,dv>={\displaystyle \frac{1}{4}}(<dx,dx>-<dy,dy>)=0$$

(9)

And $du,dv$ can be written as

$$\begin{array}{ccc}\hfill du& =& \sqrt{{\displaystyle \frac{1+\rho (u,t)}{2}}}d{w}_3\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill dv& =& \sqrt{{\displaystyle \frac{1-\rho (u,t)}{2}}}d{w}_4\hfill \end{array}$$

(11)

Note the first equation only involves u, then Fokker-Planck equation for u is a 1d-PDE:

$${\displaystyle \frac{\partial p(u,t)}{\partial t}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2}{\partial u^2}}\left({\displaystyle \frac{1+\rho (u,t)}{2}}p(u,t)\right)$$

(12)

So we can solve $p(u,t)$ first, then we look at the $v\left(t\right)$. For any given path $u\left(s\right),0\le s\le t$, the $v\left(t\right)$ is simply a sum of infinitesimal normal variables with variances ${\displaystyle \frac{1-\rho \left(u\right(s),s)}{2}}$, so we know the distrubtion of $v\left(t\right)$ condition on this path $u\left(s\right),0\le s\le t$ is a normal distribution with mean 0 and variance

$${\int }_0^t{\displaystyle \frac{1-\rho \left(u\right(s),s)}{2}}ds$$

(13)

Conditioned on a path is not easy to use for calculation, it would be more useful to condition on a value $u\left(t\right)$ instead of the whole path. Lets consider the conditional expectation

$$f(u,t)=E\left[{\int }_0^t{\displaystyle \frac{1-\rho \left(u\right(s),s)}{2}}ds\right|u\left(t\right)=u]$$

(14)

This is the path integral on all possible paths $u\left(s\right)$ that get to u at t. We have the following:

$$p(u,t+dt)f(u,t+dt)={\int }_{-\infty }^{\infty }p(x,t)[f(x,t)+dt{\displaystyle \frac{1-\rho (x,t)}{2}}]p(u,t+dt|x,t)dx$$

(15)

The $p(u,t+dt|x,t)$ is the transition probability from state $(x,t)$ to $(u,t+dt)$.

Now we follow the Fokker-Planck equation derivation technique, we will get:

$${\displaystyle \frac{\partial }{\partial t}}\left(pf\right)=p{\displaystyle \frac{1-\rho (u,t)}{2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial }{\partial u^2}}\left(pf{\displaystyle \frac{1+\rho (u,t)}{2}}\right)$$

(16)

The proof is standard derivation, readers can skip it. For completeness we include the outline below:

Outline of proof:

$$\begin{array}{ccc}& & \underset{dt\to 0}{lim}{\displaystyle \frac{p(u,t+dt)f(u,t+dt)-p(u,t)f(u,t)}{dt}}\hfill \\ & =& \underset{dt\to 0}{lim}{\displaystyle \frac{{\int }_{-\infty }^{\infty }p(x,t)f(x,t)p(u,t+dt|x,t)dx-p(u,t)f(u,t)}{dt}}+\underset{dt\to 0}{lim}{\int }_{-\infty }^{\infty }p(x,t){\displaystyle \frac{1-\rho (x,t)}{2}}p(u,t+dt|x,t)dx\hfill \end{array}$$

Note the second term comes to

$$\underset{dt\to 0}{lim}{\int }_{-\infty }^{\infty }p(x,t){\displaystyle \frac{1-\rho (x,t)}{2}}p(u,t+dt|x,t)dx=p(u,t){\displaystyle \frac{1-\rho (u,t)}{2}}$$

So we just have to prove

$$\underset{dt\to 0}{lim}{\displaystyle \frac{{\int }_{-\infty }^{\infty }p(x,t)f(x,t)p(u,t+dt|x,t)dx-p(u,t)f(u,t)}{dt}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial }{\partial u^2}}\left(pf{\displaystyle \frac{1+\rho (u,t)}{2}}\right)$$

Let $h\left(u\right)$ be a smooth function with compact support, consider

$$\begin{array}{ccc}& & {\int }_{-\infty }^{\infty }h\left(u\right){\int }_{-\infty }^{\infty }p(x,t)f(x,t)p(u,t+dt|x,t)dxdu\hfill \\ & =& {\int }_{-\infty }^{\infty }p(x,t)f(x,t){\int }_{-\infty }^{\infty }h\left(u\right)p(u,t+dt|x,t)dudx\hfill \\ & =& {\int }_{-\infty }^{\infty }p(x,t)f(x,t){\int }_{-\infty }^{\infty }(h\left(x\right)+h^{\prime }(u-x)+{\displaystyle \frac{1}{2}}{h}^{\prime \prime }{(u-x)}^2+O\left({(u-x)}^3\right))p(u,t+dt|x,t)dudx\hfill \end{array}$$

Now the integral ${\int }_{-\infty }^{\infty }{(u-x)}^{k}p(u,t+dt|x,t)du$ is the $k-th$ moment of the Brownian motion $du=\sqrt{{\displaystyle \frac{1+\rho (u,t)}{2}}}d{w}_3$, so we have

$$\begin{array}{ccc}\hfill {\int }_{-\infty }^{\infty }(u-x)p(u,t+dt|x,t)du& =& 0\hfill \\ \hfill {\int }_{-\infty }^{\infty }{(u-x)}^{2}p(u,t+dt|x,t)du& =& {\displaystyle \frac{1+\rho (x,t)}{2}}dt\hfill \\ \hfill {\int }_{-\infty }^{\infty }{(u-x)}^{k}p(u,t+dt|x,t)du& =& \mathrm{higher}\mathrm{order}\mathrm{than}dt\mathrm{when}k>2\hfill \end{array}$$

Then we have below, in the order of $dt$

$$\begin{array}{ccc}& & {\int }_{-\infty }^{\infty }h\left(u\right){\int }_{-\infty }^{\infty }p(x,t)f(x,t)p(u,t+dt|x,t)dxdu\hfill \\ & =& {\int }_{-\infty }^{\infty }p(x,t)f(x,t)h\left(x\right){\int }_{-\infty }^{\infty }p(u,t+dt|x,t)dudx+dt{\int }_{-\infty }^{\infty }{\displaystyle \frac{1}{2}}p(x,t)f(x,t)h^{\prime \prime }{\displaystyle \frac{1+\rho (x,t)}{2}}dx\hfill \\ & =& {\int }_{-\infty }^{\infty }p(x,t)f(x,t)h\left(x\right)dx+dt{\int }_{-\infty }^{\infty }{\displaystyle \frac{1}{2}}p(x,t)f(x,t)h^{\prime \prime }{\displaystyle \frac{1+\rho (x,t)}{2}}dx\hfill \end{array}$$

so

$$\begin{array}{ccc}& & \underset{dt\to 0}{lim}{\displaystyle \frac{{\int }_{-\infty }^{\infty }h\left(u\right){\int }_{-\infty }^{\infty }p(x,t)f(x,t)p(u,t+dt|x,t)dxdu-{\int }_{-\infty }^{\infty }h\left(u\right)p(u,t)f(u,t)du}{dt}}\hfill \\ & =& {\int }_{-\infty }^{\infty }{\displaystyle \frac{1}{2}}p(x,t)f(x,t)h^{\prime \prime }{\displaystyle \frac{1+\rho (x,t)}{2}}dx\hfill \\ & =& {\int }_{-\infty }^{\infty }h\left(x\right){\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left(p(x,t)f(x,t){\displaystyle \frac{1+\rho (x,t)}{2}}\right)dx\hfill \end{array}$$

The last step in above is integration by parts. Because the $h\left(u\right)$ is arbitrary smooth function so it follows that:

$$\underset{dt\to 0}{lim}{\displaystyle \frac{{\int }_{-\infty }^{\infty }p(x,t)f(x,t)p(u,t+dt|x,t)dx-p(u,t)f(u,t)}{dt}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial }{\partial u^2}}\left(pf{\displaystyle \frac{1+\rho (u,t)}{2}}\right)$$

End of Proof.

To recap, we have these 2 key equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial p(u,t)}{\partial t}}& =& {\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^2}{\partial u^2}}\left({\displaystyle \frac{1+\rho (u,t)}{2}}p(u,t)\right)\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial }{\partial t}}\left(pf\right)& =& p{\displaystyle \frac{1-\rho \left(u\right(t),t)}{2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial }{\partial u^2}}\left(pf{\displaystyle \frac{1+\rho (u,t)}{2}}\right)\hfill \end{array}$$

(18)

We can solve for p first and then solve for f (It is also possible to bundle the PDE solving for p and f together in discretization etc). Knowing $p(u,t)$ and $f(u,t)$, the whole distribution is known. The key point here is when solving p or f it is a low dimension PDE.

3. Higher Dimensions

In higher dimensions, similar technique can be applied if we assume the correlations have a dependency on one variable (though the variable might be defined as a linear combination of the base variables) and time only. A brief walk through of the idea as below:

Let $x_1,x_2,...,x_n$ be the initial Brownian motion variables with correlations ${\rho }_{ij}(x_1,x_2,...,x_n,t)$. For simplicity and avoid any singularity questions, lets assume the ${\rho }_{ij}$ all just depend on variable $M={\displaystyle \frac{1}{n}}{\sum }_i{x}_i$ and t. Now we can represent the random process by new set variables $M,x_2,...,x_n$, ($x_1$ is left out as it can be implied by others). We can do Cholesky decomposition of this set of variables and a nice property is that the Cholesky matrix elements are all just function of M and t: this is because all the ${\rho }_{ij}$ are just function of M and t and the Cholesky decomposition is a deterministic operation on those ${\rho }_{ij}$. Now we can apply the same process as before, solve for the probability density of M, and then for variance function of each of the independent Brownian motion variables coming from the Cholesky decomposition. The process will be long and tedious but there is no conceptional difference to the previous case. Note in a special case where all the ${\rho }_{ij}(M,t)$ are the same function, the change of variables is much cleaner and easier.

So in general, with the assumption that the correlation dependency degenerated to one variable only, the joint terminal distribution of dimension n can be calculated by n 1-d PDEs (1-d refer to 1 space variable).

4. Implementaion Example

We show one example of 2-d case: We discretize p and f together and solve for p first for a time step, and then solve for f. We don’t use chain rule to break out the partial derivatives of product but instead discretize on the product. With standard finite difference methods, the calculation is fast and stable. We present an example of the distribution below:

Figure 1 shows the contour of a Gaussian distribution with correlation skew. The underlying correlation function is:

$$\rho \left(u\right)=\left\{\begin{array}{cc}0.9\hfill & \mathrm{if}u<-2\sqrt{t}\hfill \\ 0.9-{\displaystyle \frac{u+2\sqrt{t}}{4\sqrt{t}}}0.4\hfill & \mathrm{if}-2\sqrt{t}\le u<2\sqrt{t}\hfill \\ 0.5\hfill & \mathrm{if}u\ge 2\sqrt{t}\hfill \end{array}\right.$$

The graph axis is in x and y. Note $u,v$ will be the two diagonal directions.

The shape of the contour is expected. As we put higher correltion when the $u={\displaystyle \frac{x+y}{2}}$ is lower, and lower correlation when u is higher, the probability is more concentrated when u is low and more dispersed when u is high. Note with u fixed, the graph also shows symmetry in the direction of v.

The following graphs shows more details on $p\left(u\right)$ in above example.

In Figure 2 the distribution of u is very close but different to a standard normal. To see the difference, we reflected the probability around center and then one can see the negative part has a fatter tail than positive part. This is expected as we correlated $x,y$ more when $x+y$ is more negative, we expect $x+y$ will have more potential to go lower in the negative direction, and as we de-correlate $x,y$ more when $x+y$ more positive, we expect the diversifying effect makes the $x+y$ less potential to go higher when $x+y$ positive.

To demonstrate this point, we can increase the skew of correlation further to see the fat tail effect. Below is the $p\left(u\right)$ for a more skewed correlation function.

The correlation function in Figure 3 is

$$\rho \left(u\right)=\left\{\begin{array}{cc}0.9\hfill & \mathrm{if}u<-2\sqrt{t}\hfill \\ 0.9-{\displaystyle \frac{u+2\sqrt{t}}{4\sqrt{t}}}1.4\hfill & \mathrm{if}-2\sqrt{t}\le u<2\sqrt{t}\hfill \\ -0.5\hfill & \mathrm{if}u\ge 2\sqrt{t}\hfill \end{array}\right.$$

The coutour in Figure 1 shows the v is concentrated when u more negative and v is spreaded when u is more positive.

Below Figure 4 shows the std dev of v conditioned on u, ie, the $\sqrt{f}$ function.

5. Copula Application

Knowing the $p\left(u\right)$ and $f\left(u\right)$ we can integrate any function on this terminal distribution. For given $u,v$, it maps to $x,y$, and the marginal distribution of x and y are still normal distribution respectively, so they can be readily used to invert CDFs.