mh_optimal_bandwidth

 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. Skold and G.O. Roberts [2003], "Density estimation for the Metropolis-Hastings algorithm". 
   [2] Silverman [1986], "Density estimation for statistics and data analysis".

0001 function optimal_bandwidth = mh_optimal_bandwidth(data,number_of_draws,bandwidth,kernel_function) 
0002 % This function gives the optimal bandwidth parameter of a kernel estimator
0003 % used to estimate a posterior univariate density from realisations of a
0004 % Metropolis-Hastings algorithm.
0005 %
0006 % INPUTS:
0007 %   data               [double]  Vector (number_of_draws*1) of draws.
0008 %   number_of_draws    [integer] Scalar, number of draws.
0009 %   bandwidth          [integer] Scalar equal to 0,-1 or -2.
0010 %                                bandwidth =  0 => Silverman [1986] rule of thumb.
0011 %                                bandwidth = -1 => Sheather and Jones [1991].
0012 %                                bandwidth = -2 => Local bandwith parameters.
0013 %   kernel_function    [string]  Name of the kernel function: 'gaussian', 'triweight',
0014 %                                'uniform', 'triangle', 'epanechnikov', 'quartic',
0015 %                                'triweight' and 'cosinus'
0016 %
0017 % OUTPUTS:
0018 %   optimal_bandwidth: [double]  Scalar or vector, optimal window width.
0019 %
0020 % SPECIAL REQUIREMENTS:
0021 %   none.
0022 %
0023 % REFERENCES:
0024 %   [1] M. Skold and G.O. Roberts [2003], "Density estimation for the Metropolis-Hastings algorithm".
0025 %   [2] Silverman [1986], "Density estimation for statistics and data analysis".
0026 
0027 % Copyright (C) 2004-2011 Dynare Team
0028 %
0029 % This file is part of Dynare.
0030 %
0031 % Dynare is free software: you can redistribute it and/or modify
0032 % it under the terms of the GNU General Public License as published by
0033 % the Free Software Foundation, either version 3 of the License, or
0034 % (at your option) any later version.
0035 %
0036 % Dynare is distributed in the hope that it will be useful,
0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 % GNU General Public License for more details.
0040 %
0041 % You should have received a copy of the GNU General Public License
0042 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0043 
0044 %% Kernel specifications.
0045 if strcmpi(kernel_function,'gaussian')
0046     % Kernel definition
0047     k    = inline('inv(sqrt(2*pi))*exp(-0.5*x.^2)');
0048     % Second derivate of the kernel function
0049     k2   = inline('inv(sqrt(2*pi))*(-exp(-0.5*x.^2)+(x.^2).*exp(-0.5*x.^2))');
0050     % Fourth derivate of the kernel function
0051     k4   = inline('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))');
0052     % Sixth derivate of the kernel function
0053     k6   = inline(['inv(sqrt(2*pi))*(-15*exp(-0.5*x.^2)+45*(x.^2).*exp(-' ...
0054                    '0.5*x.^2)-15*(x.^4).*exp(-0.5*x.^2)+(x.^6).*exp(-0.5*x.^2))']);
0055     mu02 = inv(2*sqrt(pi));
0056     mu21 = 1;
0057 elseif strcmpi(kernel_function,'uniform') 
0058     k    = inline('0.5*(abs(x) <= 1)');
0059     mu02 = 0.5;
0060     mu21 = 1/3;
0061 elseif strcmpi(kernel_function,'triangle') 
0062     k    = inline('(1-abs(x)).*(abs(x) <= 1)');
0063     mu02 = 2/3;
0064     mu21 = 1/6;
0065 elseif strcmpi(kernel_function,'epanechnikov') 
0066     k    = inline('0.75*(1-x.^2).*(abs(x) <= 1)');
0067     mu02 = 3/5;
0068     mu21 = 1/5;
0069 elseif strcmpi(kernel_function,'quartic') 
0070     k    = inline('0.9375*((1-x.^2).^2).*(abs(x) <= 1)');
0071     mu02 = 15/21;
0072     mu21 = 1/7;
0073 elseif strcmpi(kernel_function,'triweight') 
0074     k    = inline('1.09375*((1-x.^2).^3).*(abs(x) <= 1)');
0075     k2   = inline('(105/4*(1-x.^2).*x.^2-105/16*(1-x.^2).^2).*(abs(x) <= 1)');
0076     k4   = inline('(-1575/4*x.^2+315/4).*(abs(x) <= 1)');
0077     k6   = inline('(-1575/2).*(abs(x) <= 1)');
0078     mu02 = 350/429;
0079     mu21 = 1/9;
0080 elseif strcmpi(kernel_function,'cosinus') 
0081     k    = inline('(pi/4)*cos((pi/2)*x).*(abs(x) <= 1)');
0082     k2   = inline('(-1/16*cos(pi*x/2)*pi^3).*(abs(x) <= 1)');
0083     k4   = inline('(1/64*cos(pi*x/2)*pi^5).*(abs(x) <= 1)');
0084     k6   = inline('(-1/256*cos(pi*x/2)*pi^7).*(abs(x) <= 1)');
0085     mu02 = (pi^2)/16;
0086     mu21 = (pi^2-8)/pi^2;
0087 else
0088     disp('mh_optimal_bandwidth:: ');
0089     error('This kernel function is not yet implemented in Dynare!');
0090 end
0091 
0092 
0093 %% Get the Skold and Roberts' correction.
0094 if bandwidth==0 || bandwidth==-1
0095     correction = correction_for_repeated_draws(data,number_of_draws);
0096 else
0097     correction = 0;
0098 end
0099 %% Compute the standard deviation of the draws.
0100 sigma = std(data);
0101 %% Optimal bandwidth parameter.
0102 if bandwidth == 0;  % Rule of thumb bandwidth parameter (Silverman [1986].
0103     h = 2*sigma*(sqrt(pi)*mu02/(12*(mu21^2)*number_of_draws))^(1/5);
0104     h = h*correction^(1/5);
0105 elseif bandwidth == -1; % Sheather and Jones [1991] plug-in estimation of the optimal bandwidth parameter.
0106     if strcmp(kernel_function,'uniform')      || ... 
0107             strcmp(kernel_function,'triangle')     || ... 
0108             strcmp(kernel_function,'epanechnikov') || ... 
0109             strcmp(kernel_function,'quartic')
0110         error(['I can''t compute the optimal bandwidth with this kernel...' ...
0111                'Try the gaussian, triweight or cosinus kernels.']);
0112     end 
0113     Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi));
0114     g3      = abs(2*correction*k6(0)/(mu21*Itilda4*number_of_draws))^(1/9);
0115     Ihat3   = 0;
0116     for i=1:number_of_draws
0117         Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3));
0118     end     
0119     Ihat3 = - Ihat3/((number_of_draws^2)*g3^7);
0120     g2    = abs(2*correction*k4(0)/(mu21*Ihat3*number_of_draws))^(1/7);
0121     Ihat2 = 0;
0122     for i=1:number_of_draws
0123         Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2));
0124     end
0125     Ihat2 = Ihat2/((number_of_draws^2)*g2^5);
0126     h     = (correction*mu02/(number_of_draws*Ihat2*mu21^2))^(1/5); % equation (22) in Skold and Roberts [2003].
0127 elseif bandwidth == -2;     % Bump killing... I compute local bandwith parameters in order to remove
0128                             % spurious bumps introduced by long rejecting periods.
0129     if strcmp(kernel_function,'uniform')      || ... 
0130             strcmp(kernel_function,'triangle')     || ... 
0131             strcmp(kernel_function,'epanechnikov') || ... 
0132             strcmp(kernel_function,'quartic')
0133         error(['I can''t compute the optimal bandwidth with this kernel...' ...
0134                'Try the gaussian, triweight or cosinus kernels.']);
0135     end
0136     T = zeros(n,1);
0137     for i=1:n
0138         j = i;
0139         while j<= n && (data(j,1)-data(i,1))<2*eps;
0140             j = j+1;
0141         end     
0142         T(i) = (j-i);
0143         correction = correction + 2*T(i) - 1;
0144     end
0145     correction = correction/number_of_draws;
0146     Itilda4 = 8*7*6*5/(((2*sigma)^9)*sqrt(pi));
0147     g3      = abs(2*correction*k6(0)/(mu21*Itilda4*correction))^(1/9);
0148     Ihat3   = 0;
0149     for i=1:number_of_draws
0150         Ihat3 = Ihat3 + sum(k6((data(i,1)-data)/g3));
0151     end
0152     Ihat3 = -Ihat3/((n^2)*g3^7);
0153     g2    = abs(2*correction*k4(0)/(mu21*Ihat3*n))^(1/7);
0154     Ihat2 = 0;
0155     for i=1:number_of_draws;
0156         Ihat2 = Ihat2 + sum(k4((data(i)-data)/g2));
0157     end     
0158     Ihat2 = Ihat2/((number_of_draws^2)*g2^5);
0159     h = ((2*T-1)*mu02/(number_of_draws*Ihat2*mu21^2)).^(1/5); % h is a column vector (local banwidth parameters).
0160 else
0161     disp('mh_optimal_bandwidth:: ');
0162     error('Parameter bandwidth must be equal to 0, -1 or -2.');
0163 end
0164 
0165 optimal_bandwidth = h;
0166 return,
0167 
0168 
0169 function correction = correction_for_repeated_draws(draws,n)
0170 correction = 0;
0171 for i=1:n
0172     j = i;
0173     while j<=n && ( draws(j,1) - draws(i,1) )<2*eps; 
0174         j = j+1;
0175     end
0176     correction = correction + 2*(j-i) - 1;
0177 end
0178 correction = correction/n;