# Random generation of ARMA(2,2) Gaussian time series

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.

-
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

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 !

-
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

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.

-

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.

INPUTS

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):
l=max(len(phi),len(theta))
if(burnin==0):
burnin=10*l # Burn-in elements!
w=normal(0,sigma,n+burnin)
ARMA=array([])
s=0.0
l=max(len(phi),len(theta))
for i in range(n+burnin):
if(i<l):
ARMA=append(ARMA,w[i])
else:
s=0.0
for j in range(len(phi)):
s=s+phi[j]*ARMA[i-j-1]
for j in range(len(theta)):
s=s+theta[j]*w[i-j-1]
ARMA=append(ARMA,s+w[i])
if(verbose!=0):
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:

phi=[0.4,0.3]
theta=[0.1,-0.3]
sigma=2.0
n=100
burnin=500
x=ARMAgenerator(phi,theta,sigma,n,burnin)


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.
plt.show()


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