How to simulate cointegrated prices

Question

Is there any simple way to simulate cointegrated prices?

Answer 1

Consider a $T \times N$ matrix of potentially cointegrating prices $P$. Define $Y_{t}\equiv ln\left(P_{t}\right)$. In the multivariate framework, there are two basic methods to estimate the cointegrating relationships. The first is an error correction framework of the form $$\Delta Y_{t} = \beta_{0}+\beta_{1}\Delta Y_{t-1}+\beta_{2}Y_{t-1}+\varepsilon_{t}$$ that is most convenient when attempting to perform statistical tests on the coefficients. The alternate approach is a vector autoregressive model of the form $$Y_{t} = \beta_{0}+\beta_{1}Y_{t-1}+\varepsilon_{t}.$$ For the purposes of simulation, they are effectively equivalent. One must estimate $\beta_{0}$ and $\beta_{1}$ and solve for $\varepsilon_{t}$. There are many potential distribution assumptions that one could make about the behavior of $\varepsilon_{t}$, but a simple one would be that it follows a multivariate normal distribution with a mean of zero and a covariance matrix equal to the sample covariance matrix. More complicated assumptions might be that the variances and correlations are time-varying or that there are fat tails. For financial time series, these may be important to consider.

To simulate $\widetilde{Y}_{t+1}$, you thus obtain $\widetilde{\varepsilon}_{t+1}$ by whatever means are appropriate and calculate $$\widetilde{Y}_{t+1} = \beta_{0}+\beta_{1}Y_{t}+\widetilde{\varepsilon}_{t+1}$$

For $i>1$, one would need to be careful to incorporate the simulated values from the previous period so that $$\widetilde{Y}_{t+i} = \beta_{0}+\beta_{1}\widetilde{Y}_{t+i-1}+\widetilde{\varepsilon}_{t+i}$$ in order to ensure the autoregressive features in each simulated path.

After calculating the simulated values of $\widetilde{Y}_{t+i}$, one would want to convert them back to prices by calculating $\widetilde{P}_{t+i}\equiv \mathrm{exp}\left(\widetilde{Y}_{t+i}\right)$.

Answer 2

One way to construct cointegrated timeseries it to use the error-correction representation (see Engle, Granger 1987 for details of the equivalence).

To generate two timeseries that are cointegrated, start with your cointegrating vector $(\alpha_1, \alpha_2)$ so that you want $\alpha_1x_t + \alpha_2y_t$ to be stationary; choose initial values $x_0, y_0$ and a parameter $\gamma\in (0,1)$ that controls how strongly cointegrated the series are. Then generate each timestep as:

$x_{t+1} = x_t - \gamma (x_t + (\alpha_2/\alpha_1)y_t) + \epsilon_{1t}$

$y_{t+1} = y_t - \gamma (y_t + (\alpha_1/\alpha_2)x_t) + \epsilon_{2t}$

For price series, it's generally the cumulative returns that you want to be cointegrated. To generate prices, as John mentioned in his comment above, follow the above procedure for log-prices, then exponentiate.

Answer 3

I don't think the accepted answer fully answers the question. In particular simulating a non-stationary VAR is not guaranteed to give you a cointegrated system, and we may be interested in simulating $n>2$ time series. To be specific, the question is asking how to simulate a vector time series $Y_t \in \mathbb{R}^n$ with the following properties:

1. $Y_t \sim I(1)$ - the time series is integrated of order one.
2. There exists some (non-zero) cointegrating vector $\boldsymbol{a} \in \mathbb{R}^n$ s.t. $\boldsymbol{a}^TY_t \sim I(0)$, and $\boldsymbol{a}$ is unique up to a constant.

Here is a strategy based on simulating a non-stationary VAR.

VAR representation

A cointegrated system will admit a VAR(p) representation:

$$Y_t = \boldsymbol{\mu} + \Phi_1 Y_{t-1} + ... + \Phi_1 Y_{t-p} + \boldsymbol{\epsilon}_t$$

$$\Leftrightarrow \Phi(L) Y_t = \boldsymbol{\mu + \epsilon_t}$$

With $n$ time series there may be up to $n-1$ linearly independent cointegrating vectors $\boldsymbol{a}_1,...,\boldsymbol{a}_h$ which together form a basis for the space of cointegrating vectors that needs to be specified as part of the simulation. In particular letting:

$$\boldsymbol{A}^T = \begin{bmatrix} \boldsymbol{a}_1^T\\ \vdots \\ \boldsymbol{a}_h^T \end{bmatrix}$$

Then $\boldsymbol{A}^TY_t \sim I(0)$ is a stationary vector. A bit of algebra (not shown here for convenience, see 19.1 in Hamilton) gives that every row of $\Phi(1)$ is a cointegrating vector and so is part of the space spanned by $(\boldsymbol{a}_1,...,\boldsymbol{a}_h)$.

In particular there must be some $n \times h$ matrix $\boldsymbol{B}$ such that the VAR coefficients satisfy (i) $\Phi(1)=\boldsymbol{BA}^T$, and (ii) $|\Phi(z)|$ has exactly $n-h$ unit roots.

Simulation strategy

Based on the above a simulation strategy using a non-stationary VAR(1) is as follows:

1. Choose $n$ the dimension of the time series to simulate, and $h \leq n-1$ the number of cointegrating relationships.
2. Generate an $n \times h$ orthonormal matrix $\boldsymbol{A}$.
3. Set $\Phi = \boldsymbol{I} - \boldsymbol{AA}^T$
4. Simulate $Y_t = \boldsymbol{\mu} + \Phi_1 Y_{t-1} + \boldsymbol{\epsilon}_t $ using standard techniques.

In R random orthonormal matrices can be generated with the randortho command in the pracma package. This strategy can be extended to general VAR(p), but care has to be taken to ensure that the matrix of VAR coefficients in companion form does not have any eigenvalues greater than 1.

Notice also that the cointegration property is completely independent of the covariance structure of the innovations, we need only that $\boldsymbol{\epsilon}_t$ be white noise.

R code

Here is some R code for generating cointegrated time series with arbitrary $n$ and $h$. The code is used to generate 10 time series with three cointegrating relations.

library(docstring) # document function
library(pracma) # random orthonormal
library(mvtnorm) # multivariate normal

cointegrated.vector <- function(tt,n,h,C,burn = 1000){

  #'@param tt int, time points to simulate
  #'@param n int, number of time series
  #'@param h int, number of cointegrating relations
  #'@param C innovation covariance structure of innovations
  #'@param burn

  tot.t <- tt+burn
  innov <- rmvnorm(tot.t,rep(0,n),C)
  X <- matrix(0,1+tot.t,n)

  A <- randortho(n)[1:h,]
  if (h == 1) Phi <- diag(n) - outer(A,A)
  else Phi <- diag(n) - t(A)%*%A

  for (i in 1:tot.t){
    X[i+1,] <- Phi%*%X[i,] + innov[i,]
  }
  X[-burn,]
}

plot.ts(cointegrated.vector(100,10,3,diag(10)))