prior_draw

 This function generate one draw from the joint prior distribution.
 
 INPUTS 
   o init             [integer]    scalar equal to 1 (first call) or 0.
   o uniform          [integer]    scalar equal to 1 (first call) or 0.
    
 OUTPUTS 
   o pdraw            [double]     1*npar vector, draws from the joint prior density.


 SPECIAL REQUIREMENTS
   none

 NOTE 1. Input arguments 1 an 2 are only needed for initialization.
 NOTE 2. A given draw from the joint prior distribution does not satisfy BK conditions a priori.

0001 function pdraw = prior_draw(init,uniform)
0002 % This function generate one draw from the joint prior distribution.
0003 %
0004 % INPUTS
0005 %   o init             [integer]    scalar equal to 1 (first call) or 0.
0006 %   o uniform          [integer]    scalar equal to 1 (first call) or 0.
0007 %
0008 % OUTPUTS
0009 %   o pdraw            [double]     1*npar vector, draws from the joint prior density.
0010 %
0011 %
0012 % SPECIAL REQUIREMENTS
0013 %   none
0014 %
0015 % NOTE 1. Input arguments 1 an 2 are only needed for initialization.
0016 % NOTE 2. A given draw from the joint prior distribution does not satisfy BK conditions a priori.
0017 
0018 % Copyright (C) 2006-2010 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 persistent p6 p7 p3 p4 lb ub
0036 persistent uniform_index gaussian_index gamma_index beta_index inverse_gamma_1_index inverse_gamma_2_index
0037 persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws
0038 
0039 
0040 if nargin>0 && init
0041     p6 = evalin('base', 'bayestopt_.p6');
0042     p7 = evalin('base', 'bayestopt_.p7');
0043     p3 = evalin('base', 'bayestopt_.p3');
0044     p4 = evalin('base', 'bayestopt_.p4');
0045     lb = evalin('base', 'bayestopt_.lb');
0046     ub = evalin('base', 'bayestopt_.ub');
0047     number_of_estimated_parameters = length(p6);
0048     if nargin>1 && uniform
0049         prior_shape = repmat(5,number_of_estimated_parameters,1);
0050     else
0051         prior_shape = evalin('base', 'bayestopt_.pshape');
0052     end
0053     beta_index = find(prior_shape==1);
0054     if isempty(beta_index)
0055         beta_draws = 0;
0056     else
0057         beta_draws = 1;
0058     end
0059     gamma_index = find(prior_shape==2);
0060     if isempty(gamma_index)
0061         gamma_draws = 0;
0062     else
0063         gamma_draws = 1;
0064     end
0065     gaussian_index = find(prior_shape==3);
0066     if isempty(gaussian_index)
0067         gaussian_draws = 0;
0068     else
0069         gaussian_draws = 1;
0070     end
0071     inverse_gamma_1_index = find(prior_shape==4);
0072     if isempty(inverse_gamma_1_index)
0073         inverse_gamma_1_draws = 0;
0074     else
0075         inverse_gamma_1_draws = 1;
0076     end
0077     uniform_index = find(prior_shape==5);
0078     if isempty(uniform_index)
0079         uniform_draws = 0;
0080     else
0081         uniform_draws = 1;
0082     end
0083     inverse_gamma_2_index = find(prior_shape==6);
0084     if isempty(inverse_gamma_2_index)
0085         inverse_gamma_2_draws = 0;
0086     else
0087         inverse_gamma_2_draws = 1;
0088     end
0089     pdraw = zeros(number_of_estimated_parameters,1);
0090     return
0091 end
0092 
0093 if uniform_draws
0094     pdraw(uniform_index) = rand(length(uniform_index),1).*(p4(uniform_index)-p3(uniform_index)) + p3(uniform_index);  
0095 end
0096 
0097 if gaussian_draws
0098     pdraw(gaussian_index) = randn(length(gaussian_index),1).*p7(gaussian_index) + p6(gaussian_index);
0099     out_of_bound = find( (pdraw(gaussian_index)'>ub(gaussian_index)) | (pdraw(gaussian_index)'<lb(gaussian_index)));
0100     while ~isempty(out_of_bound),
0101         pdraw(gaussian_index(out_of_bound)) = randn(length(gaussian_index(out_of_bound)),1).*p7(gaussian_index(out_of_bound)) + p6(gaussian_index(out_of_bound));
0102         out_of_bound = find( (pdraw(gaussian_index)'>ub(gaussian_index)) | (pdraw(gaussian_index)'<lb(gaussian_index)));
0103     end
0104 end
0105 
0106 if gamma_draws
0107     pdraw(gamma_index) = gamrnd(p6(gamma_index),p7(gamma_index))+p3(gamma_index);
0108     out_of_bound = find( (pdraw(gamma_index)'>ub(gamma_index)) | (pdraw(gamma_index)'<lb(gamma_index)));
0109     while ~isempty(out_of_bound),
0110         pdraw(gamma_index(out_of_bound)) = gamrnd(p6(gamma_index(out_of_bound)),p7(gamma_index(out_of_bound)))+p3(gamma_index(out_of_bound));
0111         out_of_bound = find( (pdraw(gamma_index)'>ub(gamma_index)) | (pdraw(gamma_index)'<lb(gamma_index)));
0112     end
0113 end
0114 
0115 if beta_draws
0116     pdraw(beta_index) = (p4(beta_index)-p3(beta_index)).*betarnd(p6(beta_index),p7(beta_index))+p3(beta_index);
0117     out_of_bound = find( (pdraw(beta_index)'>ub(beta_index)) | (pdraw(beta_index)'<lb(beta_index)));
0118     while ~isempty(out_of_bound),
0119         pdraw(beta_index(out_of_bound)) = (p4(beta_index(out_of_bound))-p3(beta_index(out_of_bound))).*betarnd(p6(beta_index(out_of_bound)),p7(beta_index(out_of_bound)))+p3(beta_index(out_of_bound));
0120         out_of_bound = find( (pdraw(beta_index)'>ub(beta_index)) | (pdraw(beta_index)'<lb(beta_index)));
0121     end
0122 end
0123 
0124 if inverse_gamma_1_draws
0125     pdraw(inverse_gamma_1_index) = ...
0126         sqrt(1./gamrnd(p7(inverse_gamma_1_index)/2,2./p6(inverse_gamma_1_index)))+p3(inverse_gamma_1_index);
0127     out_of_bound = find( (pdraw(inverse_gamma_1_index)'>ub(inverse_gamma_1_index)) | (pdraw(inverse_gamma_1_index)'<lb(inverse_gamma_1_index)));
0128     while ~isempty(out_of_bound),
0129         pdraw(inverse_gamma_1_index(out_of_bound)) = ...
0130             sqrt(1./gamrnd(p7(inverse_gamma_1_index(out_of_bound))/2,2./p6(inverse_gamma_1_index(out_of_bound))))+p3(inverse_gamma_1_index(out_of_bound));
0131         out_of_bound = find( (pdraw(inverse_gamma_1_index)'>ub(inverse_gamma_1_index)) | (pdraw(inverse_gamma_1_index)'<lb(inverse_gamma_1_index)));
0132     end
0133 end
0134 
0135 if inverse_gamma_2_draws
0136     pdraw(inverse_gamma_2_index) = ...
0137         1./gamrnd(p7(inverse_gamma_2_index)/2,2./p6(inverse_gamma_2_index))+p3(inverse_gamma_2_index);
0138     out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'<lb(inverse_gamma_2_index)));
0139     while ~isempty(out_of_bound),
0140         pdraw(inverse_gamma_2_index(out_of_bound)) = ...
0141             sqrt(1./gamrnd(p7(inverse_gamma_2_index(out_of_bound))/2,2./p6(inverse_gamma_2_index(out_of_bound))))+p3(inverse_gamma_2_index(out_of_bound));
0142         out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'<lb(inverse_gamma_2_index)));
0143     end
0144 end