Constucting ARMA realizations

Given estimates for the AR and MA coefficients for an ARMA(p,q) model with a non zero mean, how do I construct a realization of the time series?

Basis for doing this: I have a time series. I fit an ARMA(p,q) model to it using an IMSL function call and got AR and MA estimates. Now I wish to plot a realization of the ARMA model over my original time series to identify any discrepancies.

How is the time series initialized?
Do I need to generate the normal random errors/ shocks?
If so how do I know the variance of these error?

I need lots of detail please....... HELP!