Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I get very poor replication of longitudinal parameters from my own program using the Box-Jenkins model. I had no such problem with my own program generating AR(1) Gaussian data. Is there some trick that I haven't discovered yet? I would like to use the validation framework of Martinez-Rivera & Ventosa-Santaularia (2012) to check a new time series regression method that is better than the usual on AR(1), but itself depends on ARMA random generation when applied to ARMA data.

share|improve this question
It is hard to say what the problem is without seeing what you did. ARMA(2,2) is just given by X(t) =r$_1$ X(t-1) +r$_2$ X(t-2) + e(t) + a$_1$ e(t-1) +a$_2$ e(t-2). Choose initial values for X(0) and X(1) generate e(t) for t=0,1,2,3,... using a normal random number generator with mean 0 and a fixed variance s$^2$. Then just use the formula to recursively generate the X(t)s. – Michael Chernick Sep 17 '12 at 17:24
Please show your code; that way someone can see what (if anything) you did wrong. – Peter Flom Sep 29 '12 at 19:11

To generate an ARMA(2,2) with specified coefficients just pcik two initial values for x say $x_0$ and $x_1$ and generate as many e(t) as you need from a $N(0, \sigma^2)$ distribution.

Then for each t, $$X(t) =r_1 X(t-1) +r_2 X(t-2) + e(t) + a_1 e(t-1) +a_2 e(t-2).$$ You have $X(1)$ specified as $x_1$ and $X(0)$ specified as $x_0$ take $e(2)$, $e(1)$ and $e(0)$ along with $x_1$ and $x_0$ and plug them into the equation for $X(2)$ and then continue for $X(3)$, $X(4)$,... using the recursion and the new $e(t)$ term.

There should be no problems.

share|improve this answer
How to handle negative indexes in this equation? Etc., t = 0 and then you have t-1=-1. – zygimantus Jan 10 at 14:54

As Michael said , there should be no problems. But I would add two caveats 1) discard the first couple simulations as your starting values can have an impact. I would discard the first 500 values.

2) Make sure that what you are simulating is an invertible model. Check the roots to insure they meet the invertibility requirements. Typical output from some very bad analytical engines deliver non-invertible solutions as their "automatically selected final model" and never report this to their unsuspecting users. Often errors of omission are worse than errors of comission !

share|improve this answer
Really?? I've never seen any automated model selection algorithm that would return a non-invertible model. – Rob Hyndman Sep 18 '12 at 1:43

In addition to the useful answers given, here's some Python code that I wrote that generates an $\text{ARMA}(p,q)$ gaussian time series:

Random generation of Gaussian ARMA(p,q) time series.


phi:      An array of length p with the AR coefficients (the AR part of 
          the ARMA model).

theta:    An array of length q with the MA coefficients (the MA part of 
          the ARMA model).

sigma:    Standard deviaton of the Gaussian noise.

n:        Length of the returned time-series.

burnin:   Number of datapoints that are going to be discarded (the higher 
          the better) to avoid dependence of the ARMA time-series on the 
          initial values.

from numpy import append,array
from numpy.random import normal
def ARMAgenerator(phi,theta,sigma,n,burnin=0,verbose=0):
      burnin=10*l # Burn-in elements!
    for i in range(n+burnin):
          for j in range(len(phi)):
          for j in range(len(theta)):
      print 'Measured standard deviation: '+str(sqrt(var(w[burnin:])))
    return ARMA[burnin:]

An example of the usage of the code: say you want to simulate $n=100$ datapoints of an $\text{ARMA}(p,q)$ model with AR coefficients $\phi=(0.4,0.3)$ and MA coefficients $\theta=(0.1,-0.3)$, with a zero-mean gaussian noise with $\sigma=2$. Also, say you want to simulate $\text{burnin}=500$ datapoints first in order to avoid dependencies on the initial values. You create it with the code above as follows:


And now let's plot and see! (You need to install the matplotlib library in order to do the following):

import matplotlib.pyplot as plt
plt.xlabel("Time units")
plt.ylabel("ARMA(2,2) Series")
plt.plot(x,'-r') # This plots the series as a red line.
plt.plot(x,'or') # This plots the series as red points.


ARMA(2,2) series generated with the Python code given above

share|improve this answer

AR(P) models are straightforward, you use P lags of previous values of the dependent variable for the forecast. The trick with models which have MA(Q) terms is how to seed the lagged Q error terms.

The simplest way is to get the Q error terms from the estimation routine, then supply them to your forecasting function. Alternatively, your forecasting package will most likely use model fit parameters and the P lagged dependent and independent variables to estimate the lagged error terms itself.

See the description of presampled error innovations parameter description in MATLAB arima forecast function here.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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