% generate shock decomposition and correlation of sample innovations from Dynare output
% 
% Giovanni.Lombardo@ecb.int (Jan 2009)


close all;
pause(.1);
lgy_=M_.endo_names;
dr_=oo_.dr;

% Reconstruct VAR representation
A_mx = [zeros(size(lgy_,1),dr_.nstatic),dr_.ghx,zeros(size(lgy_,1),size(lgy_,1)-dr_.npred-dr_.nstatic)];
B_mx = dr_.ghu;

nobs=options_.nobs;
% State transition matrices
% S_mx1=dr_.ghx(dr_.nstatic+1:dr_.nstatic+size(dr_.ghx,2),:);
% S_mx2=dr_.ghu(dr_.nstatic+1:dr_.nstatic+size(dr_.ghx,2),:);

%%%%%%%%%%%% get the shocks
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

variablelist=fields(oo_.SmoothedVariables);
szsh=size(eval(['oo_.SmoothedVariables.',variablelist{1}]),1);
vary_t = [];

for jj=1:size(variablelist,1);
vary_t(:,jj) = eval(['oo_.SmoothedVariables.',variablelist{jj}]);
eval([variablelist{jj},'=',['oo_.SmoothedVariables.',variablelist{jj},';']]);
end
options_.presample=max(1,options_.presample);
sh_t0=sh_t(options_.presample:end,:);
vary_t0=vary_t(options_.presample:end,:);
stst=oo_.steady_state(dr_.order_var);


