Generalised Wishart Processes

Introduction

Imagine the price of the NASDAQ composite index increased dramatically today. Will it continue to rise tomorrow? Should you invest? Perhaps this is the beginning of a trend, but it may also be an anomaly. Now suppose you discover that every major equity index – FTSE, NIKKEI, TSE, etc. – has also risen. Instinctively the rise in NASDAQ was not anomalous, because this market is correlated with other major indices. This is an example of how multivariate models which account for correlations can be better than univariate models at making univariate predictions.

In this paper, we are concerned with modelling the dynamic covariance matrix $\Sigma(t)$ of high dimensional data sets (multivariate volatility). These models are especially important in econometrics. Brownlees et al., (2009) remark that “The price of essentially every derivative security is affected by swings in volatility. Risk management models used by financial institutions and required by regulators take time-varying volatility as a key input. Poor appraisal of the risks to come can leave investors excessively exposed to market fluctuations or institutions hanging on a precipice of inadequate capital”. Indeed, Robert Engle and Clive Granger won the 2003 Nobel prize in economics “for methods of analysing economic time series with time-varying volatility”. The returns on major equity indices and currency exchanges are thought to have a time changing variance and zero mean, and GARCH (Bollerslev,, 1986), a generalization of Engle’s ARCH (Engle,, 1982), is arguably unsurpassed at predicting the volatilities of returns on these equity indices and currency exchanges (Poon and Granger,, 2005; Hansen and Lunde,, 2005; Brownlees et al.,, 2009). Multivariate volatility models can be used to understand the dynamic correlations (or co-movement) between equity indices, and can make better univariate predictions than univariate models. A good estimate of the covariance matrix $\Sigma(t)$ is also necessary for portfolio management. An optimal portfolio allocation $\bm{w}^{*}$ is said to maximise the Sharpe ratio (Sharpe,, 1966):

where $\bm{r}(t)$ are expected returns for each asset and $\Sigma(t)$ is the predicted covariance matrix for these returns. One may also wish to maximise the portfolio return $\bm{w}^{\top}\bm{r}(t)$ for a fixed level of volatility: $\sqrt{\bm{w}^{\top}\Sigma(t)\bm{w}}=\lambda$. Sharpe, Markowitz and Merton jointly received a Nobel prize for portfolio theory (Markowitz,, 1952; Merton,, 1972). Multivariate volatility models are also used to understand contagion: the transmission of a financial shock from one entity to another (Bae et al.,, 2003). And generally – in econometrics, machine learning, climate science, or otherwise – it is useful to know the dynamic correlations between multiple entities.

Despite their importance, existing multivariate volatility models suffer from tractability issues and a lack of generality. For example, multivariate GARCH (MGARCH) has a number of free parameters that scales with dimension to the fourth power, and interpretation and estimation of these parameters is difficult to impossible (Silvennoinen and Teräsvirta,, 2009; Gouriéroux,, 1997), given the constraint that $\Sigma(t)$ must be positive definite at all points in time. Thus MGARCH, and alternative multivariate stochastic volatility (MSV) modelsMSV models, pioneered by Harvey et al., (1994), assume volatility follows a random process, unlike GARCH which assumes it is a deterministic function of the past., are generally limited to studying processes with less than 5 components (Gouriéroux et al.,, 2009). Recent efforts have led to simpler but less general models, which make assumptions such as constant correlations (Bollerslev,, 1990) – leaving only the diagonal entries of $\Sigma(t)$ to vary.

We hope to unite machine learning and econometrics in an effort to solve these problems. We introduce a stochastic process with Wishart marginals: the generalised Wishart process (GWP). It is a collection of positive semi-definite random matrices indexed by any arbitrary dependent variable $\bm{z}$. We call it the generalised Wishart process, since it is a generalisation of the first Wishart process defined by Bru, (1991). Bru’s Wishart process has recently been used (Gouriéroux et al.,, 2009) in multivariate stochastic volatility (MSV) models (Philipov and Glickman,, 2006; Harvey et al.,, 1994). This prior work on Wishart processes is limited for several reasons: 1) it assumes the dependent variable is a scalar, 2) it is restricted to using an Ornstein-Uhlenbeck covariance structureAn Ornstein-Uhlenbeck process (Uhlenback and Ornstein,, 1930) was first introduced to model the velocity of a particle undergoing Brownian motion. (which means $\Sigma(t+a)$ and $\Sigma(t-a)$ are independent given $\Sigma(t)$, and complex dependencies cannot be captured), 3) it is autoregressive, and 4) there are no general learning and inference procedures. The generalised Wishart process (GWP) addresses all of these issues. Specifically, in the GWP formulation,

