1function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4, initialization) % --*-- Unitary tests --*-- 2% Computes a prior density for the structural parameters of DSGE models 3% 4% INPUTS 5% x [double] vector with n elements. 6% pshape [integer] vector with n elements (bayestopt_.pshape). 7% p6: [double] vector with n elements, first parameter of the prior distribution (bayestopt_.p6). 8% p7: [double] vector with n elements, second parameter of the prior distribution (bayestopt_.p7). 9% p3: [double] vector with n elements, lower bounds of the untruncated standard or generalized distribution 10% p4: [double] vector with n elements, upper bound of the untruncated standard or generalized distribution 11% initialization [integer] if 1: initialize persistent variables 12% 13% OUTPUTS 14% logged_prior_density [double] scalar, log of the prior density evaluated at x. 15% info [double] error code for index of Inf-prior parameter 16% 17 18% Copyright (C) 2003-2017 Dynare Team 19% 20% This file is part of Dynare. 21% 22% Dynare is free software: you can redistribute it and/or modify 23% it under the terms of the GNU General Public License as published by 24% the Free Software Foundation, either version 3 of the License, or 25% (at your option) any later version. 26% 27% Dynare is distributed in the hope that it will be useful, 28% but WITHOUT ANY WARRANTY; without even the implied warranty of 29% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 30% GNU General Public License for more details. 31% 32% You should have received a copy of the GNU General Public License 33% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 34 35persistent id1 id2 id3 id4 id5 id6 id8 36persistent tt1 tt2 tt3 tt4 tt5 tt6 tt8 37 38info=0; 39 40if nargin > 6 && initialization 41 % Beta indices. 42 tt1 = true; 43 id1 = find(pshape==1); 44 if isempty(id1) 45 tt1 = false; 46 end 47 % Gamma indices. 48 tt2 = true; 49 id2 = find(pshape==2); 50 if isempty(id2) 51 tt2 = false; 52 end 53 % Gaussian indices. 54 tt3 = true; 55 id3 = find(pshape==3); 56 if isempty(id3) 57 tt3 = false; 58 end 59 % Inverse-Gamma-1 indices. 60 tt4 = true; 61 id4 = find(pshape==4); 62 if isempty(id4) 63 tt4 = false; 64 end 65 % Uniform indices. 66 tt5 = true; 67 id5 = find(pshape==5); 68 if isempty(id5) 69 tt5 = false; 70 end 71 % Inverse-Gamma-2 indices. 72 tt6 = true; 73 id6 = find(pshape==6); 74 if isempty(id6) 75 tt6 = false; 76 end 77 % Weibull indices. 78 tt8 = true; 79 id8 = find(pshape==8); 80 if isempty(id8) 81 tt8 = false; 82 end 83end 84 85logged_prior_density = 0.0; 86dlprior = zeros(1,length(x)); 87d2lprior = dlprior; 88 89if tt1 90 logged_prior_density = logged_prior_density + sum(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1))) ; 91 if isinf(logged_prior_density) 92 if nargout ==4 93 info=id1(isinf(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1)))); 94 end 95 return 96 end 97 if nargout == 2 98 [tmp, dlprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1)); 99 elseif nargout == 3 100 [tmp, dlprior(id1), d2lprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1)); 101 end 102end 103 104if tt2 105 logged_prior_density = logged_prior_density + sum(lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2))) ; 106 if isinf(logged_prior_density) 107 if nargout ==4 108 info=id2(isinf(lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2)))); 109 end 110 return 111 end 112 if nargout == 2 113 [tmp, dlprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2)); 114 elseif nargout == 3 115 [tmp, dlprior(id2), d2lprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2)); 116 end 117end 118 119if tt3 120 logged_prior_density = logged_prior_density + sum(lpdfnorm(x(id3),p6(id3),p7(id3))) ; 121 if nargout == 2 122 [tmp, dlprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3)); 123 elseif nargout == 3 124 [tmp, dlprior(id3), d2lprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3)); 125 end 126end 127 128if tt4 129 logged_prior_density = logged_prior_density + sum(lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4))) ; 130 if isinf(logged_prior_density) 131 if nargout ==4 132 info=id4(isinf(lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4)))); 133 end 134 return 135 end 136 if nargout == 2 137 [tmp, dlprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4)); 138 elseif nargout == 3 139 [tmp, dlprior(id4), d2lprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4)); 140 end 141end 142 143if tt5 144 if any(x(id5)-p3(id5)<0) || any(x(id5)-p4(id5)>0) 145 logged_prior_density = -Inf ; 146 if nargout ==4 147 info=id5((x(id5)-p3(id5)<0) || (x(id5)-p4(id5)>0)); 148 end 149 return 150 end 151 logged_prior_density = logged_prior_density + sum(log(1./(p4(id5)-p3(id5)))) ; 152 if nargout >1 153 dlprior(id5)=zeros(length(id5),1); 154 end 155 if nargout == 3 156 d2lprior(id5)=zeros(length(id5),1); 157 end 158end 159 160if tt6 161 logged_prior_density = logged_prior_density + sum(lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6))) ; 162 if isinf(logged_prior_density) 163 if nargout ==4 164 info=id6(isinf(lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6)))); 165 end 166 return 167 end 168 if nargout == 2 169 [tmp, dlprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6)); 170 elseif nargout == 3 171 [tmp, dlprior(id6), d2lprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6)); 172 end 173end 174 175if tt8 176 logged_prior_density = logged_prior_density + sum(lpdfgweibull(x(id8),p6(id8),p7(id8))); 177 if isinf(logged_prior_density) 178 if nargout ==4 179 info=id8(isinf(log(lpdfgweibull(x(id8),p6(id8),p7(id8))))); 180 end 181 return 182 end 183 if nargout==2 184 [tmp, dlprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8)); 185 elseif nargout==3 186 [tmp, dlprior(id8), d2lprior(id8)] = lpdfgweibull(x(id8),p6(id8),p7(id8)); 187 end 188end 189 190if nargout==3 191 d2lprior = diag(d2lprior); 192end 193 194%@test:1 195%$ % Fill global structures with required fields... 196%$ prior_trunc = 1e-10; 197%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape 198%$ p1 = .4*ones(14,1); % Prior mean 199%$ p2 = .2*ones(14,1); % Prior std. 200%$ p3 = NaN(14,1); 201%$ p4 = NaN(14,1); 202%$ p5 = NaN(14,1); 203%$ p6 = NaN(14,1); 204%$ p7 = NaN(14,1); 205%$ 206%$ for i=1:14 207%$ switch p0(i) 208%$ case 1 209%$ % Beta distribution 210%$ p3(i) = 0; 211%$ p4(i) = 1; 212%$ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i)); 213%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 1); 214%$ case 2 215%$ % Gamma distribution 216%$ p3(i) = 0; 217%$ p4(i) = Inf; 218%$ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i)); 219%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 2); 220%$ case 3 221%$ % Normal distribution 222%$ p3(i) = -Inf; 223%$ p4(i) = Inf; 224%$ p6(i) = p1(i); 225%$ p7(i) = p2(i); 226%$ p5(i) = p1(i); 227%$ case 4 228%$ % Inverse Gamma (type I) distribution 229%$ p3(i) = 0; 230%$ p4(i) = Inf; 231%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false); 232%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 4); 233%$ case 5 234%$ % Uniform distribution 235%$ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i)); 236%$ p3(i) = p6(i); 237%$ p4(i) = p7(i); 238%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 5); 239%$ case 6 240%$ % Inverse Gamma (type II) distribution 241%$ p3(i) = 0; 242%$ p4(i) = Inf; 243%$ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false); 244%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 6); 245%$ case 8 246%$ % Weibull distribution 247%$ p3(i) = 0; 248%$ p4(i) = Inf; 249%$ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i)); 250%$ p5(i) = compute_prior_mode([p6(i) p7(i)], 8); 251%$ otherwise 252%$ error('This density is not implemented!') 253%$ end 254%$ end 255%$ 256%$ % Call the tested routine 257%$ try 258%$ % Initialization of priordens. 259%$ lpdstar = priordens(p5, p0, p6, p7, p3, p4, true); 260%$ % Do simulations in a loop and estimate recursively the mean and teh variance. 261%$ LPD = NaN(10000,1); 262%$ for i = 1:10000 263%$ draw = p5+randn(size(p5))*.02; 264%$ LPD(i) = priordens(p5, p0, p6, p7, p3, p4); 265%$ end 266%$ t(1) = true; 267%$ catch 268%$ t(1) = false; 269%$ end 270%$ 271%$ if t(1) 272%$ t(2) = all(LPD<=lpdstar); 273%$ end 274%$ T = all(t); 275%@eof:1