prior_draw_gsa

 Draws from the prior distributions 
 Adapted by M. Ratto from prior_draw (of DYNARE, copyright M. Juillard), 
 for use with Sensitivity Toolbox for DYNARE
 
 
 INPUTS
   o init           [integer]  scalar equal to 1 (first call) or 0.
   o rdraw          
    
 OUTPUTS 
   o pdraw          [double]   draw from the joint prior density. 

 ALGORITHM 
   ...       

 SPECIAL REQUIREMENTS
   MATLAB Statistics Toolbox
  
 Written by Marco Ratto
 Joint Research Centre, The European Commission,
 (http://eemc.jrc.ec.europa.eu/),
 marco.ratto@jrc.it 

 Reference:
 M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.

0001 function pdraw = prior_draw_gsa(init,rdraw)
0002 % Draws from the prior distributions
0003 % Adapted by M. Ratto from prior_draw (of DYNARE, copyright M. Juillard),
0004 % for use with Sensitivity Toolbox for DYNARE
0005 %
0006 %
0007 % INPUTS
0008 %   o init           [integer]  scalar equal to 1 (first call) or 0.
0009 %   o rdraw
0010 %
0011 % OUTPUTS
0012 %   o pdraw          [double]   draw from the joint prior density.
0013 %
0014 % ALGORITHM
0015 %   ...
0016 %
0017 % SPECIAL REQUIREMENTS
0018 %   MATLAB Statistics Toolbox
0019 %
0020 % Written by Marco Ratto
0021 % Joint Research Centre, The European Commission,
0022 % (http://eemc.jrc.ec.europa.eu/),
0023 % marco.ratto@jrc.it
0024 %
0025 % Reference:
0026 % M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
0027 
0028 % Copyright (C) 2012 Dynare Team
0029 %
0030 % This file is part of Dynare.
0031 %
0032 % Dynare is free software: you can redistribute it and/or modify
0033 % it under the terms of the GNU General Public License as published by
0034 % the Free Software Foundation, either version 3 of the License, or
0035 % (at your option) any later version.
0036 %
0037 % Dynare is distributed in the hope that it will be useful,
0038 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0040 % GNU General Public License for more details.
0041 %
0042 % You should have received a copy of the GNU General Public License
0043 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0044 
0045 % global M_ options_ estim_params_  bayestopt_
0046 global bayestopt_
0047 persistent npar pshape p6 p7 p3 p4 lbcum ubcum
0048   
0049 if init
0050     pshape = bayestopt_.pshape;
0051     p6 = bayestopt_.p6;
0052     p7 = bayestopt_.p7;
0053     p3 = bayestopt_.p3;
0054     p4 = bayestopt_.p4;
0055     npar = length(p6);
0056     pdraw = zeros(npar,1);
0057     lbcum = zeros(npar,1);
0058     ubcum = ones(npar,1);
0059     
0060     % set bounds for cumulative probabilities
0061     for i = 1:npar
0062       switch pshape(i)
0063         case 5% Uniform prior.
0064           p4(i) = min(p4(i),bayestopt_.ub(i));
0065           p3(i) = max(p3(i),bayestopt_.lb(i));
0066         case 3% Gaussian prior.
0067           lbcum(i) = 0.5 * erfc(-(bayestopt_.lb(i)-p6(i))/p7(i) ./ sqrt(2));;
0068           ubcum(i) = 0.5 * erfc(-(bayestopt_.ub(i)-p6(i))/p7(i) ./ sqrt(2));;
0069         case 2% Gamma prior.
0070           lbcum(i) = gamcdf(bayestopt_.lb(i)-p3(i),p6(i),p7(i));
0071           ubcum(i) = gamcdf(bayestopt_.ub(i)-p3(i),p6(i),p7(i));
0072         case 1% Beta distribution (TODO: generalized beta distribution)
0073           lbcum(i) = betainc((bayestopt_.lb(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
0074           ubcum(i) = betainc((bayestopt_.ub(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
0075         case 4% INV-GAMMA1 distribution
0076           % TO BE CHECKED
0077           lbcum(i) = gamcdf(1/(bayestopt_.ub(i)-p3(i))^2,p7(i)/2,2/p6(i));
0078           ubcum(i) = gamcdf(1/(bayestopt_.lb(i)-p3(i))^2,p7(i)/2,2/p6(i));
0079         case 6% INV-GAMMA2 distribution
0080           % TO BE CHECKED
0081           lbcum(i) = gamcdf(1/(bayestopt_.ub(i)-p3(i)),p7(i)/2,2/p6(i));
0082           ubcum(i) = gamcdf(1/(bayestopt_.lb(i)-p3(i)),p7(i)/2,2/p6(i));
0083         otherwise
0084           % Nothing to do here.
0085       end
0086     end
0087     return
0088 end
0089 
0090 
0091 for i = 1:npar
0092     rdraw(:,i) = rdraw(:,i).*(ubcum(i)-lbcum(i))+lbcum(i);
0093     switch pshape(i)
0094       case 5% Uniform prior.
0095         pdraw(:,i) = rdraw(:,i)*(p4(i)-p3(i)) + p3(i);
0096       case 3% Gaussian prior.
0097         pdraw(:,i) = norminv(rdraw(:,i),p6(i),p7(i));
0098       case 2% Gamma prior.
0099         pdraw(:,i) = gaminv(rdraw(:,i),p6(i),p7(i))+p3(i);
0100       case 1% Beta distribution (TODO: generalized beta distribution)
0101         pdraw(:,i) = betainv(rdraw(:,i),p6(i),p7(i))*(p4(i)-p3(i))+p3(i);
0102       case 4% INV-GAMMA1 distribution
0103         % TO BE CHECKED
0104         pdraw(:,i) =  sqrt(1./gaminv(rdraw(:,i),p7(i)/2,2/p6(i)))+p3(i);
0105       case 6% INV-GAMMA2 distribution
0106         % TO BE CHECKED
0107         pdraw(:,i) =  1./gaminv(rdraw(:,i),p7(i)/2,2/p6(i))+p3(i);
0108       otherwise
0109         % Nothing to do here.
0110     end
0111 end
0112 
0113