The dependent variable can come from any arbitrary index set, just as easily as it can represent time. This allows one to effortlessly condition on covariates like interest rates.

One can easily specify a range of covariance structures (periodic, smooth, Ornstein-Uhlenbeck, …).

We develop Bayesian inference procedures to make predictions, and to learn distributions over any relevant parameters. Aspects of the covariance structure are learned from data, rather than being a fixed property of the model.

Overall, the GWP is versatile and simple. It does not require any free parameters, and any optional parameters are easy to interpret. For this reason, it also scales well with dimension. Yet, the GWP provides an especially general description of multivariate volatility – more so than the most general MGARCH specifications. In the next section, we review Gaussian processes (GPs), which are used to construct the Wishart process. In the following sections we then review the Wishart distribution, present the GWP construction, introduce procedures for inference and predictions, review the main competitor, MGARCH, and present experiments that show how the GWP outperforms MGARCH on simulated and financial data. These experiments include a 5 dimensional data set, based on returns for NASDAQ, FTSE, NIKKEI, TSE, and the Dow Jones Composite, and a set of returns for 3 foreign currency exchanges. In a subsequent version we will also present a 200 dimensional experiment to show how the GWP can be used to study high dimensional problems.

Also, although it is not the focus of this paper, we show in the inference section how the GWP can additionally be used as part of a new GP based regression model that accounts for changing correlations. In other words, it can be used to predict the mean $\bm{\mu}(t)$ together with the covariance matrix $\Sigma(t)$ of a multivariate process. Alternative GP based multivariate regression models for $\bm{\mu}(t)$, which account for fixed correlations, were recently introduced by Bonilla et al., (2008), Teh et al., (2005), and Boyle and Frean, (2004).

Gaussian Processes

We briefly review Gaussian processes, since the generalised Wishart process is constructed from GPs. For more detail, see Rasmussen and Williams, (2006).

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution. Using a Gaussian process, we can define a distribution over functions $u(\bm{z})$:

where $\bm{z}$ is an arbitrary (potentially vector valued) dependent variable, and the mean $m(\bm{z})$ and kernel function $k(\bm{z},\bm{z}^{\prime})$ are respectively defined as

This means that any collection of function values has a joint Gaussian distribution:

where the $N\times N$ Gram matrix $K$ has entries $K_{ij}=k(\bm{z}_{i},\bm{z}_{j})$, and the mean $\bm{\mu}$ has entries $\bm{\mu}_{i}=m(\bm{z}_{i})$. The properties of these functions (smoothness, periodicity, etc.) are determined by the kernel function. The squared exponential kernel is popular:

Functions drawn from a Gaussian process with this kernel function are smooth, and can display long range trends. The length-scale hyperparameter $l$ is easy to interpret: it determines how much the function values $u(\bm{z})$ and $u(\bm{z}+\bm{a})$ depend on one another, for some constant $\bm{a}$.

is an example of a Gaussian process with a fixed covariance structure.

Wishart Distribution

The Wishart distribution defines a probability density function over positive definite matrices $S$:

where $V$ is a $D\times D$ positive definite scale matrix, and $\nu>0$ is the number of degrees of freedom. This distribution has mean $\nu V$ and mode $(D-\nu-1)V$ for $\nu\geq D+1$. $\Gamma_{D}(\cdot)$ is the multivariate gamma function:

The Wishart distribution is a multivariate generalisation of the Gamma distribution when $\nu$ is real valued, and the chi-square ($\chi^{2}$) distribution when $\nu$ is integer valued. The sum of squares of univariate Gaussian random variables is chi-squared distributed. Likewise, the sum of outer products of multivariate Gaussian random variables is Wishart distributed:

