Copula Processes

Introduction

Imagine measuring the distance of a rocket as it leaves Earth, and wanting to know how these measurements correlate with one another. How much does the value of the measurement at fifteen minutes depend on the measurement at five minutes? Once we’ve learned this correlation structure, suppose we want to compare it to the dependence between measurements of the rocket’s velocity. To do this, it is convenient to separate dependence from the marginal distributions of our measurements. At any given time, a rocket’s distance from Earth could have a Gamma distribution, while its velocity could have a Gaussian distribution. And separating dependence from marginal distributions is precisely what a copula function does.

While copulas have recently become popular, especially in financial applications , as Nelsen writes, “the study of copulas and the role they play in probability, statistics, and stochastic processes is a subject still in its infancy. There are many open problems…” Typically only bivariate (and recently trivariate) copulas are being used and studied. In our introductory example, we are interested in learning the correlations in different stochastic processes, and comparing them. It would therefore be useful to have a copula process, which can describe the dependencies between arbitrarily many random variables independently of their marginal distributions. We define such a process. As an example, we develop a stochastic volatility model, Gaussian Copula Process Volatility (GCPV). In doing so, we provide a Bayesian framework for the learning the marginal distributions and dependency structure of what we call a Gaussian copula process.

The volatility of a random variable is its standard deviation. Stochastic volatility models are used to predict the volatilities in a heteroscedastic sequence – a sequence of random variables with different variances, like distance measurements of a rocket as it leaves the Earth. As the rocket gets further away, the variance on the measurements increases. Heteroscedasticity is especially important in econometrics; the returns on equity indices, like the S&P 500, or on currency exchanges, are heteroscedastic. Indeed, in 2003, Robert Engle won the Nobel Prize in economics “for methods of analyzing economic time series with time-varying volatility”. GARCH , a generalized version of Engle’s ARCH, is arguably unsurpassed for predicting the volatility of returns on equity indices and currency exchanges . GCPV can outperform GARCH, and is competitive on financial data that especially suits GARCH . Before introducing GCPV, we first discuss copulas and then introduce our copula process. For a review of Gaussian processes, see Rasmussen and Williams .

Copulas

Copulas are important because they separate the dependency structure between random variables from their marginal distributions. Intuitively, we can describe the dependency structure of any multivariate joint distribution $H(x_{1},\dots,x_{n})=P(X_{1}\leq x_{1},\dots X_{n}\leq x_{n})$ through a two step process. First we take each univariate random variable $X_{i}$ and transform it through its cumulative distribution function (cdf) $F_{i}$ to get $U_{i}=F_{i}(X_{i})$, a uniform random variable. We then express the dependencies between these transformed variables through the $n$-copula $C(u_{1},\dots,u_{n})$. Formally, an $n$-copula $C:^{n}\to$ is a multivariate cdf with uniform univariate marginals: $C(u_{1},u_{2},\dots,u_{n})=P(U_{1}\leq u_{1},U_{2}\leq u_{2},\dots,U_{n}\leq u_{n})$, where $U_{1},U_{2},\dots,U_{n}$ are standard uniform random variables. Sklar precisely expressed our intuition in the theorem below.

Sklar’s theorem Let H be an n-dimensional distribution function with marginal distribution functions $F_{1},F_{2},\dots,F_{n}$. Then there exists an $n$-copula C such that for all $(x_{1},x_{2},\dots,x_{n})\in[-\infty,\infty]^{n}$,

If $F_{1},F_{2},\dots,F_{n}$ are all continuous then C is unique; otherwise C is uniquely determined on $\operatorname{Range}{F_{1}}\times\operatorname{Range}{F_{2}}\times\dots\times\operatorname{Range}{F_{n}}$. Conversely, if C is an $n$-copula and $F_{1},F_{2},\dots,F_{n}$ are distribution functions, then the function H is an n-dimensional distribution function with marginal distribution functions $F_{1},F_{2},\dots,F_{n}$.

