Description of prior

0001 function results = prior_sampler(drsave,M_,bayestopt_,options_,oo_)
0002 % This function builds a (big) prior sample.
0003 %
0004 % INPUTS
0005 %   drsave      [integer]    Scalar. If equal to 1, then dr structure is saved with each prior draw.
0006 %   M_          [structure]  Model description.
0007 %   bayestopt_  [structure]  Prior distribution description.
0008 %   options_    [structure]  Global options of Dynare.
0009 %
0010 % OUTPUTS:
0011 %   results     [structure]  Various statistics.
0012 %
0013 % SPECIAL REQUIREMENTS
0014 %   none
0015 
0016 % Copyright (C) 2009-2010 Dynare Team
0017 %
0018 % This file is part of Dynare.
0019 %
0020 % Dynare is free software: you can redistribute it and/or modify
0021 % it under the terms of the GNU General Public License as published by
0022 % the Free Software Foundation, either version 3 of the License, or
0023 % (at your option) any later version.
0024 %
0025 % Dynare is distributed in the hope that it will be useful,
0026 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0027 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0028 % GNU General Public License for more details.
0029 %
0030 % You should have received a copy of the GNU General Public License
0031 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0032 
0033 % Initialization.
0034 prior_draw(1);
0035 PriorDirectoryName = CheckPath('prior/draws',M_.dname);
0036 work = ~drsave;
0037 iteration = 0;
0038 loop_indx = 0;
0039 file_indx = [];
0040 count_bk_indeterminacy = 0;
0041 count_bk_unstability = 0;
0042 count_bk_singularity = 0;
0043 count_static_var_def = 0;
0044 count_no_steadystate = 0;
0045 count_steadystate_file_exit = 0;
0046 count_dll_problem = 0;
0047 count_complex_jacobian = 0;
0048 count_complex_steadystate = 0;
0049 count_nan_steadystate = 0;
0050 count_nan_params = 0;
0051 count_complex_params = 0;
0052 count_unknown_problem = 0;
0053 NumberOfSimulations = options_.prior_mc;
0054 NumberOfParameters = length(bayestopt_.p1);
0055 NumberOfEndogenousVariables = size(M_.endo_names,1);
0056 NumberOfElementsPerFile = ceil(options_.MaxNumberOfBytes/NumberOfParameters/NumberOfEndogenousVariables/8) ;
0057 
0058 if NumberOfSimulations <= NumberOfElementsPerFile
0059     TableOfInformations = [ 1 ,  NumberOfSimulations , 1] ;
0060 else
0061     NumberOfFiles = ceil(NumberOfSimulations/NumberOfElementsPerFile) ;
0062     NumberOfElementsInTheLastFile = NumberOfSimulations - NumberOfElementsPerFile*(NumberOfFiles-1) ;
0063     TableOfInformations = NaN(NumberOfFiles,3) ;
0064     TableOfInformations(:,1) = transpose(1:NumberOfFiles) ;
0065     TableOfInformations(1:NumberOfFiles-1,2) = NumberOfElementsPerFile*ones(NumberOfFiles-1,1) ;
0066     TableOfInformations(NumberOfFiles,2) = NumberOfElementsInTheLastFile ;
0067     TableOfInformations(1,3) = 1;
0068     TableOfInformations(2:end,3) = cumsum(TableOfInformations(1:end-1,2))+1;
0069 end
0070 
0071 pdraws = cell(TableOfInformations(1,2),drsave+1) ;
0072 sampled_prior_expectation = zeros(NumberOfParameters,1);
0073 sampled_prior_covariance  = zeros(NumberOfParameters,NumberOfParameters);
0074 
0075 file_line_number = 0;
0076 file_indx_number = 0;
0077 
0078 % Simulations.
0079 while iteration < NumberOfSimulations
0080     loop_indx = loop_indx+1;
0081     params = prior_draw();
0082     set_all_parameters(params);
0083     [dr,INFO,M_,options_,oo_] = resol(work,M_,options_,oo_);
0084     switch INFO(1)
0085       case 0
0086         file_line_number = file_line_number + 1 ;
0087         iteration = iteration + 1;
0088         pdraws(file_line_number,1) = {params};
0089         if drsave
0090             pdraws(file_line_number,2) = {dr};
0091         end
0092         [sampled_prior_expectation,sampled_prior_covariance] = ...
0093             recursive_prior_moments(sampled_prior_expectation,sampled_prior_covariance,params,iteration);
0094       case 1
0095         count_static_undefined = count_static_undefined + 1;
0096       case 2
0097         count_dll_problem = count_dll_problem + 1;
0098       case 3
0099         count_bk_unstability = count_bk_unstability + 1 ;
0100       case 4
0101         count_bk_indeterminacy = count_bk_indeterminacy + 1 ;
0102       case 5
0103         count_bk_singularity = count_bk_singularity + 1 ;
0104       case 20
0105         count_no_steadystate = count_no_steadystate + 1 ;
0106       case 19
0107         count_steadystate_file_exit = count_steadystate_file_exit + 1 ;
0108       case 6
0109         count_complex_jacobian = count_complex_jacobian + 1 ;
0110       case 21
0111         count_complex_steadystate = count_complex_steadystate + 1 ;
0112       case 22
0113         count_nan_steadystate = count_nan_steadystate + 1 ;
0114       case 23
0115         count_complex_params = count_complex_params + 1 ;
0116       case 24
0117         count_nan_params = count_nan_params + 1 ;
0118       otherwise
0119         count_unknown_problem = count_unknown_problem + 1 ;
0120     end
0121     if ( file_line_number==TableOfInformations(file_indx_number+1,2) )
0122         file_indx_number = file_indx_number + 1;
0123         save([ PriorDirectoryName '/prior_draws' int2str(file_indx_number) '.mat' ],'pdraws');
0124         if file_indx_number<NumberOfFiles
0125             pdraws = cell(TableOfInformations(file_indx_number+1,2),drsave+1);
0126         end
0127         file_line_number = 0;
0128     end
0129 end
0130 
0131 % Get informations about BK conditions and other things...
0132 results.bk.indeterminacy_share = count_bk_indeterminacy/loop_indx;
0133 results.bk.unstability_share = count_bk_unstability/loop_indx;
0134 results.bk.singularity_share = count_bk_singularity/loop_indx;
0135 results.dll.problem_share = count_dll_problem/loop_indx;
0136 results.ss.problem_share = count_no_steadystate/loop_indx;
0137 results.ss.complex_share = count_complex_steadystate/loop_indx;
0138 results.ass.problem_share = count_steadystate_file_exit/loop_indx;
0139 results.jacobian.problem_share = count_complex_jacobian/loop_indx;
0140 results.garbage_share = ...
0141     results.bk.indeterminacy_share + ...
0142     results.bk.unstability_share + ...
0143     results.bk.singularity_share + ...
0144     results.dll.problem_share + ...
0145     results.ss.problem_share + ...
0146     results.ass.problem_share + ...
0147     results.jacobian.problem_share + ...
0148     count_unknown_problem/loop_indx ;
0149 results.prior.mean = sampled_prior_expectation;
0150 results.prior.variance = sampled_prior_covariance;
0151 results.prior.mass = 1-results.garbage_share;
0152 
0153 function [mu,sigma] = recursive_prior_moments(m0,s0,newobs,iter)
0154 %  Recursive estimation of order one and two moments (expectation and
0155 %  covariance matrix). newobs should be a row vector. I do not use the
0156 %  function recursive_moments here, because this function is to be used when
0157 %  newobs is a 2D array.
0158 m1 = m0 + (newobs'-m0)/iter;
0159 qq = m1*m1';
0160 s1 = s0 + ( (newobs'*newobs-qq-s0) + (iter-1)*(m0*m0'-qq') )/iter;
0161 mu = m1;
0162 sigma = s1;
prior_sampler

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