where the $\bm{u}_{i}$ are i.i.d. $\mathcal{N}(\bm{0},V)$ $D$-dimensional random variables, and $\mathcal{W}_{D}(V,\nu)$ is a Wishart distribution with $D\times D$ scale matrix $V$, and $\nu$ degrees of freedom. $S$ is a $D\times D$ positive definite matrix. If $D=V=1$ then $\mathcal{W}$ is a chi-square distribution with $\nu$ degrees of freedom. $S^{-1}$ has the inverse Wishart distribution, $\mathcal{W}_{D}^{-1}(V^{-1},\nu)$, which is a conjugate prior for covariance matrices of zero mean Gaussian distributions. This means that for data $\mathcal{D}$ if a prior $p(R)$ is inverse Wishart, and the likelihood $p(\mathcal{D}|R)$ is Gaussian with zero mean, then the posterior $p(R|\mathcal{D})$ is also inverse Wishart.

Generalised Wishart Process Construction

We saw that the Wishart distribution is constructed from multivariate Gaussian distributions. Essentially, by replacing these Gaussian distributions with Gaussian processes, we define a process with Wishart marginals – the generalised Wishart process. It is a collection of positive semi-definite random matrices indexed by any arbitrary (potentially high dimensional) dependent variable $\bm{z}$. For clarity, we assume that time is the dependent variable, even though it takes no more effort to use a vector-valued variable $\bm{z}$ from any arbitrary set. Everything we write would still apply if we replaced $t$ with $\bm{z}$.

Suppose we have $\nu D$ independent Gaussian process functions, $u_{id}(t)\sim\mathcal{GP}(0,k)$, where $i=1,\dots,\nu$ and $d=1,\dots,D$. This means $\text{cov}(u_{id}(t),u_{id}(t^{\prime}))=k(t,t^{\prime})\delta_{ii^{\prime}}\delta_{dd^{\prime}}$, and $(u_{id}(t_{1}),u_{id}(t_{2}),\dots,u_{id}(t_{N}))^{\top}\sim\mathcal{N}(0,K)$, where $\delta_{ij}$ is the Kronecker delta, and $K$ is an $N\times N$ Gram matrix with elements $K_{ij}=k(t_{i},t_{j})$. Let $\hat{\bm{u}}_{i}(t)=(u_{i1}(t),\dots,u_{iD}(t))^{\top}$, and let $L$ be the lower Cholesky decomposition of a $D\times D$ scale matrix $V$, such that $LL^{\top}=V$. Then at each $t$ the covariance matrix $\Sigma(t)$ has a Wishart marginal distribution,

subject to the constraint that the kernel function $k(t,t)=1$.

We write $\Sigma(t)\sim\mathcal{GWP}(V,\nu,k(t,t^{\prime}))$ to mean that $\Sigma(t)$ is a collection of positive semi-definite random matrices with $\mathcal{W}_{D}(V,\nu)$ marginal distributions. Assuming the dependent variable is time, a draw from a Wishart process is a collection of matrices indexed by time (Figure 1), much like a draw from a Gaussian process is a collection of function values indexed by time.

Using this construction, we can also define a generalised inverse Wishart process (GIWP). If $\Sigma(t)\sim\mathcal{GWP}$, then inversion at each value of $t$ defines a draw $R(t)=\Sigma(t)^{-1}$ from the GIWP. The conjugacy of the GIWP with a Gaussian likelihood could be useful when doing Bayesian inference.

We can further extend this construction by replacing the Gaussian processes with copula processes (Wilson and Ghahramani,, 2010). For example, as part of Bayesian inference we could learn a mapping that would transform the Gaussian processes $u_{id}$ to Gaussian copula processes with marginals that better suit the covariance structure of our data set; the result is a Wishart copula process.

The formulation we outlined in this section is different from other multivariate volatility models in that one can specify a kernel function $k(t,t^{\prime})$ that controls how $\Sigma(t)$ varies with $t$ – for example, $k(t,t^{\prime})$ could be periodic – and $t$ need not be time: it can be an arbitrary dependent variable, including covariates like interest rates. In the next section we introduce, for the first time, general inference procedures for making predictions when using a Wishart process prior. These are based on recently developed Markov chain Monte Carlo techniques (Murray et al.,, 2010). We also introduce a new method for doing multivariate GP based regression with dynamic correlations.