As a corollary, if $F_{i}^{(-1)}(u)=\inf\{x:F(x)\geq u\}$, the quasi-inverse of $F_{i}$, then for all $u_{1},u_{2},\dots,u_{n}\in^{n}$,

In other words, (2) can be used to construct a copula. For example, the bivariate Gaussian copula is defined as

where $\Phi_{\rho}$ is a bivariate Gaussian cdf with correlation coefficient $\rho$, and $\Phi$ is the standard univariate Gaussian cdf. Li popularised the bivariate Gaussian copula, by showing how it could be used to study financial risk and default correlation, using credit derivatives as an example.

By substituting $F(x)$ for $u$ and $G(y)$ for $v$ in equation (3), we have a bivariate distribution $H(x,y)$, with a Gaussian dependency structure, and marginals $F$ and $G$. Regardless of $F$ and $G$, the resulting $H(x,y)$ can still be uniquely expressed as a Gaussian copula, so long as $F$ and $G$ are continuous. It is then a copula itself that captures the underlying dependencies between random variables, regardless of their marginal distributions. For this reason, copulas have been called dependence functions . Nelsen contains an extensive discussion of copulas.

Copula Processes

Imagine choosing a covariance function, and then drawing a sample function at some finite number of points from a Gaussian process. The result is a sample from a collection of Gaussian random variables, with a dependency structure encoded by the specified covariance function. Now, suppose we transform each of these values through a univariate Gaussian cdf, such that we have a sample from a collection of uniform random variables. These uniform random variables also have this underlying Gaussian process dependency structure. One might call the resulting values a draw from a Gaussian-Uniform process. We could subsequently put these values through an inverse beta cdf, to obtain a draw from what could be called a Gaussian-Beta process: the values would be a sample from beta random variables, again with an underlying Gaussian process dependency structure. Alternatively, we could transform the uniform values with different inverse cdfs, which would give a sample from different random variables, with dependencies encoded by the Gaussian process.

The above procedure is a means to generate samples from arbitrarily many random variables, with arbitrary marginal distributions, and desired dependencies. It is an example of how to use what we call a copula process – in this case, a Gaussian copula process, since a Gaussian copula describes the underlying dependency structure of a finite number of samples. We can now formally define a copula process.

Note that each $Q_{t_{i}}\sim\operatorname{Uniform}(0,1)$, and that $G_{t_{i}}^{(-1)}$ is the quasi-inverse of $G_{t_{i}}$, as it was previously defined.

Gaussian Copula Process $W_{t}$ is Gaussian copula process distributed if it is copula process distributed and the base measure $\mu$ is a Gaussian process. If there is a mapping $\Psi$ such that $\Psi(W_{t})\sim\mathcal{GP}(m(t),k(t,t^{\prime}))$, then we write $W_{t}\sim\mathcal{GCP}(\Psi,m(t),k(t,t^{\prime}))$.

For example, if we have $W_{t}\sim\mathcal{GCP}$ with $m(t)=0$ and $k(t,t)=1$, then in our definition of a copula process, $G_{t}=\Phi$, the standard univariate Gaussian cdf, and $H$ is the usual GP joint distribution function. Supposing this GCP is a Gaussian-Beta process, then $\Psi=\Phi^{-1}\circ F_{B}$, where $F_{B}$ is a univariate Beta cdf. One could similarly define other copula processes.

We described generally how a copula process can be used to generate samples of arbitrarily many random variables with desired marginals and dependencies. We now develop a specific and practical application of this framework. We introduce a stochastic volatility model, Gaussian Copula Process Volatility (GCPV), as an example of how to learn the joint distribution of arbitrarily many random variables, the marginals of these random variables, and to make predictions. To do this, we fit a Gaussian copula process by using a type of Warped Gaussian Process . However, our methodology varies substantially from Snelson et al. , since we are doing inference on latent variables as opposed to observations, which is a much greater undertaking that involves approximations, and we are doing so in a different context.

Gaussian Copula Process Volatility

