function z = arimasim(phi,theta,d,zinit,n)
% ARIMASIM(PHI,THETA,D,ZINIT,N) generates N observations from the
%   ARIMA model specified by PHI, THETA, D, ZINIT.
% It "burns in" with the 1st 200 observations
% ZINIT must be a d-vector with entries
%   (z(-d+1),Dz(-d+2),...,D^{d-1}z(0)).
[p m] = size(phi);
[q m] = size(theta);
n1 = 200+n;
a = normrnd(0,1,n1+q,1);
z = zeros(p,1);
for i=1:n1
   zt = z(i:(i+p-1))'*phi(p:-1:1)+a(i+q)-a(i:(i+q-1))'*theta(q:-1:1);
   z = [ z ; zt ];
end
z = z((201+p):(n1+p));
if d>0
   for k=1:d
      z = cumsum( [ zinit(d-k+1) ; z ]);
   end
   z = z(1:n);
end
