% based on sw_histdec
% Oreste Tristani
% 09/08/05
% Modified by Gianni Lombardo (latest 12-08)
% and Gianni (6-April-2012)
% Produces historical decomposition of output, inflation and the nominal
% interest rate after estimation.
% By modifying the path of the shocks in oo_.SmoothedShocks, get an 
% alternative scenario
close all
% fname=input('give names of model with no extension :','s');
% if exist([fname,'_results.mat'],'file')
%     eval(['load ',[fname,'_results.mat']]);
% else
%     disp('no results file was found');
% end


dr_=oo_.dr;
lgy_=M_.endo_names;
start_forward=10;%find(ttt==1998.75);% HERE MUST ENTER THE DATE FROM WHICH TO START THE FORECAST
forecyes='noforc'; %give a name to save with the file (e.g. if include forecast or not);
nstat=dr_.nstatic;
npred=dr_.nspred;
if exist([options_.datafile,'.mat'])
eval(['load ',options_.datafile])
else
    eval(options_.datafile)
end
nobs=(options_.nobs);
ttt=[1:nobs];%1985:0.25:2009.75;
varlist=input('Give name of file with selected variables ','s');
listofvariables=textread(varlist,'%s\n','whitespace','\n','commentstyle','matlab');%strvcat('doutput','inflobs','robs','rr_f_8q');

% [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)];

% 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,:)] ;

% Load the shocks for the historical decomposition
shockslist=fields(oo_.SmoothedShocks);
szsh=size(eval(['oo_.SmoothedShocks.',shockslist{1}]),1);
sh_t = zeros(szsh,size(dr_.ghu,2)) ;

for jj=1:size(shockslist,1);
sh_t(:,jj) = eval(['oo_.SmoothedShocks.',shockslist{jj}]);
end




% Forecastable components necessary to compute variables
inits_tt=[];
nmvar=fields(oo_.SmoothedVariables);
tmpp=struct2cell(oo_.SmoothedVariables);
for jj=1:size(tmpp,1);
    inits_tt=[inits_tt;tmpp{jj}(1)];
end

% Forecastable components
% inits_tt_for=[];
% tmpp_for=struct2cell(oo_.SmoothedVariables);
% for jj=1:size(tmpp,1);
%     inits_tt_for=[inits_tt_for;tmpp_for{jj}(start_forward)];
% end


inits_t=inits_tt([nstat+[1:npred],c']);
% inits_t_for=inits_tt_for([nstat+[1:npred],c']);
forec = zeros(size(A_var,2),szsh);
% initialize variables with first period from Smoothed variables
forec(:,1) = inits_t ;

% forec_for= zeros(size(A_var,2),szsh-start_forward);
% forec_for(:,1) = inits_t_for;
for i = 2:szsh
    forec(:,i) = A_var*forec(:,i-1) ;
    
end;

% VMA representation  HERE COULD SIMPLY RECONSTRUCT SERIES
Phi_var = zeros(size( A_var,2),size(dr_.ghu,2),szsh);
Phi_var(:,:,1) = Sig_var;
y_reconst=zeros(size(Phi_var,1),szsh);
y_reconst(:,1)=inits_t;
Phi_var(:,:,2) = A_var*Phi_var(:,:,1) ;
for i = 2:szsh
    Phi_var(:,:,i+1) = A_var*Phi_var(:,:,i) ;
    y_reconst(:,i)=A_var*y_reconst(:,i-1)+Sig_var*sh_t(i,:)';
end;
% i=i+1;
% Phi_var(:,:,i+1) = A_var*Phi_var(:,:,i) ;


%  insitialiye hist_
for jj=1:size(listofvariables,1)
    eval(['hist_',deblank(listofvariables{jj}),'=',...
        'zeros(szsh,size(dr_.ghu,2));']);
end
        
    
% for shocknum = 1:size(dr_.ghu,2)
%     sh_t_num = zeros(size(sh_t)) ;
%     sh_t_num(:,shocknum) = sh_t(:,shocknum) ;
%     for peri = 2:szsh
%         whatever = zeros(size(A_var,2),1) ;
%         for i = 2:peri % use MA to reconstruct variables out of smoothed shocks
%             whatever = whatever + squeeze(Phi_var(:,:,i-1))*sh_t_num(peri+1-i,:)' ;
%         end; 
%         for jj=1:size(listofvariables,1)
%           eval(['hist_',deblank(listofvariables{jj}),'(',num2str(peri),',',num2str(shocknum),')=',...
%         'whatever(end-size(listofvariables,1)+',num2str(jj),') ;']);
%         end
% 
%     end;
% end;

% get info about data for plotting
sampl=input('give vector with beginning and end of subsample (e.g.[year_start,quarter_start,year_end,q_end]); ');
stdate=input('give vector of year and quarter of first observation in whole sample');
if isempty(stdate)
    stdate=[1980,1];
end
stdate0=stdate(1)+(stdate(2)-1)/4;
wholesmpl=stdate0:.25:stdate0+(nobs-1)/4;
if isempty(sampl)
    sampl=[stdate(1),stdate(2),round(wholesmpl(end)),1];
end
[xticks]=setdates2(sampl(1),sampl(2),sampl(3),sampl(4)); % FIX THIS TOO
firstdate=find(wholesmpl==sampl(1)+(sampl(2)-1)/4);
lastdate=find(wholesmpl==sampl(3)+(sampl(4)-1)/4);
% reconstruct series out of each individual shock
tot_var=[];
    for jj=1:size(listofvariables,1)
        tot_var(:,jj)=sum(eval(['hist_',listofvariables{jj}])');
    end
    tot_var=tot_var';
    szsh=size(eval(['oo_.SmoothedVariables.',listofvariables{1}]),1);% size of the sample
    for kk=1:size(listofvariables,1)
        figure(200+kk);
%         plot(xticks(1:end-1),eval(['oo_.SmoothedVariables.',listofvariables{kk},'(firstdate:lastdate-1)']),'o',xticks(1:end-1),tot_var(kk,firstdate+1:lastdate)+forec(eval([listofvariables{kk},'_pos']),firstdate+1:lastdate),'-x') ;
        h=plot(xticks(1:end),eval(['oo_.SmoothedVariables.',listofvariables{kk},'(firstdate:lastdate)']),':ro',xticks(1:end),y_reconst(eval([listofvariables{kk},'_pos']),firstdate:lastdate),'-kx') ;
        set(h(1),'LineWidth',[2])
        set(h(2),'LineWidth',[2])
        ylabel(regexprep(listofvariables{kk},'\_','\\_'));
        datetick('x',27);%'keepticks',,'keeplimits'
        xticklabel_rotate
    end
%     legend('Smoothed','Counterfactual')
    cols=round(sqrt(kk));
    rws=round((kk+1)/cols);
    graph_merge_subplots(201:200+kk,[rws,cols])