Assume we have a sequence of observations $\bm{y}=(y_{1},\dots,y_{n})^{\top}$ at times $\bm{t}=(t_{1},\dots,t_{n})^{\top}$. The observations are random variables with different latent standard deviations. We therefore have $n$ unobserved standard deviations, $\sigma_{1},\dots,\sigma_{n}$, and want to learn the correlation structure between these standard deviations, and also to predict the distribution of $\sigma_{*}$ at some unrealised time $t_{*}$.

We model the standard deviation function as a Gaussian copula process:

where $g$ is a monotonic warping function, parametrized by $\bm{\omega}$. For each of the observations $\bm{y}=(y_{1},\dots,y_{n})^{\top}$ we have corresponding GP latent function values $\bm{f}=(f_{1},\dots,f_{n})^{\top}$, where $\sigma(t_{i})=g(f_{i},\bm{\omega})$, using the shorthand $f_{i}$ to mean $f(t_{i})$.

$\sigma_{t}\sim\mathcal{GCP}$, because any finite sequence $(\sigma_{1},\dots,\sigma_{p})$ is distributed as a Gaussian copula:

Here we have assumed that each observation, conditioned on knowing its variance, is normally distributed with zero mean. This is a common assumption in heteroscedastic models. The zero mean and normality assumptions can be relaxed and are not central to this paper.

Predictions with GCPV

Ultimately, we wish to infer $p(\sigma(t_{*})|\bm{y},\bm{z})$, where $\bm{z}=\{\bm{\theta},\bm{\omega}\}$, and $\bm{\theta}$ are the hyperparameters of the GP covariance function. To do this, we sample from

and then transform these samples by $g$. Letting $(C_{\bm{f}})_{ij}$ = $\delta_{ij}g(f_{i},\bm{\omega})^{2}$, where $\delta_{ij}$ is the Kronecker delta, $K_{ij}=k(t_{i},t_{j})$, $(\bm{k}_{*})_{i}=k(t_{*},t_{i})$, we have

We also wish to learn $\bm{z}$, which we can do by finding the $\bm{\hat{z}}$ that maximizes the marginal likelihood,

Unfortunately, for many functions $g$, (10) and (13) are intractable. Our methods of dealing with this can be used in very general circumstances, where one has a Gaussian process prior, but an (optionally parametrized) non-Gaussian likelihood. We use the Laplace approximation to estimate $p(\bm{f}|\bm{y},\bm{z})$ as a Gaussian. Then we can integrate (10) for a Gaussian approximation to $p(f_{*}|\bm{y},\bm{z})$, which we sample from to make predictions of $\sigma_{*}$. Using Laplace, we can also find an expression for an approximate marginal likelihood, which we maximize to determine $\bm{z}$. While we always use Laplace to determine $\bm{z}$, we compare to a full Laplace solution by also using Markov chain Monte Carlo to sample from $p(f_{*}|\bm{y},\bm{z})$.

Let us now relate the above to the Gaussian copula in (9). The prior $\Gamma_{ij}=\textrm{cov}(g^{-1}(\sigma_{i}),g^{-1}(\sigma_{j}))=\textrm{cov}(f_{i},f_{j})=k(t_{i},t_{j})$. The posterior $\Gamma_{ij}$ can be estimated as the covariance matrix of the Laplace approximation for $p(\bm{f}|\bm{y})$. Also, since each component of $\bm{f}$ is transformed separately, such that $\sigma(t_{i})=g(f(t_{i}))$, we have

One can use this to simulate from the joint distribution over the deviations.

The goal is to approximate (11) with a Gaussian, so that we can evaluate (10) and (13) and make predictions. In doing so, we follow Rasmussen and Williams in their treatment of Gaussian process classification, except we use a parametrized likelihood, and modify Newton’s method.

First, consider as an objective function the logarithm of an unnormalized (11):

where $W$ is the diagonal matrix $-\nabla\nabla\log p(\bm{y}|\bm{f},\bm{\omega})$.