Bayesian Inference

Assume we have a generalised Wishart process prior on a dynamic $D\times D$ covariance matrix:

We want to infer the posterior $\Sigma(t)$ given a $D$-dimensional data set $\mathcal{D}=\{\bm{x}(t_{n}):n=1,\dots,N\}$. We explain how to do this for a general likelihood function, $p(\mathcal{D}|\Sigma(t))$, by finding the posterior distributions over the parameters in the model, given the data $\mathcal{D}$. These parameters are: a vector of all relevant GP function values $\bm{u}$, the hyperparameters of the GP kernel function $\bm{\theta}$, the degrees of freedom $\nu$, and $L$, the lower cholesky decomposition of the scale matrix $V$ ($LL^{\top}=V$). The graphical model in Figure 2 shows all the relevant parameters and conditional dependence relationships.

We can sample from these posterior distributions using Gibbs sampling (Geman and Geman,, 1984), a Markov chain Monte Carlo algorithm where initialising $\{\bm{u},\bm{\theta},L,\nu\}$ and then sampling in cycles from

will converge to samples from $p(\bm{u},\bm{\theta},L,\nu|\mathcal{D})$. We will successively describe how to sample from the posterior distributions (14), (15), (16), and (17). In our discussion we assume there are $N$ data points (one at each time step or input), and $D$ dimensions. We then explain how to make predictions of $\Sigma(t_{*})$ at some test input $t_{*}$. Finally, we discuss a potential likelihood function, and how the GWP could also be used as a new GP based model for multivariate regression with outputs that have changing correlations.

In this section we describe how to sample from the posterior distribution (14) over the Gaussian process function values $\bm{u}$. We order the entries of $\bm{u}$ by fixing the degrees of freedom and dimension, and running the time steps from $n=1,\dots,N$. We then increment dimensions, and finally, degrees of freedom. So $\bm{u}$ is a vector of length $ND\nu$. As before, let $K$ be an $N\times N$ Gram matrix, formed by evaluating the kernel function at all pairs of training inputs. Then the prior $p(\bm{u}|\bm{\theta})$ is a Gaussian distribution with $ND\nu\times ND\nu$ block diagonal covariance matrix $K_{B}$, formed using $D\nu$ of the $K$ matrices; if the hyperparameters of the kernel function change depending on dimension or degrees of freedom, then these $K$ matrices will be different from one another. In short,

With this prior, and the likelihood formulated in terms of the other parameters, we can sample from the posterior (14). Sampling from this posterior is difficult, because the Gaussian process function values are highly correlated by the $K_{B}$ matrix. We use Elliptical Slice Sampling (Murray et al.,, 2010): it has no free parameters, jointly updates every element of $\bm{u}$, and was especially designed to sample from posteriors with correlated Gaussian priors. We found it effective.

2 Sampling the other parameters

We can similarly obtain distributions over the other parameters. The priors we use will depend on the data we are modelling. We placed a vague lognormal prior on $\bm{\theta}$ and sampled from the posterior (15) using axis aligned slice sampling if $\bm{\theta}$ was one dimensional, and Metropolis Hastings otherwise. We also used Metropolis Hastings to sample from (16), with a spherical Gaussian prior on the elements of $L$. To sample (17), one can use reversible jump MCMC (Green,, 1995; Robert and Casella,, 2004). But in our experiments we set $\nu=D+1$, and found it effective. Although learning $L$ is not expensive, one might simply wish to set it by taking the empirical covariance of the data set, dividing by the degrees of freedom, and then taking the lower cholesky decomposition.

3 Making predictions

Once we have learned the parameters $\{\bm{u},\bm{\theta},L,\nu\}$, we can find a distribution over $\Sigma(t_{*})$ at a test input $t_{*}$. To do this, we must infer the distribution over $\bm{u}_{*}$ – all the relevant GP function values at $t_{*}$:

Consider the joint distribution over $\bm{u}$ and $\bm{u}_{*}$:

