% IRFs for Dynare
% variables in y=A_var*y(-1)+ghu*epsilon; are in the order lgy2_; then they
% are reordered as in original lgy_, when IRFs are assigned to oo_


dr_=oo_.dr;
lgy_=M_.endo_names;
nstat=M_.nstatic;
if isfield(M_,'nspred');% older verisions of Dynare (<4.3.3) had this 
npred=M_.nspred;
else
    npred=M_.npred+M_.nboth;
end;

listofvariables=str2cell(lgy_);
shocks=M_.exo_names;

% [i,j]=sort(dr_.order_var);
lgy2_=lgy_(oo_.dr.order_var,:); % variables as appear in the canonical VAR rapresentation
[a,b,c]=intersect(deblank(listofvariables),str2cell(lgy2_));
c=c(:);
listofvariables=a;
%position in forecast vector and VMA representation
for jj=1:length(c);
    eval([deblank(listofvariables{jj}),'_pos=',num2str(length(nstat+[1:npred])+jj),';']);
end

%c is the order in the output; So lgy2_(c,:) should return listofvariables

% Rewrite the law of motion as a state space representation of a VAR
% The VAR matrix includes all the states plus list of variables,
% 
% 
% 
% reordering (from Dynare manual)
% % When reordered, the variables are stored in the following order: static variables, purely predetermined
% % variables (variables that appear only at the current and lagged periods in the model), variables that are
% % both predetermined and forward-looking (variables that appear at the current, future and lagged
% % periods in the model), purely forward-looking variables (variables that appear only at the current and
% % future periods in the model). In each category, the variables are arranged alphabetically.

% The dynamics is fully described by the predetermined variables only.
% so 
% A_var = [dr_.ghx([nstat+[1:npred],c'],:),zeros(size(dr_.ghx([nstat+[1:npred],c'],:),1),length(c))];%[[dr_.ghx(16:37,:);dr_.ghx(4,:);dr_.ghx(11,:);dr_.ghx(38,:)],zeros(size(dr_.ghx,2)+3,3)];
% A_var =[zeros(size(lgy_,1),M_.nstatic),dr_.ghx,zeros(size(lgy_,1),size(lgy_,1)-M_.npred-M_.nstatic)];
A_var =[zeros(size(lgy_,1),nstat),...
    dr_.ghx,zeros(size(lgy_,1),...
    size(lgy_,1)-npred-nstat)];
% Sigma matrix for VAR representation. 
% Sig_var = dr_.ghu([nstat+[1:npred],c'],:);%[dr_.ghu(16:37,:);dr_.ghu(4,:);dr_.ghu(11,:);dr_.ghu(38,:)] ;
Sig_var = dr_.ghu;%[dr_.ghu(16:37,:);dr_.ghu(4,:);dr_.ghu(11,:);dr_.ghu(38,:)] ;



% State stpace form
% y=A_var*y(-1)+Sig_var*epsilon;
Sigma_e=M_.Sigma_e;
horizons=options_.irf;
nsh=size(Sigma_e);
z=1;
if isempty(var_list_)
       var_list_=lgy2_;
end
for jj=1:size(var_list_,1)
    eval(['pos_',var_list_(jj,:),'=find_incell(lgy2_,',regexprep(['''\<',var_list_(jj,:),'\>'''],'\s*',''),');']);
end
tmpI=eye(nsh);
for jj=1:nsh;
    tot_var=zeros(size(A_var,1),horizons);
    tot_var(:,1)=Sig_var*Sigma_e(:,jj);
    for uu=2:horizons
        tot_var(:,uu)=A_var*tot_var(:,uu-1);
    end
    
    % tot_var=tot_var(oo_.dr.inv_order_var,oo_.dr.inv_order_var);
    for kk=1:size(var_list_,1);
        eval([regexprep(['oo_.irfs.',var_list_(kk,:),'_',M_.exo_names(jj,:)],'\s*',''),'=[',num2str(tot_var(eval(['pos_',var_list_(kk,:)]),:),8),'];']);
    end
end