Numerical Python code to generate artificial data from a time series process 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}
&amp;= 
\begin{bmatrix}
0.02 \\
0.03
\end{bmatrix}
+
\begin{bmatrix}
0.5 &amp; 0.1 \\
0.4 &amp; 0.5
\end{bmatrix}
\begin{bmatrix}
y_{1,t-1} \\
y_{2,t-1} 
\end{bmatrix}
+
\begin{bmatrix}
0 &amp; 0 \\
0.25 &amp; 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 &lt; 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).
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).