Supposing that $\bm{u}_{*}$ and $\bm{u}$ respectively have $p$ and $q$ elements, then $A$ is a $p\times q$ matrix of covariances between the GP function values $\bm{u}_{*}$ and $\bm{u}$ at all pairs of the training and test inputs: $A_{ij}=k_{i}(t_{*},t_{\text{mod}(N+1,j)})$ if $1+(i-1)N\leq j\leq iN$, and otherwise. The kernel function $k_{i}$ may differ from row to row, if it changes depending on the degree of freedom or dimension; for instance, we could have a different length-scale for each new dimension. $I_{p}$ is a $p\times p$ identity matrix representing the prior independence between the GP function values in $\bm{u}_{*}$. Conditioning on $\bm{u}$, we find

We can then construct $\Sigma(t_{*})$ using equation (12) and the elements of $\bm{u}_{*}$.

4 Likelihood function

So far we have avoided making the likelihood explicit; the inference procedure we described will work with a variety of likelihoods parametrized through a matrix $\Sigma(t)$, such as the multivariate $t$ distribution. However, assuming for simplicity that each of the variables $\bm{x}(t_{n})$ has a Gaussian distribution,

where $\bm{w}(t_{n})=\bm{x}(t_{n})-\bm{\mu}(t_{n})$. We can learn a distribution over $\bm{\mu}(t)$, in addition to $\Sigma(t)$. Here are three possible specifications of $\bm{\mu}$:

In (26), the mean function is directly coupled to the covariance matrix in (12), since they are both constructed using the same Gaussian processes. As the components of $\bm{\mu}$ increase in magnitude, so do the entries in $\Sigma$. This is a desirable property if we expect, for example, high returns to be associated with high volatility. This property is encouraged but not enforced in (27), where a separate vector of Gaussian processes $\bm{u}_{\nu+1}(t)$ is introduced into the expression for the mean function, but not the expression for $\Sigma(t)$. In (28), the mean function is solely this separate vector of Gaussian processes. In each of these cases, we can make mean predictions by inferring distributions over the GP function values $\bm{u}$, as outlined above. And so in each case, the GWP is being used as a GP based regression model which accounts for multiple outputs that have changing correlations.

Alternative models, which account for fixed correlations, have recently been introduced by Bonilla et al., (2008), Teh et al., (2005), and Boyle and Frean, (2004). Rather than use a GP based regression, as in (26)-(28), Gelfand et al., (2004) combine a spatial Wishart process with a parametric linear regression on the mean, to make correlated mean predictions in a spatial setting. Generally their methodology is substantially different: the correlation structure is a fixed property of their method (they do not learn the parameters of a kernel function), and they are not developing a multivariate volatility model, so are not interested in explicitly learning or evaluating the accuracy of the dynamic correlations; in fact, they do not explain how to make predictions of $\Sigma$ at a test input. Further, they do not sample $\bm{\theta},L,$ or $\nu$, and their inference is not explained except that their sampling relies solely on Metropolis Hastings with Gaussian proposals, which will not scale to high dimensions, and will not mix efficiently as the strong GP prior correlations are not accounted for. They fix $L$ as diagonal, which significantly limits the correlation structure of the dynamic covariance matrices (e.g. $\Sigma(t)$), as does fixing $\bm{\theta}$, which we have empirically found to severely affect the quality of predictions. In this paper we focus on making predictions of $\Sigma(t)$, setting $\bm{\mu}=0$.

5 Computational complexity

In contrast to the alternatives, our method scales exceptionally nicely with dimension. MGARCH, the main competitor, is limited to 5 dimensions, at which point severe assumptions – such as constant correlations – are needed for tractability in higher dimensions. We conjecture that other Wishart process models are similarly limited, with Gelfand et al., (2004) restricted to about 2 dimensions.

