function [ymean,yinit] = mean(w,alt)
% mean  Mean of VAR process.
%
% Syntax
% =======
%
%     x = mean(v)
%
% Input arguments
% ================
%
% * `v` [ VAR ] - VAR object.
%
% Output arguments
% =================
%
% * `x` [ numeric ] - Mean vector.
%
% Description
% ============
%
% Example
% ========
%

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

if ~exist('alt','var')
    alt = 1 : size(w.A,3);
end

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

[ny,p,nalt] = sizeof(w);
if p == 0
    ymean = w.K(:,alt);
    yinit = zeros([ny,0,nalt]);
    return
end

nalt = numel(alt);
realsmall = getrealsmall();
yinit = nan([ny,p,nalt]);
ymean = nan([ny,nalt]);
for i = 1 : nalt
    alti = alt(i);
    [T,R,k,Z,H,d,U] = sspace(w,alti);
    unit = abs(abs(w.eigval(1,:,alti)) - 1) <= realsmall;
    nunit = sum(unit);
    if nunit == 0
        % Stationary parameterisation.
        ymean(:,i) = sum(poly.var2poly(w.A(:,:,alti)),3) ...
            \ w.K(:,alti);
        if nargout > 1
            yinit(:,1:p,i) = ymean(:,i*ones([1,p]));
        end
    else
        % Unit-root parameterisation.
        a2 = (eye(ny*p-nunit) - T(nunit+1:end,nunit+1:end)) ...
            \ k(nunit+1:end,1);
        % Return NaNs for unit-root variables.
        dy = any(abs(U(1:ny,unit)) > realsmall,2).';
        ymean(~dy,i) = U(~dy,nunit+1:end)*a2;
        if nargout > 1
            temp = U*[zeros([nunit,1]);a2];
            temp = reshape(temp,[ny,p]);
            yinit(:,:,i) = temp(:,end:-1:1);
        end
    end
end

end