If the likelihood function $p(\bm{y}|\bm{f},\bm{\omega})$ is not log concave, then $W$ may have negative entries. Vanhatalo et al. found this to be problematic when doing Gaussian process regression with a Student-t likelihood. They instead use an expectation-maximization (EM) algorithm for finding $\bm{\hat{f}}$, and iterate ordered rank one Cholesky updates to evaluate the Laplace approximate marginal likelihood. But EM can converge slowly, especially near a local optimum, and each of the rank one updates is vulnerable to numerical instability. With a small modification of Newton’s method, we often get close to quadratic convergence for finding $\bm{\hat{f}}$, and can evaluate the Laplace approximate marginal likelihood in a numerically stable fashion, with no approximate Cholesky factors, and optimal computational requirements.

At a maximum, the negative Hessian of the objective function, $W+K^{-1}$, is positive definite. On each iteration of Newton’s method, we form $M$ by setting all negative entries of $W$ to zero. Since $K^{-1}$ is positive definite, and the eigenvalues of $M+K^{-1}$ are greater than or equal to the eigenvalues of $K^{-1}$, $M+K^{-1}$ is always positive definite. Using $M$ in place of $W$ decreases the Newton step size, and changes the direction of steps. We are always stepping towards a local maximum, and will converge, barring rare pathologies.

where $B_{\bm{\hat{f}}}$ is $B$ evaluated at $\bm{\hat{f}}$, and $\bm{a}_{\bm{\hat{f}}}$ is a numerically stable evaluation of $K^{-1}\bm{\hat{f}}$.

Given $n$ training observations, the cost of each Newton iteration is dominated by computing $chol(B)$, which takes $\mathcal{O}(n^{3})$ operations. The objective function typically changes by less than $10^{-6}$ after 3 iterations. Once Newton’s method has converged, it takes only $\mathcal{O}(1)$ operations to draw from $q(f_{*}|\bm{y},\bm{z})$ and make predictions.

2 Markov chain Monte Carlo

We use Markov chain Monte Carlo (MCMC) to sample from (11), so that we can later sample from $p(\sigma_{*}|\bm{y},\bm{z})$ to make predictions. Sampling from (11) is difficult, because the variables $\bm{f}$ are strongly coupled by a Gaussian process prior. We use a new technique, Elliptical Slice Sampling , and find it extremely effective for this purpose. It was specifically designed to sample from posteriors with correlated Gaussian priors. It has no free parameters, and jointly updates every element of $\bm{f}$. For our setting, it is over 100 times as fast as axis aligned slice sampling with univariate updates.

To make predictions, we take $J$ samples of $p(\bm{f}|\bm{y},\bm{z})$, {$\bm{f}^{1},\dots,\bm{f}^{J}$}, and then approximate (10) as a mixture of $J$ Gaussians:

Each of the Gaussians in this mixture have equal weight. So for each sample of $f_{*}|\bm{y}$, we uniformly choose a random $p(f_{*}|\bm{f}^{i},\bm{\theta})$ and draw a sample. In the limit $J\to\infty$, we are sampling from the exact $p(f_{*}|\bm{y},\bm{z})$. Mapping these samples through $g$ gives samples from $p(\sigma_{*}|\bm{y},\bm{z})$.

After one $\mathcal{O}(n^{3})$ and one $\mathcal{O}(J)$ operation, a draw from (20) takes $\mathcal{O}(1)$ operations.

3 Warping Function

This is monotonic, positive, infinitely differentiable, asymptotic towards zero as $x\to-\infty$, and tends to $(\sum_{j=1}^{K}a_{j}b_{j})x$ as $x\to\infty$. In practice, it is useful to add a small constant to (21), to avoid rare situations where the parameters $\bm{\omega}$ are trained to make $g$ extremely small for certain inputs, at the expense of a good overall fit; this can happen when the parameters $\bm{\omega}$ are learned by optimizing a likelihood. A suitable constant could be one tenth the absolute value of the smallest nonzero observation.

