function [xp,z,v] = armasim(x,phi,theta,x0,z0)
%           VALIDATE AN ARMA(p,q) PROCESS VIA ONE-STEP SIMULATION
% function [xp,z,v] = armasim(x,phi,theta,x0,z0)
% x     ... observed data set which is to be validated as ARMA(p,q)
%           process with coefficients phi and theta, n=length(x).
% phi   ... AR coefficients (estimated), p=length(phi),   default=[] <=> p=0
% theta ... MA coefficients (estimated), q=length(theta), default=[] <=> q=0
%           p=0 => x=MA(q) process, q=0 => x=AR(p) process
%           p=0 & q=0 => x=zero-mean white noise process.
% x0    ... initial data sequence relevant to ARMA(p,q), p>0, only;
%           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.
% z0    ... initial data sequence relevant to ARMA(p,q) q>0, only;
%           default=[], m0=length(z0).
%           n0<q => the leading part of z0 updated with q-n0 zeros
%           n0>q => only the trailing q values will be used.
% xp    ... xp(t) = phi(1)*x(t-1)+phi(2)*x(t-2)+...+phi(p)*x(t-p) +
%                   + theta(1)*z(t-1)+theta(2)*z(t-2)+...+theta(q)*z(t-q),
%           t=1,...,n, is one-step prediction (column vector),
%           where x(-k) = x0(n0-k) for k=0,...,p-1.
% z     ... = x-xp ... prediction (residual) error which should be
%           zero-mean white noise process, length(z)=n, column vector.
%           Do verify that actually z is white noise using trand(z)
%           and/or inspecting the behaviour of acf(z).
%           Provided that this cannot be confirmed, the model has to be
%           rejected as inappropriate.
% v     ... v=norm(z).^2/n is estimated variance of z.

% (C) 1997-99, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(1,5,nargin));
if nargin<5, z0 = []; end
if nargin<4, x0 = []; end
if nargin<3, theta = []; end
if nargin<2, phi   = []; end
%
p = length(phi);
q = length(theta);
phi = [0,reshape(phi,1,p)];
n  = length(x);
n0  = length(x0);
x0  = reshape(x0,n0,1);
x0  = [zeros(p-n0,1);x0(max(1,n0-p+1):n0)];
n0  = length(x0);
x   = [x0;reshape(x,n,1)];

if p>0
   xp=conv(x,phi); xp=xp(n0+1:n0+n);
else
   xp = zeros(n,1);
end

e = x(n0+1:n0+n) - xp;

if q>0
   m0  = length(z0);
   z0  = reshape(z0,m0,1);
   z0  = flipud([zeros(q-m0,1);z0(max(1,m0-q+1):m0)]);
   z  = zeros(n,1);
   for k=1:n
       z(k) = e(k) - dot(theta,z0);
       z0   = [z(k);z0(1:q-1)];
   end
   xp = xp + (e - z);
else
   z = e;
end

v=norm(z).^2./n;
