sarima.sim <- function (
    n = 20, 
    period = 12, 
    model, 
    seasonal
) {    
  x <- arima.sim( model, n*period )
  x <- x[1:(n*period)]
  for (i in 1:period) {
    xx <- arima.sim( seasonal, n )
    xx <- xx[1:n]
    x[i + period * 0:(n-1)] <- 
      x[i + period * 0:(n-1)] + xx
  }
  x <- ts(x, frequency=period)
  x
}
#op <- par(mfrow=c(3,1))
#x <- sarima.sim(
#  20, 
#  12, 
#  list(ar=.6, ma=.3, order=c(1,0,1)),
#  list(ar=c(.5), ma=c(1,2), order=c(1,0,2))
#)
#eda.ts(x, bands=T)