function xM = Var1SimPs(A,epsilon,Tbig,x0);
%Var1SimPs   Calculates impulse response function of a VAR(1) system.
%
%             x(t) = A * x(t-1) +  epsilon(t), where x(t) is nx1
%
%  Usage:     xM = Var1SimPs(A,epsilon,Tbig,x0) or
%                = Var1SimPs(A,epsilon,Tbig)
%
%  Input:     A             nxn VAR(1) matrix, see above
%             epsilon       nx1 or 1xn vector of shocks in inital period, or Tbigxn
%                           matrix with shocks in all periods
%             Tbig          scalar, last period to calculate for
%             x0            nx1 or 1xn vector with starting values, optional
%
%  Output:    xM            Tbig x n matrix, impulse response function
%
%
%
%  Paul.Soderlind@unisg.ch, May 2001
%-----------------------------------------------------------------------

n = size(A,1);

if (nargin == 3);                %default x0
  x0 = zeros(n,1);
end;


if size(epsilon,1)*size(epsilon,2) == n;
  epsilon2       = zeros(Tbig,n);                   %non-zero only in first period
  epsilon2(1,: ) = epsilon(:)';                     %to accomodate Tbig=1
  epsilon        = epsilon2;
end;


x1_t_1 = x0(:);                                %starting vector
xM     = zeros(Tbig,n);                        %to put results in
xM(:)  = NaN;
for t = 1:Tbig;                                %loop over time periods
  x1      = A*x1_t_1 + epsilon(t,:)';
  xM(t,:) = x1';
  x1_t_1  = x1;
end;
%-----------------------------------------------------------------------