% Code to do the variance decomposition exercise utilizing stimul_ in Dynare V4, following Christiano
%close all;
clear Etechtab
global endo_nbr ykmin_ xkmin_ it_ options_ iy_
data_primeceri10;

%% User's Control
a=8;
b=32;
startp=21;
hp_lambda=1600;

%end of User's Control
%% Get the matrix of shocks
var_names=M_.endo_names;
shocknames=M_.exo_names;
if ~isempty(oo_.SmoothedShocks);
    for iter=1:size(shocknames,1);
        name=['oo_.SmoothedShocks.' shocknames(iter,:)];
        shocks(:,iter)=eval(name); %Attention: the last row of shocks is zero, what does that mean?
    end;
end;
dr=oo_.dr;
fstobs=options_.first_obs;
num_var=length(dr.ys);
%% Get the endogenous variables
var_interested=char('dlogy','dlogc', 'dlogi', 'logHours', 'dlogRealW', 'RR_obs', 'Infl','y_t','h_t','Rn_t');
var_names=M_.endo_names;
var_obs=[num_var-7:num_var]';
for iter=1:8;
pickvar(iter)=find(dr.order_var==var_obs(iter));
end
shock=shocks(1:end-1,:);%Remove the last row, since it is always zero;
 [yy,k2] = simult_jh(zeros(num_var,1),oo_.dr,shock,1);
for iter=1:size(shock,2);
    u=zeros(size(shock));
    u(:,iter)=shock(:,iter); %Get the input of shocks
    
    structnames=fieldnames(oo_.SmoothedVariables);
    yt=zeros(size(structnames,1),1); %%????
    for ii=1:size(structnames,1);
        yt(ii)=oo_.SmoothedVariables.(char(structnames(ii,:)))(1,1); %Initial conditions
    end;
    
    yts=yt(oo_.dr.order_var)+oo_.dr.ys; %% ?? Why does he use inverse order ?
    iorder=options_.order;
    ex=u;
     yy = simult_(yt,oo_.dr,ex,iorder);
     ey=shock;
     yyall=simult_(yt,oo_.dr,ey,iorder);
 %% Plot the Historic Decomposition;
     figure;
     time=(1954.0:.25:2004.75);
     for itervar=1:length(pickvar);
subplot(3,3,itervar);
tmpind=oo_.dr.order_var(pickvar(itervar));
% startp
% size(time(startp:end))
% size(yy(tmpind,startp:end))
% pause
plot(time(startp:end),yy(tmpind,startp-fstobs+1:end))
hold on;
plot(time(startp:end),yyall(tmpind,startp-fstobs+1:end),'r-')
plot(time(startp:end)',vdata(startp:end,itervar)-mean(vdata(startp:end,itervar)),'g-');
axis tight;
title(structnames(pickvar(itervar)))
%title(structnames(itervar));
%% Frequency Decomposition
hp1    =   yy(tmpind,startp:end)'-HPF(yy(tmpind,startp:end)',1600);
hp5    =   vdata(startp:end,itervar)-HPF( vdata(startp:end,itervar),1600);
bp1    =   bpassgen(yy(tmpind,startp:end)',a,b);
bp5    =   bpassgen( vdata(startp:end,itervar),a,b);
bp1lr   =   bpassgen(yy(tmpind,startp:end)',b,80);
bp5lr   =   bpassgen( vdata(startp:end,itervar),b,80);
bp1sr   =   bpassgen(yy(tmpind,startp:end)',2,a);
bp5sr   =   bpassgen( vdata(startp:end,itervar),2,a);
hpEtech  =   100*var(hp1)./var(hp5);
bpEtech  =   100*var(bp1)./var(bp5);
bpEtechLR   =   100*var(bp1lr)./var(bp5lr);
bpEtechSR   =   100*var(bp1sr)./var(bp5sr);



Etechtab(itervar,:,iter)    =   [hpEtech bpEtech bpEtechLR bpEtechSR];
end;
suptitle(['Response to ' shocknames(iter,:)]);


end;


%If we use the historical decomposition, correlated shocks may generate
%excess volatility 
swsp_util_habit5_vardecomp_forecast