set_prior

 function [xparam1,estim_params_,bayestopt_,lb,ub]=set_prior(estim_params_)
 sets prior distributions

 INPUTS
    o estim_params_    [structure] characterizing parameters to be estimated.
    o M_               [structure] characterizing the model. 
    o options_         [structure] 
    
 OUTPUTS
    o xparam1          [double]    vector of parameters to be estimated (initial values)
    o estim_params_    [structure] characterizing parameters to be estimated
    o bayestopt_       [structure] characterizing priors
    o lb               [double]    vector of lower bounds for the estimated parameters. 
    o ub               [double]    vector of upper bounds for the estimated parameters.
    o M_               [structure] characterizing the model.
    
 SPECIAL REQUIREMENTS
    None

0001 function [xparam1, estim_params_, bayestopt_, lb, ub, M_]=set_prior(estim_params_, M_, options_)
0002 % function [xparam1,estim_params_,bayestopt_,lb,ub]=set_prior(estim_params_)
0003 % sets prior distributions
0004 %
0005 % INPUTS
0006 %    o estim_params_    [structure] characterizing parameters to be estimated.
0007 %    o M_               [structure] characterizing the model.
0008 %    o options_         [structure]
0009 %
0010 % OUTPUTS
0011 %    o xparam1          [double]    vector of parameters to be estimated (initial values)
0012 %    o estim_params_    [structure] characterizing parameters to be estimated
0013 %    o bayestopt_       [structure] characterizing priors
0014 %    o lb               [double]    vector of lower bounds for the estimated parameters.
0015 %    o ub               [double]    vector of upper bounds for the estimated parameters.
0016 %    o M_               [structure] characterizing the model.
0017 %
0018 % SPECIAL REQUIREMENTS
0019 %    None
0020 
0021 % Copyright (C) 2003-2011 Dynare Team
0022 %
0023 % This file is part of Dynare.
0024 %
0025 % Dynare is free software: you can redistribute it and/or modify
0026 % it under the terms of the GNU General Public License as published by
0027 % the Free Software Foundation, either version 3 of the License, or
0028 % (at your option) any later version.
0029 %
0030 % Dynare is distributed in the hope that it will be useful,
0031 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0032 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0033 % GNU General Public License for more details.
0034 %
0035 % You should have received a copy of the GNU General Public License
0036 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0037 
0038 nvx = size(estim_params_.var_exo,1);
0039 nvn = size(estim_params_.var_endo,1);
0040 ncx = size(estim_params_.corrx,1);
0041 ncn = size(estim_params_.corrn,1);
0042 np = size(estim_params_.param_vals,1);
0043 
0044 estim_params_.nvx = nvx;
0045 estim_params_.nvn = nvn;
0046 estim_params_.ncx = ncx;
0047 estim_params_.ncn = ncn;
0048 estim_params_.np = np;
0049 
0050 xparam1 = [];
0051 ub = [];
0052 lb = [];
0053 bayestopt_.pshape = [];
0054 bayestopt_.p1 = []; % prior mean
0055 bayestopt_.p2 = []; % prior standard deviation
0056 bayestopt_.p3 = []; % lower bound
0057 bayestopt_.p4 = []; % upper bound
0058 bayestopt_.p5 = zeros(nvx+nvn+ncx+ncn+np,1); % prior mode
0059 bayestopt_.p6 = []; % first hyper-parameter (\alpha for the BETA and GAMMA distributions, s for the INVERSE GAMMAs, expectation for the GAUSSIAN distribution, lower bound for the UNIFORM distribution).
0060 bayestopt_.p7 = []; % second hyper-parameter (\beta for the BETA and GAMMA distributions, \nu for the INVERSE GAMMAs, standard deviation for the GAUSSIAN distribution, upper bound for the UNIFORM distribution).
0061 
0062 bayestopt_.jscale = [];
0063 bayestopt_.name = [];
0064 if nvx
0065     xparam1 = estim_params_.var_exo(:,2);
0066     ub = estim_params_.var_exo(:,4); 
0067     lb = estim_params_.var_exo(:,3); 
0068     bayestopt_.pshape =  estim_params_.var_exo(:,5);
0069     bayestopt_.p1 =  estim_params_.var_exo(:,6);
0070     bayestopt_.p2 =  estim_params_.var_exo(:,7);
0071     bayestopt_.p3 =  estim_params_.var_exo(:,8);
0072     bayestopt_.p4 =  estim_params_.var_exo(:,9);
0073     bayestopt_.jscale =  estim_params_.var_exo(:,10);
0074     bayestopt_.name = cellstr(M_.exo_names(estim_params_.var_exo(:,1),:));
0075 end
0076 if nvn
0077     if isequal(M_.H,0)
0078         nvarobs = size(options_.varobs,1);
0079         M_.H = zeros(nvarobs,nvarobs);
0080     end
0081     for i=1:nvn
0082         obsi_ = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact');
0083         if isempty(obsi_)
0084             error(['The variable ' deblank(M_.endo_names(estim_params_.var_endo(i,1),:)) ' has to be declared as observable since you assume a measurement error on it.'])
0085         end
0086         estim_params_.var_endo(i,1) = obsi_;
0087     end
0088     xparam1 = [xparam1; estim_params_.var_endo(:,2)];
0089     ub = [ub; estim_params_.var_endo(:,4)]; 
0090     lb = [lb; estim_params_.var_endo(:,3)]; 
0091     bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.var_endo(:,5)];
0092     bayestopt_.p1 = [ bayestopt_.p1; estim_params_.var_endo(:,6)];
0093     bayestopt_.p2 = [ bayestopt_.p2; estim_params_.var_endo(:,7)];
0094     bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)];
0095     bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)];
0096     bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
0097     if isempty(bayestopt_.name)
0098         bayestopt_.name = cellstr(char(options_.varobs(estim_params_.var_endo(:,1),:)));
0099     else
0100         bayestopt_.name = cellstr(char(char(bayestopt_.name), options_.varobs(estim_params_.var_endo(:,1),:)));
0101     end
0102 end
0103 if ncx
0104     xparam1 = [xparam1; estim_params_.corrx(:,3)];
0105     ub = [ub; max(min(estim_params_.corrx(:,5),1),-1)];
0106     lb = [lb; max(min(estim_params_.corrx(:,4),1),-1)];
0107     bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrx(:,6)];
0108     bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrx(:,7)];
0109     bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrx(:,8)];
0110     bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrx(:,9)];
0111     bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrx(:,10)];
0112     bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrx(:,11)];
0113     if isempty(bayestopt_.name)
0114         bayestopt_.name = cellstr(char(char(strcat(cellstr(M_.exo_names(estim_params_.corrx(:,1),:)), ...
0115                                                    ',' , cellstr(M_.exo_names(estim_params_.corrx(:,2),:))))));
0116     else
0117         bayestopt_.name = cellstr(char(char(bayestopt_.name), char(strcat(cellstr(M_.exo_names(estim_params_.corrx(:,1),:)), ...
0118                                                           ',' , cellstr(M_.exo_names(estim_params_.corrx(:,2),:))))));
0119     end
0120 end
0121 if ncn
0122     if isequal(M_.H,0)
0123         nvarobs = size(options_.varobs,1);
0124         M_.H = zeros(nvarobs,nvarobs);
0125     end
0126     xparam1 = [xparam1; estim_params_.corrn(:,3)];
0127     ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
0128     lb = [lb; max(min(estim_params_.corrn(:,4),1),-1)];
0129     bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.corrn(:,6)];
0130     bayestopt_.p1 = [ bayestopt_.p1; estim_params_.corrn(:,7)];
0131     bayestopt_.p2 = [ bayestopt_.p2; estim_params_.corrn(:,8)];
0132     bayestopt_.p3 = [ bayestopt_.p3; estim_params_.corrn(:,9)];
0133     bayestopt_.p4 = [ bayestopt_.p4; estim_params_.corrn(:,10)];
0134     bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.corrn(:,11)];
0135     if isempty(bayestopt_.name)
0136         bayestopt_.name = cellstr(char(char(strcat(cellstr(M_.endo_names(estim_params_.corrn(:,1),:)),...
0137                                                    ',' ,  cellstr(M_.endo_names(estim_params_.corrn(:,2),:))))));
0138     else
0139         bayestopt_.name = cellstr(char(char(bayestopt_.name), char(strcat(cellstr(M_.endo_names(estim_params_.corrn(:,1),:)),...
0140                                                           ',' ,  cellstr(M_.endo_names(estim_params_.corrn(:,2),:))))));
0141     end
0142 end
0143 if np
0144     xparam1 = [xparam1; estim_params_.param_vals(:,2)];
0145     ub = [ub; estim_params_.param_vals(:,4)];
0146     lb = [lb; estim_params_.param_vals(:,3)];
0147     bayestopt_.pshape = [ bayestopt_.pshape; estim_params_.param_vals(:,5)];
0148     bayestopt_.p1 = [ bayestopt_.p1; estim_params_.param_vals(:,6)];
0149     bayestopt_.p2 = [ bayestopt_.p2; estim_params_.param_vals(:,7)];
0150     bayestopt_.p3 = [ bayestopt_.p3; estim_params_.param_vals(:,8)];
0151     bayestopt_.p4 = [ bayestopt_.p4; estim_params_.param_vals(:,9)];
0152     bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.param_vals(:,10)];
0153     if isempty(bayestopt_.name)
0154         bayestopt_.name = cellstr(char(M_.param_names(estim_params_.param_vals(:,1),:)));
0155     else
0156         bayestopt_.name = cellstr(char(char(bayestopt_.name),M_.param_names(estim_params_.param_vals(:,1),:)));
0157     end
0158 end
0159 
0160 bayestopt_.ub = ub;
0161 bayestopt_.lb = lb;
0162 
0163 bayestopt_.p6 = NaN(size(bayestopt_.p1)) ;
0164 bayestopt_.p7 = bayestopt_.p6 ;
0165 
0166 % generalized location parameters by default for beta distribution
0167 k = find(bayestopt_.pshape == 1);
0168 k1 = find(isnan(bayestopt_.p3(k)));
0169 bayestopt_.p3(k(k1)) = zeros(length(k1),1);
0170 k1 = find(isnan(bayestopt_.p4(k)));
0171 bayestopt_.p4(k(k1)) = ones(length(k1),1);
0172 for i=1:length(k)
0173     if (bayestopt_.p1(k(i))<bayestopt_.p3(k(i))) || (bayestopt_.p1(k(i))>bayestopt_.p4(k(i)))
0174         error(['The prior mean of ' bayestopt_.name{k(i)} ' has to be between the lower (' num2str(bayestopt_.p3(k(i))) ') and upper (' num2str(bayestopt_.p4(k(i))) ') bounds of the beta prior density!']);
0175     end
0176     mu = (bayestopt_.p1(k(i))-bayestopt_.p3(k(i)))/(bayestopt_.p4(k(i))-bayestopt_.p3(k(i)));
0177     stdd = bayestopt_.p2(k(i))/(bayestopt_.p4(k(i))-bayestopt_.p3(k(i)));
0178     bayestopt_.p6(k(i)) = (1-mu)*mu^2/stdd^2 - mu ;
0179     bayestopt_.p7(k(i)) = bayestopt_.p6(k(i))*(1/mu-1) ;
0180     m = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) , bayestopt_.p4(k(i)) ],1);
0181     if length(m)==1
0182         bayestopt_.p5(k(i)) = m;
0183     else
0184         disp(['Prior distribution for parameter ' bayestopt_.name{k(i)}  ' has two modes!'])
0185         bayestopt_.p5(k(i)) = bayestopt_.p1(k(i)) ; 
0186     end
0187 end
0188 
0189 % generalized location parameter by default for gamma distribution
0190 k =  find(bayestopt_.pshape == 2);
0191 k1 = find(isnan(bayestopt_.p3(k)));
0192 k2 = find(isnan(bayestopt_.p4(k)));
0193 bayestopt_.p3(k(k1)) = zeros(length(k1),1);
0194 bayestopt_.p4(k(k2)) = Inf(length(k2),1);
0195 for i=1:length(k)
0196     if isinf(bayestopt_.p2(k(i)))
0197         error(['Infinite prior standard deviation for parameter ' bayestopt_.name{k(i)}  ' is not allowed (Gamma prior)!'])
0198     end
0199     mu = bayestopt_.p1(k(i))-bayestopt_.p3(k(i));
0200     bayestopt_.p7(k(i)) = bayestopt_.p2(k(i))^2/mu ;
0201     bayestopt_.p6(k(i)) = mu/bayestopt_.p7(k(i)) ;  
0202     bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 2) ;
0203 end
0204 
0205 % truncation parameters by default for normal distribution
0206 k  = find(bayestopt_.pshape == 3);
0207 k1 = find(isnan(bayestopt_.p3(k)));
0208 bayestopt_.p3(k(k1)) = -Inf*ones(length(k1),1);
0209 k1 = find(isnan(bayestopt_.p4(k)));
0210 bayestopt_.p4(k(k1)) = Inf*ones(length(k1),1);
0211 for i=1:length(k)
0212     bayestopt_.p6(k(i)) = bayestopt_.p1(k(i)) ; 
0213     bayestopt_.p7(k(i)) = bayestopt_.p2(k(i)) ;
0214     bayestopt_.p5(k(i)) = bayestopt_.p1(k(i)) ;
0215 end
0216 
0217 % inverse gamma distribution (type 1)
0218 k = find(bayestopt_.pshape == 4);
0219 k1 = find(isnan(bayestopt_.p3(k)));
0220 k2 = find(isnan(bayestopt_.p4(k)));
0221 bayestopt_.p3(k(k1)) = zeros(length(k1),1);
0222 bayestopt_.p4(k(k2)) = Inf(length(k2),1);
0223 for i=1:length(k)
0224     [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
0225         inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),1,0) ;
0226     bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 4) ;
0227 end  
0228 
0229 % uniform distribution
0230 k = find(bayestopt_.pshape == 5);
0231 for i=1:length(k)
0232     [bayestopt_.p1(k(i)),bayestopt_.p2(k(i)),bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
0233         uniform_specification(bayestopt_.p1(k(i)),bayestopt_.p2(k(i)),bayestopt_.p3(k(i)),bayestopt_.p4(k(i)));
0234     bayestopt_.p3(k(i)) = bayestopt_.p6(k(i)) ;
0235     bayestopt_.p4(k(i)) = bayestopt_.p7(k(i)) ;
0236     bayestopt_.p5(k(i)) = NaN ;
0237 end
0238 
0239 % inverse gamma distribution (type 2)
0240 k = find(bayestopt_.pshape == 6);
0241 k1 = find(isnan(bayestopt_.p3(k)));
0242 k2 = find(isnan(bayestopt_.p4(k)));
0243 bayestopt_.p3(k(k1)) = zeros(length(k1),1);
0244 bayestopt_.p4(k(k2)) = Inf(length(k2),1);
0245 for i=1:length(k)
0246     [bayestopt_.p6(k(i)),bayestopt_.p7(k(i))] = ...
0247         inverse_gamma_specification(bayestopt_.p1(k(i))-bayestopt_.p3(k(i)),bayestopt_.p2(k(i)),2,0);
0248     bayestopt_.p5(k(i)) = compute_prior_mode([ bayestopt_.p6(k(i)) , bayestopt_.p7(k(i)) , bayestopt_.p3(k(i)) ], 6) ;
0249 end
0250 
0251 k = find(isnan(xparam1));
0252 if ~isempty(k)
0253     xparam1(k) = bayestopt_.p1(k);
0254 end
0255 
0256 if options_.initialize_estimated_parameters_with_the_prior_mode
0257     xparam1 = bayestopt_.p5;
0258     k = find(isnan(xparam1));% Because the uniform density do not have a mode!
0259     if ~isempty(k)
0260         xparam1(k) = bayestopt_.p1(k);
0261     end
0262     xparam1 = transpose(xparam1);
0263 end 
0264 
0265 % I create subfolder M_.dname/prior if needed.
0266 CheckPath('prior',M_.dname);
0267 
0268 % I save the prior definition if the prior has changed.
0269 if exist([ M_.dname '/prior/definition.mat'])
0270     old = load([M_.dname '/prior/definition.mat'],'bayestopt_');
0271     prior_has_changed = 0;
0272     if length(bayestopt_.p1)==length(old.bayestopt_.p1)
0273         if any(bayestopt_.p1-old.bayestopt_.p1)
0274             prior_has_changed = 1;
0275         elseif any(bayestopt_.p2-old.bayestopt_.p2)
0276             prior_has_changed = 1;
0277         elseif any(bayestopt_.p3-old.bayestopt_.p3)
0278             prior_has_changed = 1;
0279         elseif any(bayestopt_.p4-old.bayestopt_.p4)
0280             prior_has_changed = 1;
0281         elseif any(bayestopt_.p5-old.bayestopt_.p5(:))
0282             prior_has_changed = 1;
0283         elseif any(bayestopt_.p6-old.bayestopt_.p6)
0284             prior_has_changed = 1;
0285         elseif any(bayestopt_.p7-old.bayestopt_.p7)
0286             prior_has_changed = 1;
0287         end
0288     else
0289         prior_has_changed = 1;
0290     end
0291     if prior_has_changed
0292         delete([M_.dname '/prior/definition.mat']);
0293         save([M_.dname '/prior/definition.mat'],'bayestopt_');
0294     end
0295 else
0296     save([M_.dname '/prior/definition.mat'],'bayestopt_');
0297 end
0298 
0299 % initialize persistent variables in priordens()
0300 priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, ...
0301           bayestopt_.p3,bayestopt_.p4,1);
0302 
0303 % Put bayestopt_ in matlab's workspace
0304 assignin('base','bayestopt_',bayestopt_);