Numerical Python code to generate artificial data from a time series process

Question

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)
    print(data)

if __name__ == "__main__":
    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).

Answer 1

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.

Answer 2

• 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).

Answer 3

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)

---

2. 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.

3. 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.

4. 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.

5. 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.