

% Variance decomposition Dynare

clear all
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
scriptsize='\scriptsize';% typ-size for latex
% if ~exist('\scriptsize','var')
%     scriptsize='scriptsize';
% elseif isempty(scriptsize);
%     scriptsize='\scriptsize';
datafile=input('Give name of mat file with data ==> ','s');

if ~isempty(datafile);
    data_=load(datafile);
   
end
% 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 isfield([options_],'datafile')
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;
end

varlist=input('Give name of file with selected variables ','s');
if isempty(varlist)|~exist(varlist,'file');
    listofvariables=str2cell(lgy_);
else
listofvariables=textread(varlist,'%s\n','whitespace','\n','commentstyle','matlab');;%strvcat('doutput','inflobs','robs','rr_f_8q');
end
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)];

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


% State stpace form
% y=A_var*y(-1)+Sig_var*epsilon;
Sigma_e=M_.Sigma_e;
horizons=[1,4,20,40,100,500];
% ergo=500;
vardec=zeros(size(A_var,1),size(Sig_var,2));
zervar=zeros(size(Sigma_e));
z=1;
for jj=horizons
    
    tot_var=zeros(size(A_var));
    for uu=1:jj
        tot_var=A_var*tot_var*A_var'+Sig_var*Sigma_e*Sig_var';
    end
   tot_vard=diag(tot_var);
   tot_vard(tot_vard<1e-10)=NaN;
for kk=1:size(Sig_var,2);
    zervar=zeros(size(Sigma_e));
    zervar(kk,kk)=Sigma_e(kk,kk);
    tmp=zeros(size(A_var));
    
    for uu=1:jj
    tmp=A_var*tmp*A_var'+Sig_var*zervar*Sig_var';
    end
    vardec(:,kk,z)=diag(tmp)./(tot_vard);
    
end
z=z+1;
end 

vardec_sub=vardec(length([nstat+[1:npred]])+1:end,:,:);

z=1;
for jj=horizons
dyntable(['Variance decomposition ',num2str(jj),' periods ahead'],...
    strvcat(' ',shocks),strvcat(listofvariables),vardec_sub(:,:,z)*100,10,6, ...
	       1);
       z=z+1;
end



prvdec=input(upper('LATEX VARIANCE DECOMPOSITIONS? (y/n(default)) :'),'s');
if strcmpi(prvdec,'y');
if ~exist('tex_stuff','dir');
    mkdir('tex_stuff');
end
    z=1;
