1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
function pdraw = prior_draw(init,uniform)
% 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.
% Copyright (C) 2006-2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
persistent p6 p7 p3 p4 lb ub
persistent uniform_index gaussian_index gamma_index beta_index inverse_gamma_1_index inverse_gamma_2_index
persistent uniform_draws gaussian_draws gamma_draws beta_draws inverse_gamma_1_draws inverse_gamma_2_draws
if nargin>0 && init
p6 = evalin('base', 'bayestopt_.p6');
p7 = evalin('base', 'bayestopt_.p7');
p3 = evalin('base', 'bayestopt_.p3');
p4 = evalin('base', 'bayestopt_.p4');
lb = evalin('base', 'bayestopt_.lb');
ub = evalin('base', 'bayestopt_.ub');
number_of_estimated_parameters = length(p6);
if nargin>1 && uniform
prior_shape = repmat(5,number_of_estimated_parameters,1);
else
prior_shape = evalin('base', 'bayestopt_.pshape');
end
beta_index = find(prior_shape==1);
if isempty(beta_index)
beta_draws = 0;
else
beta_draws = 1;
end
gamma_index = find(prior_shape==2);
if isempty(gamma_index)
gamma_draws = 0;
else
gamma_draws = 1;
end
gaussian_index = find(prior_shape==3);
if isempty(gaussian_index)
gaussian_draws = 0;
else
gaussian_draws = 1;
end
inverse_gamma_1_index = find(prior_shape==4);
if isempty(inverse_gamma_1_index)
inverse_gamma_1_draws = 0;
else
inverse_gamma_1_draws = 1;
end
uniform_index = find(prior_shape==5);
if isempty(uniform_index)
uniform_draws = 0;
else
uniform_draws = 1;
end
inverse_gamma_2_index = find(prior_shape==6);
if isempty(inverse_gamma_2_index)
inverse_gamma_2_draws = 0;
else
inverse_gamma_2_draws = 1;
end
pdraw = zeros(number_of_estimated_parameters,1);
return
end
if uniform_draws
pdraw(uniform_index) = rand(length(uniform_index),1).*(p4(uniform_index)-p3(uniform_index)) + p3(uniform_index);
end
if gaussian_draws
pdraw(gaussian_index) = randn(length(gaussian_index),1).*p7(gaussian_index) + p6(gaussian_index);
out_of_bound = find( (pdraw(gaussian_index)'>ub(gaussian_index)) | (pdraw(gaussian_index)'<lb(gaussian_index)));
while ~isempty(out_of_bound),
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));
out_of_bound = find( (pdraw(gaussian_index)'>ub(gaussian_index)) | (pdraw(gaussian_index)'<lb(gaussian_index)));
end
end
if gamma_draws
pdraw(gamma_index) = gamrnd(p6(gamma_index),p7(gamma_index))+p3(gamma_index);
out_of_bound = find( (pdraw(gamma_index)'>ub(gamma_index)) | (pdraw(gamma_index)'<lb(gamma_index)));
while ~isempty(out_of_bound),
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));
out_of_bound = find( (pdraw(gamma_index)'>ub(gamma_index)) | (pdraw(gamma_index)'<lb(gamma_index)));
end
end
if beta_draws
pdraw(beta_index) = (p4(beta_index)-p3(beta_index)).*betarnd(p6(beta_index),p7(beta_index))+p3(beta_index);
out_of_bound = find( (pdraw(beta_index)'>ub(beta_index)) | (pdraw(beta_index)'<lb(beta_index)));
while ~isempty(out_of_bound),
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));
out_of_bound = find( (pdraw(beta_index)'>ub(beta_index)) | (pdraw(beta_index)'<lb(beta_index)));
end
end
if inverse_gamma_1_draws
pdraw(inverse_gamma_1_index) = ...
sqrt(1./gamrnd(p7(inverse_gamma_1_index)/2,2./p6(inverse_gamma_1_index)))+p3(inverse_gamma_1_index);
out_of_bound = find( (pdraw(inverse_gamma_1_index)'>ub(inverse_gamma_1_index)) | (pdraw(inverse_gamma_1_index)'<lb(inverse_gamma_1_index)));
while ~isempty(out_of_bound),
pdraw(inverse_gamma_1_index(out_of_bound)) = ...
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));
out_of_bound = find( (pdraw(inverse_gamma_1_index)'>ub(inverse_gamma_1_index)) | (pdraw(inverse_gamma_1_index)'<lb(inverse_gamma_1_index)));
end
end
if inverse_gamma_2_draws
pdraw(inverse_gamma_2_index) = ...
1./gamrnd(p7(inverse_gamma_2_index)/2,2./p6(inverse_gamma_2_index))+p3(inverse_gamma_2_index);
out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'<lb(inverse_gamma_2_index)));
while ~isempty(out_of_bound),
pdraw(inverse_gamma_2_index(out_of_bound)) = ...
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));
out_of_bound = find( (pdraw(inverse_gamma_2_index)'>ub(inverse_gamma_2_index)) | (pdraw(inverse_gamma_2_index)'<lb(inverse_gamma_2_index)));
end
end
|
