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.
 A: 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.
A: 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!

A: 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 !
A: 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.
A: Use the rGARMA function in the ts.extend package
You can generate random vectors from any stationary Gaussian ARMA model using the ts.extend package.  This package generates random vectors directly form the multivariate normal distribution using the computed autocorrelation matrix for the random vector, so it gives random vectors from the exact distribution and does not require "burn-in" iterations.  Here is an example from an $\text{ARMA}(2,2)$ model.
#Load the package
library(ts.extend)

#Set parameters
AR    <- c(0.4,  0.3)
MA    <- c(0.1, -0.3)
SIGMA <- 2
m     <- 100

#Generate n = 12 random vectors from this model
set.seed(1)
SERIES1 <- rGARMA(n = 12, m = m, 
                  ar = AR, ma = MA, mean = 0, errorvar = SIGMA^2)

#Plot the series using ggplot2 graphics
library(ggplot2)
plot(SERIES1)


Another nice thing about this function is that you can also generate values from the conditional distribution where you specify some of the values in the vector.  Again, this give rando vectors from the exact conditional multivariate normal distribution; it does not require "burn-in" iterations.  Here is an example where we specify some conditioning values.
#Specify vector of conditioning values (y_1 = 2, y_34 = -2, y_90 = 3)
CONDVALS     <- rep(NA, m)
CONDVALS[1]  <- 2
CONDVALS[34] <- -2
CONDVALS[90] <- 3

#Generate n = 12 random vectors from this conditional model
set.seed(1)
SERIES2 <- rGARMA(n = 12, m = m, condvals = CONDVALS, 
                  ar = AR, ma = MA, mean = 0, errorvar = SIGMA^2)

#Plot the series using ggplot2 graphics
plot(SERIES2)


