% Program sw_histdec
% Oreste Tristani
% 09/08/05
% Modified by Gianni Lombardo (latest 12-08)
% Produces historical decomposition of output, inflation and the nominal
% interest rate after estimation of the SW model
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
% sh_t(:,2) = oo_.SmoothedShocks.ub ;
% sh_t(:,3) = oo_.SmoothedShocks.ug ;
% sh_t(:,4) = oo_.SmoothedShocks.ui ;
% sh_t(:,5) = oo_.SmoothedShocks.up ;
% sh_t(:,6) = oo_.SmoothedShocks.ur ;
% sh_t(:,7) = oo_.SmoothedShocks.uw ;
% VMA representation
Phi_var = zeros(size( A_var,2),size(dr_.ghu,2),szsh);
Phi_var(:,:,1) = Sig_var;
for i = 1:szsh
    Phi_var(:,:,i+1) = A_var*Phi_var(:,:,i) ;
end;


% Historical decomp
% hist_doutput = zeros(szsh,size(dr_.ghu,2)) ;
% hist_robs = zeros(szsh,size(dr_.ghu,2)) ;
% hist_rr_f_8q=hist_robs;
% hist_inflobs = zeros(szsh,size(dr_.ghu,2)) ;
% 

%%%%%%%%%%%% was here

% Forecastable components
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);
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;

% for i = 2:szsh+1-start_forward; %this is the forward looking bit
%     forec_for(:,i)=A_var*forec_for(:,i-1) ;
% end
%%%%%%%%%%%%%%%%%%% was above
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
            whatever = whatever + squeeze(Phi_var(:,:,i-1))*sh_t_num(peri+1-i,:)' ;
        end; % this is not made general yet (gianni)
        for jj=1:size(listofvariables,1)
          eval(['hist_',deblank(listofvariables{jj}),'(',num2str(peri),',',num2str(shocknum),')=',...
        'whatever(end-size(listofvariables,1)+',num2str(jj),') ;']);
        end
%         hist_doutput(peri,shocknum) = whatever(size(A_var,2)-3) ;
%         hist_inflobs(peri,shocknum) = whatever(size(A_var,2)-2) ;
%         hist_robs(peri,shocknum) = whatever(size(A_var,2)-1) ;
%         hist_rr_f_8q(peri,shocknum) = whatever(size(A_var,2)) ;
    end;
end;

%%%%%%%%%%%%% plot shocks decomposition