By inferring the parameters of the warping function, or distributions of these parameters, we are learning a transformation which will best model $\sigma_{t}$ with a Gaussian process. The more flexible the warping function, the more potential there is to improve the GCPV fit – in other words, the better we can estimate the ‘perfect’ transformation. To test the importance of this flexibility, we also try a simple unparametrized warping function, $g(x)=e^{x}$. In related work, Goldberg et al. place a GP prior on the log noise level in a standard GP regression model on observations, except for inference they use Gibbs sampling, and a high level of ‘jitter’ for conditioning.

Experiments

We then sample from $p(f_{*}|\bm{y})$ using the Laplace approximation (19). We also do this using MCMC (20) with $J=10000$, after discarding a previous $10000$ samples of $p(\bm{f}|\bm{y})$ as burn-in. We pass these samples of $f_{*}|\bm{y}$ through $g$ and $g^{2}$ to draw from $p(\sigma_{*}|\bm{y})$ and $p(\sigma_{*}^{2}|\bm{y})$, and compute the sample mean and variance of $\sigma_{*}|\bm{y}$. We use the sample mean as a point predictor, and the sample variance for error bounds on these predictions, and we use $10000$ samples to compute these quantities. For GCPV we use Laplace and MCMC for inference, but for GP-EXP we only use Laplace. We compare predictions to GARCH(1,1), which has been shown in extensive and recent reviews to be competitive with other GARCH variants, and more sophisticated models . We use the Matlab Econometrics Toolbox implementation of GARCH.

We make forecasts of volatility, and we predict historical volatility. By ‘historical volatility’ we mean the volatility at observed time points, or between these points. Uncovering historical volatility is important. It could, for instance, be used to study what causes fluctuations in the stock market, or to understand physical systems.

To evaluate our model, we use the Mean Squared Error (MSE) between the true variance, or proxy for the truth, and the predicted variance. Although likelihood has advantages, we are limited in space, and we wish to harmonize with the econometrics literature, and other assessments of volatility models, where MSE is the standard. In a similar assessment of volatility models, Brownlees et al. found that MSE and quasi-likelihood rankings were comparable.

We simulate observations from $\mathcal{N}(0,\sigma^{2}(t))$, using $\sigma(t)=\sin(t)\cos(t^{2})+1$, at $\bm{t}=(0,0.02,0.04,\dots,4)^{\top}$. We call this data set TRIG. We also simulate using a standard deviation that jumps from $0.1$ to $7$ and back, at times $\bm{t}=(0,0.1,0.2,\dots,6)^{\top}$. We call this data set JUMP. To forecast, we use all observations up until the current time point, and make 1, 7, and 30 step ahead predictions. So, for example, in TRIG we start by observing $t=0$, and make forecasts at $t=0.02,0.14,0.60$. Then we observe $t=0,0.02$ and make forecasts at $t=0.04,0.16,0.62$, and so on, until all data points have been observed. For historical volatility, we predict the latent $\sigma_{t}$ at the observation times, which is safe since we are comparing to the true volatility, which is not used in training; the results are similar if we interpolate. Figure 1 panels a) and b) show the true volatility for TRIG and JUMP respectively, alongside GCPV Laplace, GCPV MCMC, GP-EXP Laplace, and GARCH(1,1) predictions of historical volatility. Table 1 shows the results for forecasting and historical volatility.

In panel a) we see that GCPV more accurately captures the dependencies between $\sigma$ at different times points than GARCH: if we manually decrease the lengthscale in the GCPV covariance function, we can replicate the erratic GARCH behaviour, which inaccurately suggests that the covariance between $\sigma_{t}$ and $\sigma_{t+a}$ decreases quickly with increases in $a$. We also see that GCPV with an unparametrized exponential warping function tends to overestimates peaks and underestimate troughs. In panel b), the volatility is extremely difficult to reconstruct or forecast – with no warning it will immediately and dramatically increase or decrease. This behaviour is not suited to a smooth squared exponential covariance function. Nevertheless, GCPV outperforms GARCH, especially in regions of low volatility. We also see this in panel a) for $t\in(1.5,2)$. GARCH is known to respond slowly to large returns, and to overpredict volatility . In JUMP, the greater the peaks, and the smaller the troughs, the more GARCH suffers, while GCPV is mostly robust to these changes.

