function X = arf(X,A,Z,RANGE,varargin)
% arf  Run autoregressive function on a tseries object.
%
% Syntax
% =======
%
%     X = arf(X,A,Z,RANGE,...)
%
% Input arguments
% ================
%
% * `X` [ tseries ] - Input data from which initial condition will be
% taken.
%
% * `A` [ numeric ] - Vector of coefficients of the autoregressive
% polynomial.
%
% * `Z` [ numeric | tseries ] - Exogenous input or constantn in the
% autoregressive process.
%
% * `RANGE` [ numeric | Inf ] - Date range on which the new time series
% observations will be computed; `RANGE` does not include pre-sample
% initial condition. `Inf` means the entire possible range will be used
% (taking into account the length of pre-sample initial condition needed).
%
% Output arguments
% =================
%
% * `X` [ tseries ] - Output data with new observations created by running
% an autoregressive process described by `A` and `Z`.
%
% Description
% ============
%
% The autoregressive process has one of the following forms:
%
%     a1*x + a2*x(-1) + ... + an*x(-n) = z,
%
% or
%
%     a1*x + a2*x(+1) + ... + an*x(+n) = z,
%
% depending on whether the range is increasing (running forward in time),
% or decreasing (running backward in time). The coefficients `a1`,...`an`
% are gathered in the `A` vector,
%
%     A = [a1,a2,...,an].
%
% Example
% ========
%
% The following two lines create an autoregressive process constructed from
% normally distributed residuals,
%
% $$ x_t = \rho x_{t-1} + \epsilon_t $$
%
%     rho = 0.8;
%     X = tseries(1:20,@randn);
%     X = arf(X,[1,-rho],X,2:20);
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

if nargin < 4
    RANGE = Inf;
end

% Parse input arguments.
P = inputParser();
P.addRequired('X',@istseries);
P.addRequired('A',@isnumeric);
P.addRequired('Z',@(x) isnumericscalar(x) || istseries(x));
P.addRequired('RANGE',@isnumeric);
P.parse(X,A,Z,RANGE);

%**************************************************************************

A = A(:).';
order = length(A) - 1;

% Work out range (includes pre/post-sample initial condition).
xfirst = X.start;
xlast = X.start + size(X.data,1) - 1;
if isequal(RANGE,Inf)
    RANGE = xfirst : xlast;
end
if RANGE(1) <= RANGE(end)
    time = 'forward';
    RANGE = RANGE(1)-order : RANGE(end);
else
    time = 'backward';
    RANGE = RANGE(end) : RANGE(1)+order;
end

% Get endogenous data.
xdata = rangedata(X,RANGE);
xsize = size(xdata);
xdata = xdata(:,:);

% Do noting if effective range is empty.
nper = length(RANGE);
if nper <= order
    return
end

% Get exogenous (z) data.
if istseries(Z)
    zdata = mygetdata(Z,RANGE);
    zdata = zdata(:,:);
    % expand zdata in 2nd dimension if needed
else
    if isempty(Z)
        Z = 0;
    end
    zdata = Z(ones([1,nper]),:);
end
if size(zdata,2) == 1 && size(xdata,2) > 1
    zdata = zdata(:,ones([1,size(xdata,2)]));
end

% Normalise polynomial vector.
if A(1) ~= 1
    zdata = zdata / A(1);
    A = A / A(1);
end

% Run AR.
if strcmp(time,'forward')
    shifts = -1 : -1 : -order;
    timevec = 1+order : nper;
else
    shifts = 1 : order;
    timevec = nper-order : -1 : 1;
end

for i = 1 : size(xdata,2)
    for t = timevec
        xdata(t,i) = -A(2:end)*xdata(t+shifts,i) + zdata(t,i);
    end
end

% Update output series.
X = subsasgn(X,RANGE,reshape(xdata,[size(xdata,1),xsize(2:end)]));

end
