% back out shocks from dynare state space solution
% giannilmbd@gmail.com
% USE:
% 1) run dynare on your model of interest
% 2) Make sure you save in a file:
%         a) the list of conditioning variables
%         b) the list of conditioning shocks
%         c) the path of the conditioning variables (with periods in colums)
% 3) run this program
% OUTPUT: conditional IRFS are in
% oo_.cond_path.<name_of_variable_name_of_innovation>

dr_=oo_.dr;
lgy_=M_.endo_names;
nstat=M_.nstatic;
npred=M_.nspred;

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

% % [i,j]=sort(dr_.order_var);
% lgy2_=lgy_(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),dr_.nstatic),dr_.ghx,zeros(size(lgy_,1),size(lgy_,1)-dr_.npred-dr_.nstatic)];
% % 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,:)] ;
[A_var,Sig_var]=dynare_VAR_representation;
lgy2_=lgy_;
% load name of variables to condition on
disp('=================================================================================')
name_file_var=input('Give name of file with list of variables for conditioning = ','s');
if ~exist(name_file_var,'file');
    error('File does not exist')

end
disp('=================================================================================')
name_file_path=input('Give name of file with paths of variables for conditioning = ','s');
if ~exist(name_file_path,'file');
    error('File does not exist')
    
end
disp('=================================================================================')
name_file_shocks=input('Give name of file with list of shocks for conditioning = ','s');
if ~exist(name_file_shocks,'file');
    error('File does not exist')
    
end
disp('=================================================================================')
% retrieve information on path,conditioning variable list and list of
% shocks
matrix_paths=dlmread(name_file_path);
fid=fopen(name_file_var,'r');
varlist=textscan(fid,'%s');
varlist=varlist{1};
fclose(fid);
fid=fopen(name_file_shocks,'r');
shocklist=textscan(fid,'%s');
shocklist=shocklist{1};
fclose(fid);
%%% NOTE NEED TO CHECK ORDERING OF BELOW
pos_var_vec=[];
for jj=1:size(varlist,1);
    pos_var_vec(jj,1)=find_incell(lgy2_,varlist{jj});
end

pos_shock_vec=[];
for jj=1:size(shocklist,1);
    pos_shock_vec(jj,1)=find_incell(M_.exo_names,shocklist{jj});
end

% invertibility condition
M=Sig_var(pos_var_vec,pos_shock_vec);

if rcond(M)<1e-16;
    error('Invertibility condition failed');
end
% extract shocks
shock_path=zeros(size(pos_shock_vec,1),size(matrix_paths,2));
Ma_expr=zeros(size(pos_shock_vec,1),1);
MXA=[];
for jj=1:size(matrix_paths,2);
    % solve for shocks
    shock_path(:,jj)=M\(matrix_paths(:,jj)-Ma_expr);
    % create [A^jj*C+A^(jj-1)*C+...]
    MXA=[A_var^(jj)*Sig_var(:,pos_shock_vec),MXA];
    % multiply this by the vectorized path of shocks
    Vmx=MXA*vec(shock_path(:,1:jj));
% extract conditioning variables.
    Ma_expr=Vmx(pos_var_vec,:);
end
Sigma_e=M_.Sigma_e;
% check results using VAR
C=eye(size(Sigma_e,1));
C=C(:,pos_shock_vec);
Y=zeros(size(A_var,1),size(shock_path,2));
for jj=2:size(shock_path,2)+1
    Y(:,jj)=A_var*Y(:,jj-1)+Sig_var*C*shock_path(:,jj-1);
end
% extract path of conditioning variables for testing
CondY=Y(pos_var_vec,:);
% test
test=CondY(:,2:end)-matrix_paths;
if max(max(abs(test)))>1e-11;
    error('discrepancy in the solution')
end
% extend path beyond conditioning
for jj=size(shock_path,2)+2:size(shock_path,2)+12
    Y(:,jj)=A_var*Y(:,jj-1);
end
% reorder results
YY=Y;

% Assign response to variable names to structure oo_.cond_path
oo_.cond_path=[];
for jj=1:size(shocklist,1)
    for kk=1:size(M_.endo_names,1);
        eval([regexprep(['oo_.cond_path.',M_.endo_names(kk,:),'_',shocklist{jj}],'\s*',''),'=[',num2str(YY(kk,:),8),'];']);
    end
end

% example of plotting results: pick the variables you want!
figure;plot(oo_.cond_path.TRc_e_trc)
% figure;
% plot(oo_.cond_path.BLT_OBS_US_RES_BLT_NEWS_US,'k-o')
% hold on
% plot(oo_.cond_path.BLT_OBS_EA_RES_BLT_NEWS_US,'b-x')
% legend('US','EA')
% title('BLT')
% figure;
% plot(oo_.cond_path.RP_HAT_US_RES_BLT_NEWS_US,'k-o')
% hold on
% plot(oo_.cond_path.RP_HAT_EA_RES_BLT_NEWS_US,'b-x')
% legend('US','EA')
% title('Risk-Premium')
% figure;
% plot(oo_.cond_path.Y_HAT_US_RES_BLT_NEWS_US,'k-o')
% hold on
% plot(oo_.cond_path.Y_HAT_EA_RES_BLT_NEWS_US,'b-x')
% legend('US','EA')
% title('GDP')
% figure;
% plot(oo_.cond_path.PI_HAT_US_RES_BLT_NEWS_US,'k-o')
% hold on
% plot(oo_.cond_path.PI_HAT_EA_RES_BLT_NEWS_US,'b-x')
% legend('US','EA')
% title('Inflation')