2 Financial Data

The returns on the daily exchange rate between the Deutschmark (DM) and the Great Britain Pound (GBP) from 1984 to 1992 have become a benchmark for assessing the performance of GARCH models . This exchange data, which we refer to as DMGBP, can be obtained from www.datastream.com, and the returns are calculated as $r_{t}=\log(P_{t+1}/P_{t})$, where $P_{t}$ is the number of DM to GBP on day $t$. The returns are assumed to have a zero mean function.

We use a rolling window of the previous 120 days of returns to make 1, 7, and 30 day ahead volatility forecasts, starting at the beginning of January 1988, and ending at the beginning of January 1992 (659 trading days). Every 7 days, we retrain the parameters of GCPV and GARCH. Every time we retrain parameters, we predict historical volatility over the past 120 days. The average MSE for these historical predictions is given in Table 1, although they should be observed with caution; unlike with the simulations, the DMGBP historical predictions are trained using the same data they are assessed on. In Figure 1c), we see that the GARCH one day ahead forecasts are lifted above the GCPV forecasts, but unlike in the simulations, they are now operating on a similar lengthscale. This suggests that GARCH could still be overpredicting volatility, but that GCPV has adapted its estimation of how $\sigma_{t}$ and $\sigma_{t+a}$ correlate with one another. Since GARCH is suited to this financial data set, it is reassuring that GCPV predictions have a similar time varying structure. Overall, GCPV and GARCH are competitive with one another for forecasting currency exchange returns, as seen in Table 1. Moreover, a learned warping function $g$ outperforms an unparametrized one, and a full Laplace solution is comparable to using MCMC for inference, in accuracy and speed. This is also true for the simulations. Therefore we recommend whichever is more convenient to implement.

Discussion

We defined a copula process, and as an example, developed a stochastic volatility model, GCPV, which can outperform GARCH. With GCPV, the volatility $\sigma_{t}$ is distributed as a Gaussian Copula Process, which separates the modelling of the dependencies between volatilities at different times from their marginal distributions – arguably the most useful property of a copula. Further, GCPV fits the marginals in the Gaussian copula process by learning a warping function. If we had simply chosen an unparametrized exponential warping function, we would incorrectly be assuming that the log volatilities are marginally Gaussian distributed. Indeed, for the DMGBP data, we trained the warping function $g$ over a 120 day period, and mapped its inverse through the univariate standard Gaussian cdf $\Phi$, and differenced, to estimate the marginal probability density function (pdf) of $\sigma_{t}$ over this period. The learned marginal pdf, shown in Figure 1d), is similar to a Gamma(4.15,0.00045) distribution. However, in using a rolling window to retrain the parameters of $g$, we do not assume that the marginals of $\sigma_{t}$ are stationary; we have a time changing warping function.

While GARCH is successful, and its simplicity is attractive, our model is also simple and has a number of advantages. We can effortlessly handle missing data, we can easily incorporate covariates other than time (like interest rates) in our covariance function, and we can choose from a rich class of covariance functions – squared exponential, Brownian motion, Matérn, periodic, etc. In fact, the volatility of high frequency intradaily returns on equity indices and currency exchanges is cyclical , and GCPV with a periodic covariance function is uniquely well suited to this data. And the parameters of GCPV, like the covariance function lengthscale, or the learned warping function, provide insight into the underlying source of volatility, unlike the parameters of GARCH.

Finally, copulas are rapidly becoming popular in applications, but often only bivariate copulas are being used. We introduced a copula process where one can learn the dependency structure between arbitrarily many random variables independently of their marginal distributions. We hope the Gaussian Copula Process Volatility model will encourage other applications of copula processes. More generally, we hope our work will help bring together the machine learning and econometrics communities.

Acknowledgments: Thanks to Carl Edward Rasmussen and Ferenc Huszár for helpful conversations. AGW is supported by an NSERC grant.

References