function [x,phir,thetar] = armagen(phi,theta,n,sigma,x0)
%                   GENERATOR OF A GAUSSIAN ARMA(p,q) PROCESS
% function [x,phir,thetar] = armagen(phi,theta,n,sigma,x0)
% phi   ... AR coefficients, p=length(phi), default=[] <=> p=0
% theta ... MA coefficients, q=length(theta), default=[] <=> q=0
%           p=0 => x=MA(q) process, q=0 => x=AR(p) process
%           p=0 & q=0 => x=gaussian zero-mean white noise WN(0,sigma^2).
% n     ... number of samples to generate, default=[] <=> n=100
% sigma ... standard deviation of the generating white noise process,
%           default=[] <=> sigma=1.
% x0    ... initial data sequence relevant to AR only (p>0),
%           default=[], n0=length(x0).
%           n0<p => the leading part of x0 updated with p-n0 zero samples
%           n0>p => only the trailing p samples will be used.
% x     ... x(t) = phi(1)*x(t-1)+phi(2)*x(t-2)+...+phi(p)*x(t-p) +
%                  + z(t)+theta(1)*z(t-1)+theta(2)*z(t-2)+...+theta(q)*z(t-q),
%                  where z is gaussian zero-mean white noise with
%                  variance sigma^2.
%                  t=1,2,...,n; x(1-j)=x0(n0+1-j) for j=1,2,...,p.
%           column vector of length n.
% phir  ... roots of polynomial 1-phi(1)*z-...-phi(p)*z^p
% thetar... roots of polynomial 1+theta(1)*z+...+theta(q)*z^q

% (C) 1996-97, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(0,5,nargin));
if nargin<5,      x0   = []; end
if nargin<4,      sigma= []; end
if nargin<3,      n    = []; end
if nargin<2,      theta= []; end
if nargin<1,      phi  = []; end
%
if isempty(sigma), sigma=1; end
if isempty(n), n=100; end
p   = length(phi);
q   = length(theta);
phi = reshape(phi,1,p);
phir  = roots([-fliplr(phi),1]);
if ~all(abs(phir)>1)
   disp('WARNING: Process is not causal');
end
theta  = [1,reshape(theta,1,q)];
thetar = roots([fliplr(theta)]);
if ~all(abs(thetar)>1)
   disp('WARNING: Process is not invertible');
end
n0  = length(x0);
x0  = reshape(x0,1,n0);
x0  = [zeros(1,p-n0),x0(max(1,n0-p+1):n0)];
x0  = fliplr(x0);
if p==0
   x = sigma*randn(n+q,1);
   if q>0
      x = conv(x,theta);
      x = x(q+1:q+n);
%      x = x(1:n);
   end
   return
end
% z   = sigma*randn(1,q+1);
z   = [sigma*randn(1),zeros(1,q)];
x   = zeros(n,1);
for k=1:n
    x(k) = dot(phi,x0) + dot(theta,z);
    x0   = [x(k),x0(1:p-1)];
    z    = [sigma*randn(1),z(1:q)];
end