Our method is mainly limited by taking the cholesky decomposition of the block diagonal $K_{B}$, a $ND\nu\times ND\nu$ matrix. However, $\text{chol}(\text{blkdiag}(A,B,\dots))=\text{blkdiag}(\text{chol}(A),\text{chol}(B),\dots)$. So in the case with equal length-scales for each dimension, we only need to take the cholesky of an $N\times N$ matrix $K$, an $\mathcal{O}(N^{3})$ operation, independent of dimension! In the more general case with $D$ different length-scales, it is an $\mathcal{O}(DN^{3})$ operation. Taking into account the likelihood, and other operations, the total training complexity is either $\mathcal{O}(N^{3}+\nu D^{2})$ for equal length-scales, or $\mathcal{O}(DN^{3}+\nu D^{2})$ using separate length-scales for each dimension. In the latter more general case, we could in principle go to about 1000 dimensions, assuming $\nu$ is $\mathcal{O}(D)$, and for instance, a couple years worth of financial training data, which is typical for GARCH (Brownlees et al.,, 2009). Using sparse GP techniques we could go to even higher dimensions. In practice, MCMC may be infeasible for very high $D$, but we have found Elliptical Slice Sampling incredibly robust: we will shortly include a 200 dimensional experiment in an updated version. Overall, this is an impressive scaling – without making further assumptions in our model, we can go well beyond 5 dimensions with full generality.

Multivariate GARCH

We compare predictions of $\Sigma(t)$ made by the generalised Wishart process to those made by multivariate GARCH (MGARCH), since GARCH (Bollerslev,, 1986) is extremely popular and arguably unsurpassed at predicting the volatility of returns on equity indices and currency exchanges (Poon and Granger,, 2005; Hansen and Lunde,, 2005; Brownlees et al.,, 2009).

Consider a zero mean $D$ dimensional vector stochastic process $\bm{x}(t)$ with a time changing covariance matrix $\Sigma(t)$ as in (24). In the general MGARCH framework,

The first and most general MGARCH model, the VEC model of Bollerslev et al., (1988), specifies $\Sigma_{t}$ as

$A_{i}$ and $B_{j}$ are $D(D+1)/2\times D(D+1)/2$ matrices of parameters, and $\bm{a}_{0}$ is a $D(D+1)/2\times 1$ vector of parameters.We use $\Sigma(t)$ and $\Sigma_{t}$ interchangeably. The vech operator stacks the columns of the lower triangular part of a $D\times D$ matrix into a vector of length $D(D+1)/2$. For example, $\text{vech}(\Sigma)=(\Sigma_{11},\Sigma_{21},\dots,\Sigma_{D1},\Sigma_{22},\dots,\Sigma_{D2},\dots,\Sigma_{DD})^{\top}$. This model is general, but difficult to use. There are $(p+q)(D(D+1)/2)^{2}+D(D+1)/2$ parameters! These parameters are hard to interpret, and there are no conditions under which $\Sigma_{t}$ is positive definite for all $t$. Gouriéroux, (1997) discusses the challenging (and sometimes impossible) problem of keeping $\Sigma_{t}$ positive definite. Training is done by a constrained maximum likelihood, where the log likelihood is given by

supposing that $\bm{\eta}_{t}\sim\mathcal{N}(0,I)$, and that there are $N$ training points.

Subsequent efforts have led to simpler but less general models. We can let $A_{j}$ and $B_{j}$ be diagonal matrices. This model has notably fewer (though still $(p+q+1)D(D+1)/2$) parameters, and there are conditions under which $\Sigma_{t}$ is positive definite for all $t$ (Engle et al.,, 1994). But now there are no interactions between the different conditional variances and covariances. A popular variant assumes constant correlations between the $D$ components of $\bm{x}$, and only lets the marginal variances – the diagonal entries of $\Sigma(t)$ – vary (Bollerslev,, 1990).

We compare to the ‘full’ BEKK variant of Engle and Kroner, (1995), as implemented by Kevin Shepphard in the UCSD GARCH Toolbox.http://www.kevinsheppard.com/wiki/UCSD_GARCH We chose BEKK because it is the most general MGARCH variant in widespread use. We use the first order model:

where $A,B$ and $C$ are $D\times D$ matrices of parameters. $C$ is lower triangular to ensure that $\Sigma_{t}$ is positive definite during maximum likelihood training. For a full review of multivariate GARCH models, see Silvennoinen and Teräsvirta, (2009).

Experiments

To make these predictions we learn distributions over the GWP parameters through the Gibbs sampling procedure outlined in section 5. The kernel functions we use are solely parametrized by a one dimensional length-scale $l$, which indicates how dependent $\Sigma(t)$ and $\Sigma(t+a)$ are on one another. We place a lognormal prior on the length-scale, and sample from the posterior with axis-aligned slice sampling.

