function [s,c,s1,c1] = srf(THIS,TIME,varargin)
% srf  Shock (impulse) response function.
%
% Syntax
% =======
%
%     [R,CUM] = srf(V,NPER)
%     [R,CUM] = srf(V,RANGE)
%
% Input arguments
% ================
%
% * `V` [ SVAR ] - SVAR object for which the impulse response function will
% be computed.
%
% * `NPER` [ numeric ] - Number of periods.
%
% * `RANGE` [ numeric ] - Date range.
%
% Output arguments
% =================
%
% * `R` [ tseries | struct ] - Shock response functions.
%
% * `CUM` [ tseries | struct ] - Cumulative shock response functions.
%
% Options
% ========
%
% * `'presample='` [ `true` | *`false`* ] - Include zeros for pre-sample
% initial conditions in the output data.
%
% * `'select='` [ cellstr | char | logical | numeric | *`Inf`* ] - Selection
% of shocks to which the responses will be simulated.
%
% Description
% ============
%
% For backward compatibility, the following calls into the `srf` function is
% also possible:
%
%     [~,~,s,c] = srf(this,nper)
%     [~,~,s,c] = srf(this,range)
%
% Example
% ========
%

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

opt = passvalopt('SVAR.srf',varargin{:});

% Tell if `time` is `nper` or `range`.
if length(TIME) == 1 && round(TIME) == TIME && TIME > 0
    range = 1 : TIME;
else
    range = TIME(1) : TIME(end);
end
nper = length(range);

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

[ny,p,nalt] = sizeof(THIS);

select = utils.nameselection(ny,THIS.enames,opt.select);
ne = sum(select);

% Compute VMA matrices.
A = THIS.A;
B = THIS.B;
Phi = timedom.var2vma(A,B,nper,select);

% Create shock paths.
Eps = zeros(ny,ne,nper,nalt);
for ialt = 1 : nalt
    E = eye(ny);
    E = E(:,select);
    Eps(:,:,1,ialt) = E;
end

% Permute dimensions so that time runs along the 2nd dimension.
Phi = permute(Phi,[1,3,2,4]);
Eps = permute(Eps,[1,3,2,4]);

% Add a total of `p` zero initial conditions.
if opt.presample
    Phi = [zeros(ny,p,ne,nalt),Phi];
    Eps = [nan(ny,p,ne,nalt),Eps];
    xrange = range(1)-p : range(end);
else
    xrange = range;
end

s = VAR.outputdata(THIS,'auto',xrange, ...
    [Phi;Eps],[],[THIS.ynames,THIS.enames]);
% For bkw compatibility.
if nargout > 2
    s1 = VAR.outputdata(THIS,'tseries',xrange, ...
        [Phi;Eps],[],[THIS.ynames,THIS.enames]);
end

if nargout > 1
    Psi = cumsum(Phi,2);
    c = VAR.outputdata(THIS,'auto',xrange, ...
        [Psi;Eps],[],[THIS.ynames,THIS.enames]);
    % For bkw compatibility.
    if nargout > 3
        c1 = VAR.outputdata(THIS,'tseries',xrange, ...
            [Psi;Eps],[],[THIS.ynames,THIS.enames]);
    end
end

end