kernel_density_estimate

 Estimates a continuous density. 
 
 INPUTS
   data                  [double]  Vector (number_of_draws*1) of draws.
   number_of_grid_points [integer] Scalar, number of points where the density is estimated.
                                   This (positive) integer must be a power of two.  
   number_of_draws       [integer] Scalar, number of draws.
   bandwidth             [double]  Real positive scalar.                                  
   kernel_function       [string]  Name of the kernel function: 'gaussian', 'triweight',
                                   'uniform', 'triangle', 'epanechnikov', 'quartic', 
                                   'triweight' and 'cosinus'

 OUTPUTS
    abscissa             [double] Vector (number_of_grid_points*1) of values on the abscissa axis.
    f:                   [double] Vector (number_of_grid_points*1) of values on the ordinate axis, 
                                  (density estimates).
        
 SPECIAL REQUIREMENTS
    none.

 REFERENCES
    A kernel density estimator is used (see Silverman [1986], "Density estimation for statistics and data analysis")
    The code is adapted from Anders Holtsberg's matlab toolbox (stixbox).

0001 function [abscissa,f] = kernel_density_estimate(data,number_of_grid_points,number_of_draws,bandwidth,kernel_function) 
0002 % Estimates a continuous density.
0003 %
0004 % INPUTS
0005 %   data                  [double]  Vector (number_of_draws*1) of draws.
0006 %   number_of_grid_points [integer] Scalar, number of points where the density is estimated.
0007 %                                   This (positive) integer must be a power of two.
0008 %   number_of_draws       [integer] Scalar, number of draws.
0009 %   bandwidth             [double]  Real positive scalar.
0010 %   kernel_function       [string]  Name of the kernel function: 'gaussian', 'triweight',
0011 %                                   'uniform', 'triangle', 'epanechnikov', 'quartic',
0012 %                                   'triweight' and 'cosinus'
0013 %
0014 % OUTPUTS
0015 %    abscissa             [double] Vector (number_of_grid_points*1) of values on the abscissa axis.
0016 %    f:                   [double] Vector (number_of_grid_points*1) of values on the ordinate axis,
0017 %                                  (density estimates).
0018 %
0019 % SPECIAL REQUIREMENTS
0020 %    none.
0021 %
0022 % REFERENCES
0023 %    A kernel density estimator is used (see Silverman [1986], "Density estimation for statistics and data analysis")
0024 %    The code is adapted from Anders Holtsberg's matlab toolbox (stixbox).
0025 %
0026 
0027 % Copyright (C) 2004-2008 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 if min(size(data))>1
0045     error('kernel_density_estimate:: data must be a one dimensional array.');
0046 else
0047     data = data(:);
0048 end
0049 
0050 test = log(number_of_grid_points)/log(2);
0051 if (abs(test-round(test)) > 1e-12)
0052     error('kernel_density_estimate:: The number of grid points must be a power of 2.');
0053 end
0054 
0055 %% Kernel specification.
0056 if strcmpi(kernel_function,'gaussian') 
0057     kernel = @(x) inv(sqrt(2*pi))*exp(-0.5*x.^2);
0058 elseif strcmpi(kernel_function,'uniform') 
0059     kernel = @(x) 0.5*(abs(x) <= 1); 
0060 elseif strcmpi(kernel_function,'triangle') 
0061     kernel = @(x) (1-abs(x)).*(abs(x) <= 1);
0062 elseif strcmpi(kernel_function,'epanechnikov') 
0063     kernel = @(x) 0.75*(1-x.^2).*(abs(x) <= 1);
0064 elseif strcmpi(kernel_function,'quartic') 
0065     kernel = @(x) 0.9375*((1-x.^2).^2).*(abs(x) <= 1);
0066 elseif strcmpi(kernel_function,'triweight') 
0067     kernel = @(x) 1.09375*((1-x.^2).^3).*(abs(x) <= 1);
0068 elseif strcmpi(kernel_function,'cosinus') 
0069     kernel = @(x) (pi/4)*cos((pi/2)*x).*(abs(x) <= 1);
0070 end
0071 
0072 %% Non parametric estimation (Gaussian kernel should be used (FFT)).
0073 lower_bound  = min(data) - (max(data)-min(data))/3;
0074 upper_bound  = max(data) + (max(data)-min(data))/3;
0075 abscissa = linspace(lower_bound,upper_bound,number_of_grid_points)';
0076 inc = abscissa(2)-abscissa(1); 
0077 xi  = zeros(number_of_grid_points,1);
0078 xa  = (data-lower_bound)/(upper_bound-lower_bound)*number_of_grid_points; 
0079 for i = 1:number_of_draws
0080     indx = floor(xa(i));
0081     temp = xa(i)-indx;
0082     xi(indx+[1 2]) = xi(indx+[1 2]) + [1-temp,temp]';
0083 end
0084 xk = [-number_of_grid_points:number_of_grid_points-1]'*inc;
0085 kk = kernel(xk/bandwidth);
0086 kk = kk / (sum(kk)*inc*number_of_draws);
0087 f = ifft(fft(fftshift(kk)).*fft([ xi ; zeros(number_of_grid_points,1) ]));
0088 f = real(f(1:number_of_grid_points));