for jj=horizons
    tb=tabtex_highlight(strvcat('variables ',strvcat(regexprep(listofvariables,'\_','\\_'))),['Variance decomposition ',num2str(jj),' period(s) ahead'],3,...
        ones(size(vardec_sub,2),1)*20,'abs','>',(strvcat(regexprep(str2cell(shocks),'\_','\\_'))),vardec_sub(:,:,z)'*100);
tabtex2texpdf(['vardec_',num2str(jj),'_pa'],1,'',tb,0,scriptsize);
z=z+1;
end
end
%%%%%%%%%%

data=[];
obs_=[];
zz=1;
if exist('data_','var')
for uu=1:length(listofvariables)
    if isfield(data_,listofvariables{uu});
        data(zz,:)=eval(['(data_.',listofvariables{uu},')']);
        obs_{zz,1}=listofvariables{uu};
        zz=zz+1;
        
    end
end
% moments
ergo=1000;

tot_var=zeros(size(A_var));
    for uu=1:ergo
        tot_var=A_var*tot_var*A_var'+Sig_var*Sigma_e*Sig_var';
    end
    
totvar_sub=tot_var(length([nstat+[1:npred]])+1:end,length([nstat+[1:npred]])+1:end);
correl=totvar_sub./sqrt(diag(totvar_sub)*diag(totvar_sub)');
correl_data=corr(data');
table(['Correlation '],...
    strvcat(' ',strvcat(listofvariables)),strvcat(listofvariables),correl,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(listofvariables,'\_','\\_'))),['Correlations Model'],3,...
        (strvcat(regexprep(listofvariables,'\_','\\_'))),triu(correl));
tabtex2texpdf(['correlations_model'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\tiny');

tb=tabtex(strvcat('variables ',strvcat(regexprep(obs_,'\_','\\_'))),['Correlations Data'],3,...
        (strvcat(regexprep(obs_,'\_','\\_'))),triu(correl_data));
tabtex2texpdf(['correlations_data'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\tiny');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nrows1=round(length(listofvariables)/2);
nrows2=length(listofvariables)-nrows1;
% divide the corr-matrix in three blocks


tmp1=triu(correl(1:nrows1,1:nrows1));
tmp2=(correl(1:nrows1,nrows1+1:end));
tmp3=triu(correl(nrows1+1:end,nrows1+1:end));
varlist1=listofvariables(1:nrows1,:);
varlist2=listofvariables(nrows1+1:end,:);
%1)
table(['Correlation Model'],...
    strvcat(' ',strvcat(varlist1)),strvcat(varlist1),tmp1,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(varlist1,'\_','\\_'))),['Correlations Model'],1,...
        (strvcat(regexprep(varlist1,'\_','\\_'))),tmp1);
tabtex2texpdf(['correlations_model1'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');


table(['Correlation Model'],...
    strvcat(' ',strvcat(varlist2)),strvcat(varlist1),tmp2,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(varlist2,'\_','\\_'))),['Correlations Model'],1,...
        (strvcat(regexprep(varlist1,'\_','\\_'))),tmp2);
tabtex2texpdf(['correlations_model2'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');


table(['Correlation Model'],...
    strvcat(' ',strvcat(varlist2)),strvcat(varlist2),tmp3,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(varlist2,'\_','\\_'))),['Correlations Model'],1,...
        (strvcat(regexprep(varlist2,'\_','\\_'))),tmp3);
tabtex2texpdf(['correlations_model3'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nrows1=round(length(obs_)/2);
nrows2=length(obs_)-nrows1;
% divide the corr-matrix in three blocks


tmp1=triu(correl_data(1:nrows1,1:nrows1));
tmp2=(correl_data(1:nrows1,nrows1+1:end));
tmp3=triu(correl_data(nrows1+1:end,nrows1+1:end));
varlist1=obs_(1:nrows1,:);
varlist2=obs_(nrows1+1:end,:);
%1)
table(['Correlation Data'],...
    strvcat(' ',strvcat(varlist1)),strvcat(varlist1),tmp1,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(varlist1,'\_','\\_'))),['Correlations Data'],1,...
        (strvcat(regexprep(varlist1,'\_','\\_'))),tmp1);
tabtex2texpdf(['correlations_data1'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');


table(['Correlation Data'],...
    strvcat(' ',strvcat(varlist2)),strvcat(varlist1),tmp2,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(varlist2,'\_','\\_'))),['Correlations Data'],1,...
        (strvcat(regexprep(varlist1,'\_','\\_'))),tmp2);
tabtex2texpdf(['correlations_data2'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');


table(['Correlation Data'],...
    strvcat(' ',strvcat(varlist2)),strvcat(varlist2),tmp3,8,6, ...
	       2);
 tb=tabtex(strvcat('variables ',strvcat(regexprep(varlist2,'\_','\\_'))),['Correlations Data'],1,...
        (strvcat(regexprep(varlist2,'\_','\\_'))),tmp3);
tabtex2texpdf(['correlations_data3'],1,'',cell2str(regexprep(str2cell(tb),'\<0&\>','--&')),0,'\scriptsize');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


standdev=sqrt(diag(totvar_sub));
standdev_data=zeros(length(standdev),1)*NaN;

for uu=1:length(listofvariables)
    if isfield(data_,listofvariables{uu});
        standdev_data(uu,1)=eval(['std(data_.',listofvariables{uu},')']);
    end
end

table(['Standard Deviations '],...
    strvcat('Variables ','std'),strvcat(strvcat(listofvariables)),standdev,8,6, ...
	       2);
 tb=tabtex(strvcat('Variables ','std-model','std-data'),['Standard Deviations'],3,...
        (strvcat(regexprep(listofvariables,'\_','\\_'))),[standdev,standdev_data]);
tabtex2texpdf(['standard_dev'],1,'',tb,0,scriptsize);

end %if exist data_