function [history] = mymkirf2_stoch(dr_,ys_,nperiods,shock,init)

%  dr_ structure returned by dynare (AR form for the model)
%  ys_ steady state vector
%  nperiods number of periods for IRFs
%  shock  vector to be used to generated irfs
%  init vector with initial positions for the endogenous variables

order =2;

if ~exist('init','var')
    init = 0*ys_;
end

if ~exist('burn_in','var')
    burn_in = 1000;
end


     
history = getpath(dr_,ys_,nperiods,shock,order,init);
history = history(:,2:end);
    
         


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBROUTINES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function history = getpath(dr_,ys_,nperiods,shock,order,h0)

nvars = size(dr_.ghx,1);
nshocks = size(dr_.ghu,2);
statevar_pos = (dr_.nstatic +1):(nvars-dr_.nfwrd);
%statevar_pos = (dr_.nstatic +1):(dr_.nstatic+size(dr_.ghx,2));

% if (max(size(shock)) > nshocks) | (max(size(shock))<1) 
%     error('erroneous shock vector as argument')
% end
% 
% if ( size(shock,1)<size(shock,2) )
%     shock = shock';
% end


history = zeros(nvars,nperiods+1);
shock_length = size(shock,1);

for i = 2:nperiods+1 
    if i-1 <= shock_length
        history(:,i) = dr_.ghx*history(statevar_pos,i-1) + dr_.ghu*shock(i-1,:)';
    else 
        history(:,i) = dr_.ghx*history(statevar_pos,i-1);
    end;
end
   


if order>1
     history2 = zeros(nvars,nperiods+1);
     if nargin > 5
         history2(:,1) = h0;
     end
     for i = 2:nperiods+1
         if i-1 <= shock_length
           history2(:,i) = 0.5*dr_.ghs2 + dr_.ghx*history2(statevar_pos,i-1) + dr_.ghu*shock(i-1,:)' + ... 
                     0.5*dr_.ghxx*kron(history(statevar_pos,i-1),history(statevar_pos,i-1)) + ...
                     0.5*dr_.ghuu*kron(shock(i-1,:)',shock(i-1,:)') + ...
                     dr_.ghxu*kron(history(statevar_pos,i-1),shock(i-1,:)');
         else
         history2(:,i) = 0.5*dr_.ghs2 + dr_.ghx*history2(statevar_pos,i-1) + ... 
                         0.5*dr_.ghxx*kron(history(statevar_pos,i-1),history(statevar_pos,i-1));
         end;
     end
     history = history2;
 end