For each experiment, we choose a kernel function we want to use with the GWP. We then compare to a GWP that uses an Ornstein-Uhlenbeck (OU) kernel function, $k(t,t^{\prime})=\exp(-|t-t^{\prime}|/l)$. Even though we are still taking advantage of the inference procedures in the GWP formulation, we refer to this variant of GWP as a simple Wishart process (WP), since the classic Bru, (1991) construction is like a special case of our generalised Wishart process restricted to using a one dimensional Gaussian process with an OU covariance structure.

We begin by generating a $2\times 2$ time varying covariance matrix $\Sigma_{p}(t)$ with periodic components, and simulating data at 291 time steps from a Gaussian distribution:

Periodicity is especially common to financial and climate data, where daily trends repeat themselves. For example, the intraday volatility on equity indices and currency exchanges has a periodic covariance structure. Andersen and Bollerslev, (1997) discuss the lack of – and critical need for – models that account for this periodicity. In the GWP formulation, we can easily account for this by using a periodic kernel function. We reconstruct $\Sigma_{p}(t)$ using the kernel $k(t,t^{\prime})=\exp(-2\sin((t-t^{\prime})^{2})/l^{2})$. We reconstructed the historical $\Sigma_{p}$ at all 291 data points, and after having learned the parameters for each of the models from the first 200 data points, made one step forecasts for the last 91 points. Table 1 and Figure 3 show the results. We call this data set PERIODIC. The GWP outperforms the competition on all error measures. It identifies the periodicity and underlying smoothness of $\Sigma_{p}$ that neither the WP nor MGARCH accurately discern: both are too erratic. And MGARCH is especially poor at learning the time changing covariance (off-diagonal entry of $\Sigma_{p}$) in this data set.

For our next experiment, we predict $\Sigma(t)$ for the returns on three currency exchanges – the Canadian to US Dollar, the Euro to US Dollar, and the US Dollar to the Great Britain Pound – in the period 15/7/2008-15/2/2010; this encompasses the recent financial crisis and so is of particular interest to economists. We call this data set EXCHANGE. We use the proxy $S_{ij}(t)=x_{i}(t)x_{j}(t)$. With the GWP, we use the squared exponential kernel $k(t,t^{\prime})=\exp(-0.5(t-t^{\prime})^{2}/l^{2})$. We make 200 one step ahead forecasts, having learned the parameters for each of the models on the previous 200 data points. We also make 200 historical predictions for the same data points as the forecasts. The results are in Table 1.

To make forecasts and historical predictions on this data set (EQUITY), we used a GWP with a squared exponential kernel, $k(t,t^{\prime})=\exp(-0.5(t-t^{\prime})^{2}/l^{2})$. We follow the same procedure as before and make 200 forecasts and historical predictions; results are in Table 1.

Both the EXCHANGE and EQUITY data sets are especially suited to GARCH (Poon and Granger,, 2005; Hansen and Lunde,, 2005; Brownlees et al.,, 2009; Bollerslev and Ghysels,, 1996; McCullough and Renfro,, 1998; Brooks et al.,, 2001). However, the generalised Wishart process outperforms GARCH on both of these data sets. Based on our experiments, there is evidence that the GWP is particularly good at capturing the co-variances (off-diagonal elements of $\Sigma(t)$) as compared to GARCH. The GWP also outperforms the WP, which has a fixed Ornstein-Uhlenbeck covariance structure, even though in our experiments the WP takes advantage of the new inference procedures we have derived. Thus the difference in performance is likely because the GWP is capable of capturing complex interdependencies in volatility, whereas the WP is not.

Discussion

We introduced a stochastic process – the generalised Wishart process (GWP) – which we used to model time-varying covariance matrices $\Sigma(t)$. Unlike the alternatives, the GWP can easily model a diverse class of covariance structures for $\Sigma(t)$. In the future, the GWP could be applied to study how these covariance matrices depend on other variables, like interest rates, in addition to time. Future research could also apply the GWP to extremely high dimensional problems.

We hope to unify efforts in machine learning and econometrics to inspire new multivariate volatility models that are simultaneously general, easy to interpret, and tractable in high dimensions.

Thanks to Carl Edward Rasmussen and John Patrick Cunningham for helpful discussions. AGW is supported by an NSERC grant.

References