1function [nvar,vartan,NumberOfConditionalDecompFiles] = ...
2    dsge_simulated_theoretical_conditional_variance_decomposition(SampleSize,Steps,M_,options_,oo_,type)
3% function [nvar,vartan,NumberOfConditionalDecompFiles] = ...
4%     dsge_simulated_theoretical_conditional_variance_decomposition(SampleSize,Steps,M_,options_,oo_,type)
5% This function computes the posterior or prior distribution of the conditional variance
6% decomposition of the endogenous variables (or a subset of the endogenous variables).
7%
8% INPUTS
9%   SampleSize   [integer]       scalar, number of simulations.
10%   Steps        [integers]      horizons at which to conduct decomposition
11%   M_           [structure]     Dynare structure describing the model.
12%   options_     [structure]     Dynare structure defining global options.
13%   oo_          [structure]     Dynare structure where the results are saved.
14%   type         [string]        'prior' or 'posterior'
15%
16%
17% OUTPUTS
18%   nvar                             [integer]  nvar is the number of stationary variables.
19%   vartan                           [char]     array of characters (with nvar rows).
20%   NumberOfConditionalDecompFiles   [integer]  scalar, number of prior or posterior data files (for covariance).
21
22% Copyright (C) 2009-2020 Dynare Team
23%
24% This file is part of Dynare.
25%
26% Dynare is free software: you can redistribute it and/or modify
27% it under the terms of the GNU General Public License as published by
28% the Free Software Foundation, either version 3 of the License, or
29% (at your option) any later version.
30%
31% Dynare is distributed in the hope that it will be useful,
32% but WITHOUT ANY WARRANTY; without even the implied warranty of
33% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
34% GNU General Public License for more details.
35%
36% You should have received a copy of the GNU General Public License
37% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
38
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_conditional_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 '_PosteriorConditionalVarianceDecomposition*'])
56    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecompositionME*'])
57else
58    delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition*'])
59    delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecompositionME*'])
60end
61
62% Set varlist (vartan)
63if ~posterior
64    if isfield(options_,'varlist')
65        temp = options_.varlist;
66    end
67    options_.varlist = options_.prior_analysis_endo_var_list;
68end
69[ivar,vartan ] = get_variables_list(options_, M_);
70if ~posterior
71    if exist('temp','var')
72        options_.varlist = temp;
73    end
74end
75nvar = length(ivar);
76
77% Set the size of the auto-correlation function to zero.
78nar = options_.ar;
79options_.ar = 0;
80
81NumberOfSavedElementsPerSimulation = nvar*M_.exo_nbr*length(Steps);
82MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
83
84ME_present=0;
85if ~all(diag(M_.H)==0)
86    if (isoctave && octave_ver_less_than('6')) || (~isoctave && matlab_ver_less_than('8.1'))
87        [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
88    else
89        [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
90    end
91    if ~isempty(observable_pos_requested_vars)
92        ME_present=1;
93        nobs_ME=length(observable_pos_requested_vars);
94        NumberOfSavedElementsPerSimulation_ME = nobs_ME*(M_.exo_nbr+1)*length(Steps);
95        MaXNumberOfConditionalDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
96    end
97end
98
99if SampleSize<=MaXNumberOfConditionalDecompLines
100    Conditional_decomposition_array = zeros(nvar,length(Steps),M_.exo_nbr,SampleSize);
101    NumberOfConditionalDecompFiles = 1;
102else
103    Conditional_decomposition_array = zeros(nvar,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines);
104    NumberOfLinesInTheLastConditionalDecompFile = mod(SampleSize,MaXNumberOfConditionalDecompLines);
105    NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfConditionalDecompLines);
106end
107
108if ME_present
109    if SampleSize<=MaXNumberOfConditionalDecompLines_ME
110        Conditional_decomposition_array_ME = zeros(nobs_ME,length(Steps),M_.exo_nbr+1,SampleSize);
111        NumberOfConditionalDecompFiles_ME = 1;
112    else
113        Conditional_decomposition_array_ME = zeros(nobs_ME,length(Steps),M_.exo_nbr+1,SampleSize);
114        NumberOfLinesInTheLastConditionalDecompFile_ME = mod(SampleSize,MaXNumberOfConditionalDecompLines_ME);
115        NumberOfConditionalDecompFiles_ME = ceil(SampleSize/MaXNumberOfConditionalDecompLines_ME);
116    end
117    NumberOfConditionalDecompLines_ME = size(Conditional_decomposition_array_ME,4);
118    ConditionalDecompFileNumber_ME = 0;
119end
120
121NumberOfConditionalDecompLines = size(Conditional_decomposition_array,4);
122ConditionalDecompFileNumber = 0;
123
124
125StateSpaceModel.number_of_state_equations = M_.endo_nbr;
126StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
127
128first_call = 1;
129
130linea = 0;
131linea_ME = 0;
132for file = 1:NumberOfDrawsFiles
133    if posterior
134        temp=load([M_.dname '/metropolis/' M_.fname '_' type '_draws' num2str(file) ]);
135    else
136        temp=load([M_.dname '/prior/draws/' type '_draws' num2str(file) ]);
137    end
138    isdrsaved = columns(temp.pdraws)-1;
139    NumberOfDraws = rows(temp.pdraws);
140    for linee = 1:NumberOfDraws
141        linea = linea+1;
142        linea_ME = linea_ME+1;
143        if isdrsaved
144            M_=set_parameters_locally(M_,temp.pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
145            dr = temp.pdraws{linee,2};
146        else
147            M_=set_parameters_locally(M_,temp.pdraws{linee,1});
148            [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
149        end
150        if first_call
151            endo_nbr = M_.endo_nbr;
152            nstatic = M_.nstatic;
153            nspred = M_.nspred;
154            iv = (1:endo_nbr)';
155            ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
156            StateSpaceModel.number_of_state_equations = M_.endo_nbr;
157            StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
158            StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal;
159            StateSpaceModel.order_var = dr.order_var;
160            StateSpaceModel.observable_pos=options_.varobs_id;
161            first_call = 0;
162            clear('endo_nbr','nstatic','nspred','k');
163        end
164        [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,M_.exo_nbr);
165        StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
166        M_ = set_measurement_errors(temp.pdraws{linee,1},temp.estim_params_,M_);
167        StateSpaceModel.measurement_error=M_.H;
168        clear('dr');
169        [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]=conditional_variance_decomposition(StateSpaceModel, Steps, ivar);
170        Conditional_decomposition_array(:,:,:,linea) =ConditionalVarianceDecomposition;
171        if ME_present
172            Conditional_decomposition_array_ME(:,:,:,linea) =ConditionalVarianceDecomposition_ME;
173        end
174        if linea == NumberOfConditionalDecompLines
175            ConditionalDecompFileNumber = ConditionalDecompFileNumber + 1;
176            linea = 0;
177            if posterior
178                save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
179                     'Conditional_decomposition_array');
180            else
181                save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition' int2str(ConditionalDecompFileNumber) '.mat' ], ...
182                     'Conditional_decomposition_array');
183            end
184            if (ConditionalDecompFileNumber==NumberOfConditionalDecompFiles-1)% Prepare last round.
185                Conditional_decomposition_array = zeros(nvar, length(Steps),M_.exo_nbr,NumberOfLinesInTheLastConditionalDecompFile) ;
186                NumberOfConditionalDecompLines = NumberOfLinesInTheLastConditionalDecompFile;
187            elseif ConditionalDecompFileNumber<NumberOfConditionalDecompFiles-1
188                Conditional_decomposition_array = zeros(nvar,length(Steps),M_.exo_nbr,MaXNumberOfConditionalDecompLines);
189            else
190                clear('Conditional_decomposition_array');
191            end
192        end
193        %with measurement error
194        if ME_present
195            if linea_ME == NumberOfConditionalDecompLines_ME
196                ConditionalDecompFileNumber_ME = ConditionalDecompFileNumber_ME + 1;
197                linea_ME = 0;
198                if posterior
199                    save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
200                         'Conditional_decomposition_array_ME');
201                else
202                    save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
203                         'Conditional_decomposition_array_ME');
204                end
205                if (ConditionalDecompFileNumber_ME==NumberOfConditionalDecompFiles_ME-1)% Prepare last round.
206                    Conditional_decomposition_array_ME = zeros(nobs_ME, length(Steps),M_.exo_nbr+1,NumberOfLinesInTheLastConditionalDecompFile_ME) ;
207                    NumberOfConditionalDecompLines_ME = NumberOfLinesInTheLastConditionalDecompFile_ME;
208                elseif ConditionalDecompFileNumber_ME<NumberOfConditionalDecompFiles_ME-1
209                    Conditional_decomposition_array_ME = zeros(nobs_ME,length(Steps),M_.exo_nbr+1,MaXNumberOfConditionalDecompLines_ME);
210                else
211                    clear('Conditional_decomposition_array_ME');
212                end
213            end
214        end
215    end
216end
217
218options_.ar = nar;
219