1function [nvar,vartan,NumberOfDecompFiles] = ...
2    dsge_simulated_theoretical_variance_decomposition(SampleSize,M_,options_,oo_,type)
3% function [nvar,vartan,NumberOfDecompFiles] = ...
4%     dsge_simulated_theoretical_variance_decomposition(SampleSize,M_,options_,oo_,type)
5% This function computes the posterior or prior distribution of the variance
6% decomposition of the observed endogenous variables.
7%
8% INPUTS
9%   SampleSize   [integer]       scalar, number of simulations.
10%   M_           [structure]     Dynare structure describing the model.
11%   options_     [structure]     Dynare structure defining global options.
12%   oo_          [structure]     Dynare structure where the results are saved.
13%   type         [string]        'prior' or 'posterior'
14%
15%
16% OUTPUTS
17%   nvar              [integer]  nvar is the number of stationary variables.
18%   vartan            [char]     array of characters (with nvar rows).
19%   CovarFileNumber   [integer]  scalar, number of prior or posterior data files (for covariance).
20
21% Copyright (C) 2007-2020 Dynare Team
22%
23% This file is part of Dynare.
24%
25% Dynare is free software: you can redistribute it and/or modify
26% it under the terms of the GNU General Public License as published by
27% the Free Software Foundation, either version 3 of the License, or
28% (at your option) any later version.
29%
30% Dynare is distributed in the hope that it will be useful,
31% but WITHOUT ANY WARRANTY; without even the implied warranty of
32% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
33% GNU General Public License for more details.
34%
35% You should have received a copy of the GNU General Public License
36% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
37
38nodecomposition = 0;
39
40% Get informations about the _posterior_draws files.
41if strcmpi(type,'posterior')
42    NumberOfDrawsFiles = length(dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]));
43    posterior = 1;
44elseif strcmpi(type,'prior')
45    NumberOfDrawsFiles = length(dir([M_.dname '/prior/draws/' type '_draws*' ]));
46    CheckPath('prior/moments',M_.dname);
47    posterior = 0;
48else
49    disp('dsge_simulated_theoretical_variance_decomposition:: Unknown type!')
50    error()
51end
52
53%delete old stale files before creating new ones
54if posterior
55    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecomposition*']);
56    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME*']);
57else
58    delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecomposition*']);
59    delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME*']);end
60
61% Set varlist (vartan)
62if ~posterior
63    if isfield(options_,'varlist')
64        temp = options_.varlist;
65    end
66    options_.varlist = options_.prior_analysis_endo_var_list;
67end
68[ivar,vartan,options_] = get_variables_list(options_,M_);
69if ~posterior
70    if exist('temp','var')
71        options_.varlist = temp;
72    end
73end
74nvar = length(ivar);
75
76% Set the size of the auto-correlation function to zero.
77nar = options_.ar;
78options_.ar = 0;
79
80
81
82nexo = M_.exo_nbr;
83
84NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
85MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
86
87ME_present=0;
88if ~all(diag(M_.H)==0)
89    if (isoctave && octave_ver_less_than('6')) || (~isoctave && matlab_ver_less_than('8.1'))
90        [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
91    else
92        [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
93    end
94    if ~isempty(observable_pos_requested_vars)
95        ME_present=1;
96        nobs_ME=length(observable_pos_requested_vars);
97        NumberOfSavedElementsPerSimulation_ME = nobs_ME*(nexo+1);
98        MaXNumberOfDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
99    end
100end
101
102if SampleSize<=MaXNumberOfDecompLines
103    Decomposition_array = zeros(SampleSize,nvar*nexo);
104    NumberOfDecompFiles = 1;
105else
106    Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*nexo);
107    NumberOfLinesInTheLastDecompFile = mod(SampleSize,MaXNumberOfDecompLines);
108    NumberOfDecompFiles = ceil(SampleSize/MaXNumberOfDecompLines);
109end
110
111NumberOfDecompLines = rows(Decomposition_array);
112DecompFileNumber = 1;
113
114if ME_present
115    if SampleSize<=MaXNumberOfDecompLines_ME
116        Decomposition_array_ME = zeros(SampleSize,nobs_ME*(nexo+1));
117        NumberOfDecompFiles_ME = 1;
118    else
119        Decomposition_array_ME = zeros(MaXNumberOfDecompLines_ME,nobs_ME*(nexo+1));
120        NumberOfLinesInTheLastDecompFile_ME = mod(SampleSize,MaXNumberOfDecompLines_ME);
121        NumberOfDecompFiles_ME = ceil(SampleSize/MaXNumberOfDecompLines_ME);
122    end
123    NumberOfDecompLines_ME = rows(Decomposition_array_ME);
124    DecompFileNumber_ME = 1;
125end
126% Compute total variances (covariances are not saved) and variances
127% implied by each structural shock.
128linea = 0;
129linea_ME = 0;
130only_non_stationary_vars=0;
131for file = 1:NumberOfDrawsFiles
132    if posterior
133        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
134    else
135        temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
136    end
137    isdrsaved = columns(temp.pdraws)-1;
138    NumberOfDraws = rows(temp.pdraws);
139    for linee = 1:NumberOfDraws
140        linea = linea+1;
141        linea_ME = linea_ME+1;
142        if isdrsaved
143            M_=set_parameters_locally(M_,temp.pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
144            dr = temp.pdraws{linee,2};
145        else
146            M_=set_parameters_locally(M_,temp.pdraws{linee,1});
147            [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
148        end
149        if file==1 && linee==1
150            [tmp, stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
151            if isempty(stationary_vars)
152                fprintf('\ndsge_simulated_theoretical_variance_decomposition:: All requested endogenous variables have a unit root and thus infinite variance.\n')
153                fprintf('dsge_simulated_theoretical_variance_decomposition:: No decomposition is performed.\n')
154                only_non_stationary_vars=1;
155            end
156        end
157        if only_non_stationary_vars
158            Decomposition_array(linea,:) = NaN;
159        else
160            tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
161            for i=1:nvar
162                for j=1:nexo
163                    Decomposition_array(linea,(i-1)*nexo+j) = tmp{2}(i,j);
164                end
165            end
166            if ME_present
167                M_ = set_measurement_errors(temp.pdraws{linee,1},temp.estim_params_,M_);
168                ME_Variance=diag(M_.H);
169                tmp_ME=NaN(nobs_ME,nexo+1);
170                tmp_ME(:,1:end-1)=tmp{2}(index_subset,:).*repmat(diag(tmp{1}(index_subset,index_subset))./(diag(tmp{1}(index_subset,index_subset))+ME_Variance(index_observables)),1,nexo);
171                tmp_ME(:,end)=1-sum(tmp_ME(:,1:end-1),2);
172                for i=1:nobs_ME
173                    for j=1:nexo+1
174                        Decomposition_array_ME(linea,(i-1)*(nexo+1)+j) = tmp_ME(i,j);
175                    end
176                end
177            end
178        end
179        if linea == NumberOfDecompLines
180            if posterior
181                save([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecomposition' int2str(DecompFileNumber) '.mat' ],'Decomposition_array');
182            else
183                save([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecomposition' int2str(DecompFileNumber) '.mat' ],'Decomposition_array');
184            end
185            DecompFileNumber = DecompFileNumber + 1;
186            linea = 0;
187            test = DecompFileNumber-NumberOfDecompFiles;
188            if ~test% Prepare the last round...
189                Decomposition_array = zeros(NumberOfLinesInTheLastDecompFile,nvar*nexo);
190                NumberOfDecompLines = NumberOfLinesInTheLastDecompFile;
191            elseif test<0
192                Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*nexo);
193            else
194                clear('Decomposition_array');
195            end
196        end
197        if ME_present
198            if linea_ME == NumberOfDecompLines_ME
199                if posterior
200                    save([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME');
201                else
202                    save([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME');
203                end
204                DecompFileNumber_ME = DecompFileNumber_ME + 1;
205                linea_ME = 0;
206                test = DecompFileNumber_ME-NumberOfDecompFiles_ME;
207                if ~test% Prepare the last round...
208                    Decomposition_array_ME = zeros(NumberOfLinesInTheLastDecompFile_ME,nobs_ME*(nexo+1));
209                    NumberOfDecompLines_ME = NumberOfLinesInTheLastDecompFile_ME;
210                elseif test<0
211                    Decomposition_array_ME = zeros(MaXNumberOfDecompLines_ME,nobs_ME*(nexo+1));
212                else
213                    clear('Decomposition_array_ME');
214                end
215            end
216        end
217    end
218end
219
220options_.ar = nar;
221