function [ Y ] = Simula( U, E, theta, xo )
%SIMULA simulate the output of the following ARMAX model
% y(k+1) = a*y(k) + b*u(k) + c*e(k) + e(k+1)
%
%   [ Y ] = Simula( U, E, theta, xo )
%
%Input values
%   U       : column vector storing the input data
%   E       : column vector storing the error realization
%   theta   : column vector storing the paramethers [a b c]'
%   xo      : initial state [y(0) U(0) e(0)]

N = length(U);
Y = [xo(1); NaN*ones(N,1)]; 
U = [xo(2); U];
E = [xo(3); E];

for i = 1:N
    Y(i+1) = [Y(i) U(i) E(i)]*theta + E(i+1);
end
Y(1) =[];