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

Question

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.

Answer 1

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 !

Answer 2

To generate an ARMA(2,2) with specified coefficents just pcik two initial values for x say x$_)$ and x$_1$ and generate as many e(t) as you need from a N(0, σ$^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.

Answer 3

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

And...tada!