listofvariables=variablelist;
% for jj=1:size(listofvariables,1)
%     eval(['hist_',deblank(listofvariables{jj}),'=',...
%         'zeros(szsh,size(dr_.ghu,2));']);
% end
%%%%%%%%%% generate series  using smoothed shocks
yfo=(vary_t0(1,:)');
yfo_sub(:,1,size(sh_t0,2))=yfo*0;
yfo_exp=(vary_t0(1,:)');
yfo_sub_exp=yfo_sub;
% yfo_sub_exp_tot=[];


    for jj=2:size(sh_t0,1)

        yfo(:,jj)=A_mx*yfo(:,jj-1)+B_mx*sh_t0(jj,:)';
        yfo_exp(:,jj)=A_mx*yfo(:,jj-1);
        for kk=1:size(sh_t0,2)
%             yfo_sub(:,jj,kk)=A_mx*yfo(:,jj-1)+B_mx(:,kk)*sh_t0(jj,kk)';
            yfo_sub(:,jj,kk)=A_mx*yfo_sub(:,jj-1,kk)+B_mx(:,kk)*sh_t0(jj,kk)';
%             yfo_sub(:,jj,kk)=yfo_sub(:,jj,kk)+A_mx^(jj-1)*yfo_sub(:,1,kk);
            yfo_sub_exp(:,jj,kk)=yfo_sub(:,jj,kk)-yfo_exp(:,jj);
            
        end
    end
    
    
    yfo_sub_exp_tot=yfo_sub_exp+repmat(yfo_exp,[1,1,size(yfo_sub,3)])/size(yfo_sub,3);
    
    % to check that things sum up run the following: notice that there is a constant of discrepancy vs. data
    % due to non-zero steady-states for drifting variables
    % NOTICE that yfo returns the true series. The discrepancy wrt the
    % decomposed sequence is:
    %max([yfo(end,:)'-sum(yfo_sub_exp(end,:,:),3)'-yfo_exp(end,:)'])
    %     figure;plot([(yfo(end,:))',sum(yfo_sub_exp(end,:,:),3)'+yfo_exp(end,:)'])
    % legend('true','smoothed')
            
%     yfo=yfo(dr_.inv_order_var,:); %turn order into lgy_ order but here I'm following the existing order

%%%%%%%%%%%%%%%%%%% plots of desired variables/shocks
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');
tempif=0;
if size(listofvariables,1)>size(variablelist,1)
    disp('Too many variables declared in file');
end
for kk=1:size(listofvariables,1);
   eval(['pos_',deblank(listofvariables{kk}),'=[',...
       num2str(find_incell(variablelist,['\<',deblank(listofvariables{kk}),'\>'])),'];']);
end




while tempif==0
selectsh=input('Give name of file with list of selected shocks (with extension) ','s');
if exist(selectsh,'file');
    tempif=1;
end
end
selsh=strvcat(textread(selectsh,'%s\n','whitespace','\n','commentstyle','matlab'));
[xx,ia]=intersect(shockslist,deblank(selsh));
shockslist_sh=shockslist(ia,:);


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)-options_.presample+1;
lastdate=find(wholesmpl==sampl(3)+(sampl(4)-1)/4)-options_.presample+1;
% datetick('x',27,'keepticks');
% hText = xticklabel_rotate([1:4:size(pii,1)]',45,dates(1:4:end,:),'FontSize',14); 

% The following plots show the contribution of each shock through history (i.e. as if only that shock hit the 
% economy. Also it shows the sum of the contributions. The latter is not the same as the true data. To obtain this one
% should subtract the mean process (conditional expected value) of the economy under one shock at the time and add
% the mean process under all shocks together.
% With reference to previously computed variables, yfo_sub_exp_tot is the variable that would sum up to the observed
% data.
tmp=[];
if ~exist('figures/shockdec','dir')
    mkdir('figures/shockdec');
end
shks=regexprep(regexprep(cell2str(regexprep(shockslist_sh,'\_','\-'),0),'\[|\]','"'),',','", "');
for kk=1:size(listofvariables,1)
    figure(kk)
%     for jj=1:size(shockslist_sh,1);
%         bar(eval(['hist_',listofvariables{kk},'(:,',num2str(ia(jj)),')']));
%         hold on
%     end
%     hold off

%   area(xticks,squeeze(yfo_sub(eval(['pos_',deblank(listofvariables{kk})]),firstdate:lastdate,ia)))
tmp_yfo=yfo_sub(eval(['pos_',deblank(listofvariables{kk})]),firstdate:lastdate,ia);
% bar(xticks,squeeze(yfo_sub(eval(['pos_',deblank(listofvariables{kk})]),firstdate:lastdate,ia)),'stacked')
% figure;
nrw=round(sqrt(length(ia)+1));
cls=ceil(((length(ia)+1)/nrw));
maxy=max(max(abs(tmp_yfo)));%./max(abs(sum(tmp_yfo,3)));
for count_ia=1:length(ia)
    subplot(cls,nrw,count_ia);
    colr=count_ia/length(ia);
    if sum(tmp_yfo(:,:,count_ia))<0;
        colr='r';
    else
        colr='b';
    end
h=ciplot(zeros(size(tmp_yfo,2),1),tmp_yfo(:,:,count_ia),xticks,colr);
% set(gca,'YLim', [-maxy maxy])
datetick('x',27,'keepticks','keeplimits');
title(regexprep(shockslist_sh{(count_ia)},'\_','-'));

end
subplot(cls,nrw,count_ia+1);
h=plot(xticks,sum(tmp_yfo,3));

datetick('x',27,'keepticks','keeplimits');
title('Sum of shocks');


[a3,h3]=suplabel(['Shock: ',regexprep([listofvariables{kk}],'\_','-')],'t');

% csvwrite(['figures/shockdec/',deblank(listofvariables{kk}),'.csv'],squeeze(yfo_sub(eval(['pos_',deblank(listofvariables{kk})]),firstdate:lastdate,ia)))
dlmwrite(['figures/shockdec/',deblank(listofvariables{kk}),'.csv'],[shks],'delimiter',' ');
dlmwrite(['figures/shockdec/',deblank(listofvariables{kk}),'.csv'],squeeze(yfo_sub(eval(['pos_',deblank(listofvariables{kk})]),firstdate:lastdate,ia)),'-append','precision','%10.10f','delimiter',',');
dlmwrite(['figures/shockdec/',deblank(listofvariables{kk}),'.csv'],[repmat('"',length(xticks),1),datestr(xticks,'yyQQ'),repmat('"',length(xticks),1)],'-append','delimiter','');

%     legend(shockslist_sh);
%     title([listofvariables{kk}],'FontSize',20)
%     hold on
%     plot(xticks,sum(yfo_sub(eval(['pos_',deblank(listofvariables{kk})]),firstdate:lastdate,ia),3),'--o','LineWidth',3);%[',num2str(ia'),']
% %         plot(xticks,(yfo(eval(['pos_',deblank(listofvariables{kk})]),:)),'--ok','LineWidth',3);%[',num2str(ia'),']
% 
%     hold off
%     datetick('x',27,'keepticks','keeplimits');
end

%%%%%%%%%%%% correlation of shocks (innovatios)
sh_t0=sh_t(max(options_.presample+1,2):end,:);
shc=corr(sh_t0);
[xx,yy]=find_incell(shockslist,'measurement');
uu=setdiff([1:size(shockslist,1)]',yy);
shockslist=shockslist(uu,:);
shc=shc(uu,uu);
nrows1=round(length(uu)/2);
nrows2=length(uu)-nrows1;
% divide the corr-matrix in three blocks


tmp1=triu(shc(1:nrows1,1:nrows1));
tmp2=triu(shc(1:nrows1,nrows1+1:end));
tmp3=triu(shc(nrows1+1:end,nrows1+1:end));
shockslist1=shockslist(1:nrows1,:);
shockslist2=shockslist(nrows1+1:end,:);
%1)
table(['Correlation of innovations'],...
    strvcat(' ',strvcat(shockslist1)),strvcat(shockslist1),tmp1,8,6, ...
	       2);
 tb=tabtex(strvcat('shocks ',strvcat(regexprep(shockslist1,'\_','\\_'))),['Correlations'],1,...
        (strvcat(regexprep(shockslist1,'\_','\\_'))),tmp1);
tabtex2texpdf(['correlations_shocks1'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');


table(['Correlation of innovations'],...
    strvcat(' ',strvcat(shockslist1)),strvcat(shockslist2),tmp2,8,6, ...
	       2);
 tb=tabtex(strvcat('shocks ',strvcat(regexprep(shockslist1,'\_','\\_'))),['Correlations'],1,...
        (strvcat(regexprep(shockslist2,'\_','\\_'))),tmp2);
tabtex2texpdf(['correlations_shocks2'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');


table(['Correlation of innovations'],...
    strvcat(' ',strvcat(shockslist2)),strvcat(shockslist2),tmp3,8,6, ...
	       2);
 tb=tabtex(strvcat('shocks ',strvcat(regexprep(shockslist2,'\_','\\_'))),['Correlations'],1,...
        (strvcat(regexprep(shockslist2,'\_','\\_'))),tmp3);
tabtex2texpdf(['correlations_shocks3'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');

% [i,j]=find(abs(shc)==max(max(abs(shc-eye(size(shc))))));

return