tempif=0;
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);
lastdate=find(wholesmpl==sampl(3)+(sampl(4)-1)/4);
% datetick('x',27,'keepticks');
% hText = xticklabel_rotate([1:4:size(pii,1)]',45,dates(1:4:end,:),'FontSize',14); 
for kk=1:size(listofvariables,1)
    figure(100+kk)
%     for jj=1:size(shockslist_sh,1);
%         bar(eval(['hist_',listofvariables{kk},'(:,',num2str(ia(jj)),')']));
%         hold on
%     end
%     hold off
% area(1:nobs,eval(['hist_',listofvariables{kk},'(:,[',num2str(ia'),'])']))
area(xticks,eval(['hist_',listofvariables{kk},'(firstdate:lastdate,[',num2str(ia'),'])']))
    legend(regexprep(shockslist_sh,'\_','\\_'));
    title([regexprep(listofvariables{kk},'\_','\\_')],'FontSize',20)
    hold on
    plot(xticks,sum(eval(['hist_',listofvariables{kk},'(firstdate:lastdate,:)']),2),'--o','LineWidth',3);%[',num2str(ia'),']
    hold off
    datetick('x',27,'keepticks','keeplimits');
end












% Check that things add up
% Surprise components
if 0
    tot_var=[];
    for jj=1:size(listofvariables,1)
        tot_var(:,jj)=sum(eval(['hist_',listofvariables{jj}])');
    end
    tot_var=tot_var';
    % tot_dout = sum(hist_doutput')' ;
    % tot_robs = sum(hist_robs')' ;
    % tot_inflobs = sum(hist_inflobs')' ;
    % tot_rr_f_8q=sum(hist_rr_f_8q')';
    szsh=size(eval(['oo_.SmoothedVariables.',listofvariables{1}]),1);% size of the sample
    for kk=1:size(listofvariables,1)
        figure(200+kk);
        plot(ttt(2:szsh)',eval(['oo_.SmoothedVariables.',listofvariables{kk},'(1:end-1)']),'o',ttt(2:szsh)',tot_var(kk,2:end)+forec(eval([listofvariables{kk},'_pos']),2:end),'x') ;
    end
    % figure(11);
    % plot(ttt(1:szsh)',4*oo_.SmoothedVariables.robs,'o',ttt(1:szsh)',4*(tot_robs+forec(robs_pos,:)'),'.') ;
    % figure(12);
    % plot(ttt(1:szsh)',4*oo_.SmoothedVariables.inflobs,'o',ttt(1:szsh)',4*(tot_inflobs+forec(inflobs_pos,:)'),'.') ;
    % figure(120)
    %
    % plot(ttt(1:szsh)',4*oo_.SmoothedVariables.rr_f_8q,'o',ttt(1:szsh)',4*(tot_rr_f_8q+forec(rr_f_8q_pos,:)'),'.') ;

    % ==================================================================

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);
dyntable(['Correlation of innovations'],...
    strvcat(' ',strvcat(shockslist)),strvcat(shockslist),shc,8,6, ...
	       2);
 tb=tabtex(strvcat('shocks ',strvcat(regexprep(shockslist,'\_','\\_'))),['Correlations'],1,...
        (strvcat(regexprep(shockslist,'\_','\\_'))),shc);
tabtex2texpdf(['correlations_shocks'],1,'',tb,0);

return
hist_supply_dout = hist_doutput(:,1) + hist_doutput(:,5) + hist_doutput(:,7) ;
hist_demand_dout = hist_doutput(:,2) + hist_doutput(:,3) + hist_doutput(:,4) ;
hist_policy_dout = hist_doutput(:,6) ;

hist_supply_robs = hist_robs(:,1) + hist_robs(:,5) + hist_robs(:,7) ;
hist_demand_robs = hist_robs(:,2) + hist_robs(:,3) + hist_robs(:,4) ;
hist_policy_robs = hist_robs(:,6) ;
hist_supply_rr_f_8q = hist_rr_f_8q(:,1) + hist_rr_f_8q(:,5) + hist_rr_f_8q(:,7) ;
hist_demand_rr_f_8q = hist_rr_f_8q(:,2) + hist_rr_f_8q(:,3) + hist_rr_f_8q(:,4) ;
hist_policy_rr_f_8q = hist_rr_f_8q(:,6) ;
hist_supply_inflobs = hist_inflobs(:,1) + hist_inflobs(:,5) + hist_inflobs(:,7) ;
hist_demand_inflobs = hist_inflobs(:,2) + hist_inflobs(:,3) + hist_inflobs(:,4) ;
hist_policy_inflobs = hist_inflobs(:,6) ;

% figure(20);
% plot(ttt(1:szsh)',4*(oo_.SmoothedVariables.doutput-forec(23,:)'),'-',ttt(1:szsh)',4*hist_supply_dout,'.',ttt(1:szsh)',4*hist_demand_dout,'o',ttt(1:szsh)',4*hist_policy_dout,'-.') ;
% title('Historical decomposition of output growth');
% xlabel('time');
% legend('Actual (demeaned)','Supply','Demand','Policy',0);
% 
% figure(21);
% plot(ttt(1:szsh)',4*(oo_.SmoothedVariables.robs-forec(24,:)'),'-',ttt(1:szsh)',4*hist_supply_robs,'.',ttt(1:szsh)',4*hist_demand_robs,'o',ttt(1:szsh)',4*hist_policy_robs,'-.') ;
% title('Historical decomposition of nominal interest rate');
% xlabel('time');
% legend('Actual (demeaned)','Supply','Demand','Policy',0);
% 
% figure(22);
% plot(ttt(1:szsh)',4*(oo_.SmoothedVariables.inflobs-forec(25,:)'),'-',ttt(1:szsh)',4*hist_supply_inflobs,'.',ttt(1:szsh)',4*hist_demand_inflobs,'o',ttt(1:szsh)',4*hist_policy_inflobs,'-.') ;
% title('Historical decomposition of inflation');
% xlabel('time');
% legend('Actual (demeaned)','Supply','Demand','Policy',0);
y_hist_all_inflobs=zeros(size(hist_inflobs));
for i = 4:szsh
    y_dout(i) = oo_.SmoothedVariables.doutput(i) ;%+ oo_.SmoothedVariables.doutput(i-1) + oo_.SmoothedVariables.doutput(i-2) + oo_.SmoothedVariables.doutput(i-3) ;
    y_forec_dout(i) = forec(doutput_pos,i) ;%+ forec(23,i-1) + forec(23,i-2) + forec(23,i-3) ;
    y_hist_supply_doutput(i) = hist_supply_dout(i);% + hist_supply_dout(i-1) + hist_supply_dout(i-2) + hist_supply_dout(i-3) ;
    y_hist_demand_doutput(i) = hist_demand_dout(i);% + hist_demand_dout(i-1) + hist_demand_dout(i-2) + hist_demand_dout(i-3) ;
    y_hist_policy_doutput(i) = hist_policy_dout(i);% + hist_policy_dout(i-1) + hist_policy_dout(i-2) + hist_policy_dout(i-3) ;
    y_infl(i) = oo_.SmoothedVariables.inflobs(i) + oo_.SmoothedVariables.inflobs(i-1) + oo_.SmoothedVariables.inflobs(i-2) + oo_.SmoothedVariables.inflobs(i-3) ;
    y_forec_infl(i) = forec(inflobs_pos,i) + forec(inflobs_pos,i-1) + forec(inflobs_pos,i-2) + forec(inflobs_pos,i-3) ;
    y_hist_supply_inflobs(i) = hist_supply_inflobs(i) + hist_supply_inflobs(i-1) + hist_supply_inflobs(i-2) + hist_supply_inflobs(i-3) ;
    y_hist_demand_inflobs(i) = hist_demand_inflobs(i) + hist_demand_inflobs(i-1) + hist_demand_inflobs(i-2) + hist_demand_inflobs(i-3) ;
    y_hist_policy_inflobs(i) = hist_policy_inflobs(i) + hist_policy_inflobs(i-1) + hist_policy_inflobs(i-2) + hist_policy_inflobs(i-3) ;
    y_hist_all_inflobs(i,:)=hist_inflobs(i,:)+hist_inflobs(i-1,:)+hist_inflobs(i-2,:)+hist_inflobs(i-3,:);
end;

% figure(20)
% plot(ttt(56:szsh)',4*(oo_.SmoothedVariables.doutput(56:end)-forec(23,56:end)'),'-',ttt(56:szsh)',4*hist_supply_dout(56:end),'.',ttt(56:szsh)',4*hist_demand_dout(56:end),'o',ttt(56:szsh)',4*hist_policy_dout(56:end),'-.') ;
% title('Historical decomposition of output growth');
% xlabel('time');
% legend('Actual (demeaned)','Supply','Demand','Policy',0);
 
figure(21)
plot(ttt(56:szsh)',4*(oo_.SmoothedVariables.robs(56:end)-forec(robs_pos,56:szsh)'),'-',ttt(56:szsh)',4*hist_supply_robs(56:end),'.',ttt(56:szsh)',4*hist_demand_robs(56:end),'o',ttt(56:szsh)',4*hist_policy_robs(56:end),'-.') ;
title('Historical decomposition of nominal interest rate');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);
%%
figure(22)
plot(ttt(56:szsh-1)',diff([4*(oo_.SmoothedVariables.robs(56:end)-0*forec(robs_pos,56:szsh)')]),'k','LineWidth',3)
hold on
h=bar(ttt(56:szsh-1)',diff([4*hist_supply_robs(56:end),4*hist_demand_robs(56:end),...
    4*hist_policy_robs(56:end)]),1.5,'grouped');%,'-',ttt(56:szsh)',4*hist_supply_robs(56:end),'.',ttt(56:szsh)',4*hist_demand_robs(56:end),'o',ttt(56:szsh)',4*hist_policy_robs(56:end),'-.') ;
% bar(ttt(56:szsh)',[4*hist_demand_robs(56:end)],'g');%
% bar(ttt(56:szsh)',[4*hist_policy_robs(56:end)],'r');%
set(h(1),'FaceColor','g')
set(h(2),'FaceColor','b')
set(h(3),'FaceColor','r')

title('Historical decomposition of nominal interest rate');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);
hold off
axis tight
%%
figure(27)
plot(ttt(56:szsh-1)',diff([4*(oo_.SmoothedVariables.rr_f_8q(56:end)-0*forec(rr_f_8q_pos,56:szsh)')]),'k','LineWidth',3)
hold on
h=bar(ttt(56:szsh-1)',diff([4*hist_supply_rr_f_8q(56:end),4*hist_demand_rr_f_8q(56:end),...
    4*hist_policy_rr_f_8q(56:end)]),1.5,'grouped');%,'-',ttt(56:szsh)',4*hist_supply_robs(56:end),'.',ttt(56:szsh)',4*hist_demand_robs(56:end),'o',ttt(56:szsh)',4*hist_policy_robs(56:end),'-.') ;
% bar(ttt(56:szsh)',[4*hist_demand_robs(56:end)],'g');%
% bar(ttt(56:szsh)',[4*hist_policy_robs(56:end)],'r');%
set(h(1),'FaceColor','g')
set(h(2),'FaceColor','b')
set(h(3),'FaceColor','r')

title('Historical decomposition of natural interest rate');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);
hold off
axis tight

%%
% figure(22)
% plot(ttt(56:szsh)',4*(oo_.SmoothedVariables.inflobs(56:end)-forec(inflobs_pos,56:szsh)'),'-',ttt(56:szsh)',4*hist_supply_inflobs(56:end),'.',ttt(56:szsh)',4*hist_demand_inflobs(56:end),'o',ttt(56:szsh)',4*hist_policy_inflobs(56:end),'-.') ;
% title('Historical decomposition of inflation');
% xlabel('time');
% legend('Actual (demeaned)','Supply','Demand','Policy',0);
%%
figure(23)
plot(ttt(56:szsh)',(y_dout(56:end)-y_forec_dout(56:szsh)),'-',ttt(56:szsh)',y_hist_supply_doutput(56:end),'.',ttt(56:szsh)',y_hist_demand_doutput(56:end),'o',ttt(56:szsh)',y_hist_policy_doutput(56:end),'-.') ;
title('Historical decomposition of annualised output growth');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);

figure(24)
plot(ttt(56:szsh)',[y_dout(56:end)'-0*y_forec_dout(56:szsh)'],'k','LineWidth',3);
hold on
h=bar(ttt(56:szsh)',[y_hist_supply_doutput(56:end)',y_hist_demand_doutput(56:end)',...
    y_hist_policy_doutput(56:end)'],1.5,'grouped');%
set(h(1),'FaceColor','g')
set(h(2),'FaceColor','b')
set(h(3),'FaceColor','r')
hold off
title('Historical decomposition of annualised output growth');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);
axis tight
%%

figure(25)
plot(ttt(56:szsh)',(y_infl(56:end)-y_forec_infl(56:szsh)),'-',ttt(56:szsh)',y_hist_supply_inflobs(56:end),'.',ttt(56:szsh)',y_hist_demand_inflobs(56:end),'o',ttt(56:szsh)',y_hist_policy_inflobs(56:end),'-.') ;
title('Historical decomposition of annualised inflation');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);


figure(26)


plot(ttt(56:szsh-1)',diff([y_infl(56:end)'-0*y_forec_infl(56:szsh)']),'k','LineWidth',3);
hold on
h=bar(ttt(56:szsh-1)',diff([y_hist_supply_inflobs(56:end)',y_hist_demand_inflobs(56:end)',...
    y_hist_policy_inflobs(56:end)']),1.5,'grouped');%
set(h(1),'FaceColor','g')
set(h(2),'FaceColor','b')
set(h(3),'FaceColor','r')
hold off
title('Historical decomposition of annualised inflation');
xlabel('time');
legend('Actual (demeaned)','Supply','Demand','Policy',0);
axis tight
%% plot forward
figure; plot(ttt(start_forward:szsh),forec_for([doutput_pos,inflobs_pos,robs_pos,rr_f_8q_pos],:));
title('Forecast')
legend(listofvariables)
%% saving series

%%Interest
interest_rate_true_a=[4*(oo_.SmoothedVariables.robs(56:end))];
interest_rate_forec_a=(4*forec(robs_pos,56:szsh)')    ;
interest_rate_forec_fwd=(4*forec_for(robs_pos,:)')    ;
true_robs=robs(60:end);
interest_rate_constant=-mean(oo_.SmoothedVariables.robs-[robs(60:end)]);
interest_rate_supply=4*hist_supply_robs(56:end);
interest_rate_demand=4*hist_demand_robs(56:end);
interest_rate_policy=4*hist_policy_robs(56:end);
interest=struct('interest_constant',interest_rate_constant,...
'interest_rate_true_a',interest_rate_true_a,...
'interest_rate_forec_fwd_a',interest_rate_forec_fwd,...
    'interest_rate_forec_a',interest_rate_forec_a,...
    'interest_rate_supply',interest_rate_supply,...
    'interest_rate_demand',interest_rate_demand,...
    'interest_rate_policy',interest_rate_policy,...
    'interest_single_shocks',4*hist_robs(56:end,:));
% natural rate
natural_rate_true_a=[4*(oo_.SmoothedVariables.rr_f_8q(56:end))];
natural_rate_constant=conster;
natural_rate_forec_a=(4*forec(rr_f_8q_pos,56:szsh)')    ;
natural_rate_forec_fwd_a=(4*forec_for(rr_f_8q_pos,:)')    ;

natural_rate_supply=4*hist_supply_rr_f_8q(56:end);
natural_rate_demand=4*hist_demand_rr_f_8q(56:end);
natural_rate_policy=4*hist_policy_rr_f_8q(56:end);
natural_rate=struct('natural_rate_true_a',natural_rate_true_a,...
'natural_rate_constant',natural_rate_constant,...    
'natural_rate_forec_fwd_a',natural_rate_forec_fwd_a,...
'natural_rate_forec_a',natural_rate_forec_a,...
'natural_rate_supply',natural_rate_supply,...
        'natural_rate_demand',natural_rate_demand,...
    'natural_rate_policy',natural_rate_policy,...
    'natural_rate_single_shocks',4*hist_rr_f_8q(56:end,:));


%% Output
output_growth_true_q=y_dout(56:end)';
output_growth_forec_q=forec(doutput_pos,56:szsh)';
output_growth_forec_fwd_q=forec_for(doutput_pos,:)';
trueout=doutput(60:end);
output_growth_constant=-mean(oo_.SmoothedVariables.doutput-doutput(60:end));
output_growth_supply=y_hist_supply_doutput(56:end)';
output_growth_demand=y_hist_demand_doutput(56:end)';
output_growth_policy=y_hist_policy_doutput(56:end)';
output=struct('output_growth_true_q',output_growth_true_q,...
    'output_growth_constant',output_growth_constant,...
    'output_growth_forec_q',output_growth_forec_q,...
    'output_growth_forec_fwd_q',output_growth_forec_fwd_q,...
    'output_growth_supply',output_growth_supply,...
    'output_growth_demand',output_growth_demand,...
    'output_growth_policy',output_growth_policy,...
    'output_growth_single_shocks',hist_doutput(56:end,:));
%% inflation

inflation_year_on_year_true=y_infl(56:end)';
inflation_year_on_year_forec=y_forec_infl(56:szsh)';
inflation_q_on_q_forec_fwd=forec_for(inflobs_pos,:)';

trueinfl=inflobs(60:end);
inflation_yoy_constant=-4*mean(oo_.SmoothedVariables.inflobs-inflobs(60:end));
inflation_yoy_supply=y_hist_supply_inflobs(56:end)';
inflation_yoy_demand=y_hist_demand_inflobs(56:end)';
inflation_yoy_policy=y_hist_policy_inflobs(56:end)';
inflation=struct('inflation_year_on_year_true',inflation_year_on_year_true,...
    'inflation_yoy_constant',inflation_yoy_constant,...
    'inflation_year_on_year_forec',inflation_year_on_year_forec,...
    'inflation_q_on_q_forec_fwd',inflation_q_on_q_forec_fwd,...
    'inflation_yoy_supply',inflation_yoy_supply,...
    'inflation_yoy_demand',inflation_yoy_demand,...
    'inflation_yoy_policy',inflation_yoy_policy,...
    'inflation_yoy_single_shocks',y_hist_all_inflobs(56:end,:));
decomposition=struct('inflation',inflation,'output',output,'interest',interest,'natural',natural_rate);
savename=['decomposition_gianni_',forecyes,datestr(datenum(now),'yymmdd',2),'.mat'];

decomposition.time=ttt(56:szsh);
nameshocks={'technology','preference','government','investment','costpush','policy','wage'};
decomposition.nameshocks=nameshocks;

eval(['save ',savename, ' decomposition -MAT'])
%% Historical shocks
listshocks=strvcat('ea','eb','eg','ei','ep','er','ew');
clear hist_shocks
for jj=1:size(lgx_,1);
    eval(['hist_shocks.',deblank(char(nameshocks{jj})),'=oo_.SmoothedVariables.',deblank(listshocks(jj,:)),'(56:end);']);
end
savename=['hist_shocks_gianni_',datestr(datenum(now),'yymmdd',2),'.mat'];
eval(['save ',savename, ' hist_shocks -MAT'])

% allseries=[interest_rate_true_a,interest_rate_forec_a,interest_rate_supply,interest_rate_demand,...
%     interest_rate_policy,output_growth_true_q,output_growth_forec_q,output_growth_supply,...
%     output_growth_demand,output_growth_policy,inflation_year_on_year_true,inflation_year_on_year_forec,...
%     inflation_yoy_supply,inflation_yoy_demand,inflation_yoy_policy,...
%     natural_rate_true_a,natural_rate_supply,natural_rate_demand,...
%     natural_rate_policy];
% dlmwrite('decomposition.dat',allseries,',');

%%%%%%%%%%%% correlation of shocks (innovatios)
shc=corr(sh_t);
