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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
function optimal_bandwidth = mh_optimal_bandwidth(data,number_of_draws,bandwidth,kernel_function)
% This function gives the optimal bandwidth parameter of a kernel estimator
% used to estimate a posterior univariate density from realisations of a
% Metropolis-Hastings algorithm.
%
% INPUTS:
% data [double] Vector (number_of_draws*1) of draws.
% number_of_draws [integer] Scalar, number of draws.
% bandwidth [integer] Scalar equal to 0,-1 or -2.
% bandwidth = 0 => Silverman [1986] rule of thumb.
% bandwidth = -1 => Sheather and Jones [1991].
% bandwidth = -2 => Local bandwith parameters.
% kernel_function [string] Name of the kernel function: 'gaussian', 'triweight',
% 'uniform', 'triangle', 'epanechnikov', 'quartic',
% 'triweight' and 'cosinus'
%
% OUTPUTS:
% optimal_bandwidth: [double] Scalar or vector, optimal window width.
%
% SPECIAL REQUIREMENTS:
% none.
%
% REFERENCES:
% [1] M. Skoeld and G.O. Roberts [2003], "Density estimation for the Metropolis-Hastings algorithm".
% [2] Silverman [1986], "Density estimation for statistics and data analysis".
% Copyright (C) 2004-2017 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/>.
% Kernel specifications.
if strcmpi(kernel_function,'gaussian')
% Kernel definition
k = @(x)inv(sqrt(2*pi))*exp(-0.5*x.^2);
% Second derivate of the kernel function
k2 = @(x)inv(sqrt(2*pi))*(-exp(-0.5*x.^2)+(x.^2).*exp(-0.5*x.^2));
% Fourth derivate of the kernel function
k4 = @(x)inv(sqrt(2*pi))*(3*exp(-0.5*x.^2)-6*(x.^2).*exp(-0.5*x.^2)+(x.^4).*exp(-0.5*x.^2));
% Sixth derivate of the kernel function
k6 = @(x)inv(sqrt(2*pi))*(-15*exp(-0.5*x.^2)+45*(x.^2).*exp(-0.5*x.^2)-15*(x.^4).*exp(-0.5*x.^2)+(x.^6).*exp(-0.5*x.^2));
mu02 = inv(2*sqrt(pi));
mu21 = 1;
elseif strcmpi(kernel_function,'uniform')
k = @(x)0.5*(abs(x) <= 1);
mu02 = 0.5;
mu21 = 1/3;
elseif strcmpi(kernel_function,'triangle')
k = @(x)(1-abs(x)).*(abs(x) <= 1);
mu02 = 2/3;
mu21 = 1/6;
elseif strcmpi(kernel_function,'epanechnikov')
k = @(x)0.75*(1-x.^2).*(abs(x) <= 1);
mu02 = 3/5;
mu21 = 1/5;
elseif strcmpi(kernel_function,'quartic')
k = @(x)0.9375*((1-x.^2).^2).*(abs(x) <= 1);
mu02 = 15/21;
mu21 = 1/7;
elseif strcmpi(kernel_function,'triweight')
k = @(x)1.09375*((1-x.^2).^3).*(abs(x) <= 1);
k2 = @(x)(105/4*(1-x.^2).*x.^2-105/16*(1-x.^2).^2).*(abs(x) <= 1);
k4 = @(x)(-1575/4*x.^2+315/4).*(abs(x) <= 1);
k6 = @(x)(-1575/2).*(abs(x) <= 1);
mu02 = 350/429;
mu21 = 1/9;
elseif strcmpi(kernel_function,'cosinus')
k = @(x)(pi/4)*cos((pi/2)*x).*(abs(x) <= 1);
k2 = @(x)(-1/16*cos(pi*x/2)*pi^3).*(abs(x) <= 1);
k4 = @(x)(1/64*cos(pi*x/2)*pi^5).*(abs(x) <= 1);
k6 = @(x)(-1/256*cos(pi*x/2)*pi^7).*(abs(x) <= 1);
mu02 = (pi^2)/16;
mu21 = (pi^2-8)/pi^2;
else
disp('mh_optimal_bandwidth:: ');
error('This kernel function is not yet implemented in Dynare!');
end
% Get the Skold and Roberts' correction.
if bandwidth==0 || bandwidth==-1
correction = correction_for_repeated_draws(data,number_of_draws);
else
correction = 0;
end
% Compute the standard deviation of the draws.
sigma = std(data);
% Optimal bandwidth parameter.
if bandwidth == 0 % Rule of thumb bandwidth parameter (Silverman [1986].
h = 2*sigma*(sqrt(pi)*mu02/(12*(mu21^2)*number_of_draws))^(1/5);
h = h*correction^(1/5);
elseif bandwidth == -1 % Sheather and Jones [1991] plug-in estimation of the optimal bandwidth parameter.
if strcmp(kernel_function,'uniform') || ...
strcmp(kernel_function,'triangle') || ...
strcmp(kernel_function,'epanechnikov') || ...
strcmp(kernel_function,'quartic')
error(['I can''t compute the optimal bandwidth with this kernel...' ...
'Try the gaussian, triweight or cosinus kernels.']);
end
Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi));
g3 = abs(2*correction*k6(0)/(mu21*Itilda4*number_of_draws))^(1/9);
Ihat3 = 0;
for i=1:number_of_draws
Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3));
end
Ihat3 = - Ihat3/((number_of_draws^2)*g3^7);
g2 = abs(2*correction*k4(0)/(mu21*Ihat3*number_of_draws))^(1/7);
Ihat2 = 0;
for i=1:number_of_draws
Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2));
end
Ihat2 = Ihat2/((number_of_draws^2)*g2^5);
h = (correction*mu02/(number_of_draws*Ihat2*mu21^2))^(1/5); % equation (22) in Skold and Roberts [2003].
elseif bandwidth == -2 % Bump killing... I compute local bandwith parameters in order to remove
% spurious bumps introduced by long rejecting periods.
if strcmp(kernel_function,'uniform') || ...
strcmp(kernel_function,'triangle') || ...
strcmp(kernel_function,'epanechnikov') || ...
strcmp(kernel_function,'quartic')
error(['I can''t compute the optimal bandwidth with this kernel...' ...
'Try the gaussian, triweight or cosinus kernels.']);
end
T = zeros(number_of_draws, 1);
for i=1:number_of_draws
j = i;
while j<=number_of_draws && (data(j,1)-data(i,1))<2*eps
j = j+1;
end
T(i) = (j-i);
correction = correction + 2*T(i) - 1;
end
correction = correction/number_of_draws;
Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi));
g3 = abs(2*correction*k6(0)/(mu21*Itilda4*number_of_draws))^(1/9);
Ihat3 = 0;
for i=1:number_of_draws
Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3));
end
Ihat3 = -Ihat3/((n^2)*g3^7);
g2 = abs(2*correction*k4(0)/(mu21*Ihat3*number_of_draws))^(1/7);
Ihat2 = 0;
for i=1:number_of_draws
Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2));
end
Ihat2 = Ihat2/((number_of_draws^2)*g2^5);
h = ((2*T-1)*mu02/(number_of_draws*Ihat2*mu21^2)).^(1/5); % h is a column vector (local banwidth parameters).
else
disp('mh_optimal_bandwidth:: ');
error('Parameter bandwidth must be equal to 0, -1 or -2.');
end
optimal_bandwidth = h;
return
function correction = correction_for_repeated_draws(draws,n)
correction = 0;
for i=1:n
j = i;
while j<=n && ( draws(j,1) - draws(i,1) )<2*eps
j = j+1;
end
correction = correction + 2*(j-i) - 1;
end
correction = correction/n;
|
