# How to simulate two correlated time series

My question is similar to How to simulate two correlated AR(1) time series? but more complex as I want to simulate time series with various AR and MA parts. I have one time series with the historic online sales of a certain product and a second time series that shows the number of times the websites of that product was opened (i.e., site visits) in the same time period. Naturally, the historic sales and site visits correlate. I now want to simulate two new time series based on this. If i simulate them separately, I lose the information that they are correlated. So how do I simulate this, taking into account correlation? The answers on How to simulate two correlated AR(1) time series? only partially help as they assume an AR(1) process. My historic sales time series is, however, an ARIMA(0,1,1) process and the site visits an ARIMA(3,1,1) with drift process (when using (auto.arima)).

The data is provided below (each time series consists of 100 observations). I would very much appreciate any answer that outlines how to code this in R. In additional challenge is that the simulated data cannot have negative values.

Data$Sales [1] 2 1 2 3 0 1 0 6 6 2 0 0 0 1 1 0 2 0 2 1 3 0 2 0 0 0 3 1 0 1 0 1 4 1 4 0 2 2 1 0 1 0 1 0 0 0 1 1 0 0 0 0 [53] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 Data$Site visits
[1] 1 2 4 1 1 3 0 5 3 0 0 0 2 4 0 1 1 3 4 1 6 0 0 0 1 0 2 2 2 0 1 4 3 3 3 1 1 2 3 2 0 1 2 2 0 0 3 2 0 0 0 2
[53] 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0

as.numeric(cor(Data$$Sales,Data$$Site visits))
[1] 0.6155173
$$$$
`
• The information in your question suggests you aren't using appropriate models. It looks like you need a GLM version: that is, you can model an underlying bivariate process $(Y_1(t),Y_2(t))$ in any way you like and that process determines the observations through two link functions $h_i$ and a response distribution, such as a Poisson$(\lambda)$ distribution, so that what you observe is $(X_1(t),X_2(t))$ where $X_i(t)\sim\text{Poisson}(h_i^{-1}(Y_i(t))).$ Often $Y(t)$ is modeled as a Gaussian process. – whuber Sep 14 '20 at 14:57
• to follow @whuber comment, you need a fully specified joint model to generate from it. – Xi'an Sep 14 '20 at 15:01
• I unfortunately have no statistical or mathematical background so I honestly don't know how I should transfer any of this into an R code. While ARIMA might not even be the most appropriate model for this, i have to do the simulation for 3000 different products, for which I cannot separately assess what works best, which is why I simply use auto.arima. The other post that I linked to outlines very nicely how to do the simulation in R for an AR(1) process, so I'm looking for an answer in a similar fashion – Jennifer Weingarten Sep 14 '20 at 15:03