I'm writing code to generate artificial data from a bivariate time series process, i.e. a vector autoregression. This is my first foray into numerical Python, and it seemed like a good place to start.

The specification is of this form:

\begin{align} \begin{bmatrix} y_{1,t} \\ y_{2,t} \end{bmatrix} &= \begin{bmatrix} 0.02 \\ 0.03 \end{bmatrix} + \begin{bmatrix} 0.5 & 0.1 \\ 0.4 & 0.5 \end{bmatrix} \begin{bmatrix} y_{1,t-1} \\ y_{2,t-1} \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0.25 & 0 \end{bmatrix} \begin{bmatrix} y_{1,t-2} \\ y_{2,t-2} \end{bmatrix} + \begin{bmatrix} u_{1,t} \\ u_{2,t} \end{bmatrix} \end{align}

where \$u\$ has a multivariate normal distribution, i.e. \begin{align} \begin{bmatrix} u_{1,t} \\ u_{2,t} \end{bmatrix} \sim N \left( \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}, \Sigma \right) \end{align}

My basic algorithm generates enough observations to get past the effects of the (arbitrary) initial conditions, and only returns the number of observations asked for. In this case, it generates ten times the number of observations asked for and discards the first 90%.

import scipy

def generate_data(y0, A0, A1, A2, mu, sigma, T):
    """ Generate sample data from VAR(2) process with multivariate
        normal errors

    :param y0: Vector of initial conditions
    :param A0: Vector of constant terms
    :param A1: Array of coefficients on first lags
    :param A2: Array of coefficients on second lags
    :param mu: Vector of means for error terms
    :param sigma: Covariance matrix for error terms
    :param T: Number of observations to generate

    if y0.ndim != 1:
        raise ValueError("Vector of initial conditions must be 1 dimensional")

    K = y0.size

    if A0.ndim != 1 or A0.shape != (K,):
        raise ValueError("Vector of constant coefficients must be 1 dimensional and comformable")

    if A1.shape != (K, K) or A2.shape != (K, K):
        raise ValueError("Coefficient matrices must be conformable")

    if mu.shape != (K,):
        raise ValueError("Means of error distribution must be conformable")

    if sigma.shape != (K, K):
        raise ValueError("Covariance matrix of error distribution must be conformable")

    if T < 3:
        raise ValueError("Cannot generate less than 3 observations")

    N = 10*T

    errors = scipy.random.multivariate_normal(mu, sigma, size = N-1)
    data = scipy.zeros((N, 2))
    data[0, :] = y0
    data[1, :] = y0
    for i in range(2, N):
        data[i, :] = A0 + scipy.dot(A1, data[i-1, :]) + scipy.dot(A2, data[i-1, :]) + errors[i-1,:]
    return(data[-T:, :])

def main():
    y0 = scipy.array([0, 0])
    A0 = scipy.array([0.2, 0.3])
    A1 = scipy.array([[0.5, 0.1], [0.4, 0.5]])
    A2 = scipy.array([[0, 0], [0.25, 0]])
    mu = scipy.array([0, 0])
    sigma = scipy.array([[0.09, 0], [0, 0.04]])
    T = 30

    data = generate_data(y0, A0, A1, A2, mu, sigma, T)

if __name__ == "__main__":

As I said, this code is very simple, but since it's my first attempt at a numerical Python project, I'm looking for any feedback before I slowly start building a larger code project that incorporates this. Does any of my code (base Python, SciPy, etc.) need improvement?

This specific model comes from:

Brüggemann, Ralf, and Helmut Lütkepohl. "Lag Selection in Subset VAR Models with an Application to a US Monetary System." (2000).

  • \$\begingroup\$ I don't understand why you pass in size = N-1 to scipy.random.multivariate_normal instead of size=N when defining errors, and relatedly, I also don't understand why errors[i-1,:] is added to data[i, :] intead of errors[i, :]. Most likely I'm missing something, but is that code consistent with the spec? \$\endgroup\$
    – Curt F.
    Aug 16, 2015 at 4:14

3 Answers 3


If you are going to generate 10 times the amount of data needed, there's no need to keep it around. You can allocate the array you intend to return, and overwrite it several times. Similarly, you don't need to precompute and store all the error terms, but can generate them on the fly as you need them:

    data = scipy.zeros((T, 2))
    rmvn = scipy.random.multivariate_normal
    prev_1 = prev_2 = y0
    for samples in range(10*T):
        idx = samples % T
        data[idx, :] =  A0 + A1.dot(prev_1) + A2.dot(prev_2) + rmvn(mu, sigma)
        prev_1, prev_2 = data[idx, :], prev_1
    return data

The above approach, holding the previous entries in auxiliary variables, to avoid indexing complications, also makes it kind of obvious that your requirement for T >= 3 is easy to work around.

I have also used the .dot method, rather than the dot function, as I believe it makes the long expression a little more readable.

  • \$\begingroup\$ This is great advice, since you're right that I don't need to store all the previously calculated observations. The `.dot' method does make the code easier to read, too. \$\endgroup\$
    – Michael A
    Aug 25, 2015 at 19:59
  • Is

    scipy.dot(A1, data[i-1, :]) + scipy.dot(A2, data[i-1, :])

    a bug or a copy-paste error? According to the spec, the second term should have i - 2.

  • A 90% lead-in seems arbitrary. In fact, any lead-in would be arbitrary, as long as you didn't quantify how much the initial condition should fade. In any case, I recommend to pass it as a parameter (either percentage to discard, or - preferably - the desired decay of the initial condition).

  • \$\begingroup\$ Thanks; your first bullet point is a bug that I missed when I copy-pasted. I corrected that in the original question, and I'll think some more about the second point (because off the top of my head, I'm not exactly sure how to calculate the number of observations needed to decay the initial conditions by a certain amount. It's easy for a univariate AR(1) or something like that, but I'm not sure for a bivariate AR(2)). \$\endgroup\$
    – Michael A
    Aug 14, 2015 at 19:32
  • 1
    \$\begingroup\$ It is against the CR charter to edit the code after the review, because it invalidates the answer. Please roll the edit back. \$\endgroup\$
    – vnp
    Aug 14, 2015 at 19:44
  • \$\begingroup\$ I rolled back the edit; thanks for the heads up. Can you point me to a source that talks about the decay of initial conditions? Thinking more about this, I'm not sure how to interpret it. Is the decay of the initial condition the same as saying that the expected value of \$y_t\$ differs from the initial condition by some percentage, or what? \$\endgroup\$
    – Michael A
    Aug 25, 2015 at 19:57
  • \$\begingroup\$ This is a good answer, but without additional references, clarification, or even any comments on the decay of the initial conditions, I can't accept it. I asked for some references here, if you have any to contribute. \$\endgroup\$
    – Michael A
    Sep 4, 2015 at 21:43

I think this is very good code! A few small recommendations:

  1. Consider wrapping up your parameters into dictionaries. For your code, this is really a matter of style rather than best practices or anything. Your generate_data takes in 7 parameters, which is getting close to making it hard to use your function, mostly because remembering the order to pass those parameters is starting to get tough.

One way to do it:

def generate_data(initial_conditions, model_params, error_params, n_points=30):
     Generate data from a VAR(2) process with multivariate normal errors
    y0 = initial_conditions['y0']
    A0 = model_params['A0']
    A1 = model_params['A1']
    A2 = model_params['A2']
    mu = error_params['mu']
    sigma = error_params['sigma']
    # ...
    # the rest of your code...

def main():
    initial_conditions = {'y0': scipy.array([0, 0])}
    model_params = {'A0': scipy.array([0.2, 0.3]),
                    'A1': scipy.array([[0.5, 0.1], [0.4, 0.5]]),
                    'A2': scipy.array([[0, 0], [0.25, 0]]),
    error_params = {'mu': scipy.array([0, 0]),
                    'sigma': scipy.array([[0.09, 0], [0, 0.04]]),

    data = generate_data(initial_conditions, model_params, error_params)

  1. Defining default parameter values is often a good idea. Above, rather than set T to 30 in main() I just set the default number of points in the def for generate_data.

  2. I also renamed that variable because (at least to me) n_points is more descriptive than T, but admittedly I don't work very heavily with time series.

  3. The error structure you are using right now has no cross-correlations, i.e. the sigma matrix is diagonal. If that will always be the case, it could be much faster to avoid multivariate_normal() in favor of doing separate, 1d calls to normal() instead. Of course, if you will later use this code with non-diagonal sigma matrices, keeping it the way it is is probably better.

  4. I agree with the other answers that you shouldn't keep the data you won't be returning, and also that you should pass in some kind of parameter to offer control over how many points to compute before starting to "keep" the data.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.