baxter_king_filter_.m 4.74 KB
Newer Older
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
function o = baxter_king_filter_(o, high_frequency, low_frequency, K) % --*-- Unitary tests --*--

% Implementation of Baxter and King (1999) band pass filter for dseries objects. The code is adapted from
% the one provided by Baxter and King. This filter isolates business cycle fluctuations with a period of length 
% ranging between high_frequency to low_frequency (quarters).
%
% INPUTS 
%  - o                  dseries object.
%  - high_frequency     positive scalar, period length (default value is 6).
%  - low_frequency      positive scalar, period length (default value is 32).
%  - K                  positive scalar integer, truncation parameter (default value is 12).
%
% OUTPUTS 
%  - o                  dseries object.
%
% REMARKS 
% This filter use a (symmetric) moving average smoother, so that K observations at the beginning and at the end of the 
% sample are lost in the computation of the filter.

% Copyright (C) 2013-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/>.

if nargin<4 || isempty(K)
    K = 12;
    if nargin<3 || isempty(low_frequency)
        % Set default number of periods corresponding to the lowest frequency.
        low_frequency = 32;
        if nargin<2 || isempty(high_frequency)
            % Set default number of periods corresponding to the highest frequency.
            high_frequency = 6;
            if nargin<1
                error('dseries::baxter_king_filter: I need at least one argument')
            end
        else
            if high_frequency<2
                error('dseries::baxter_king_filter: Second argument must be greater than 2!')
            end
            if high_frequency>low_frequency
                error('dseries::baxter_king_filter: Second argument must be smaller than the third argument!')
            end
        end
    end
end
       
% translate periods into frequencies.
hf=2.0*pi/high_frequency;
lf=2.0*pi/low_frequency;

% Set weights for the band-pass filter's lag polynomial.
weights = zeros(K+1,1); lpowers = transpose(1:K);
weights(2:K+1) = (sin(lpowers*hf)-sin(lpowers*lf))./(lpowers*pi);
weights(1) = (hf-lf)/pi;

% Set the constraint on the sum of weights.
if low_frequency>1000
    % => low pass filter.
    sum_of_weights_constraint = 1.0;
else
    sum_of_weights_constraint = 0.0;
end

% Compute the sum of weights.
sum_of_weights = weights(1) + 2*sum(weights(2:K+1));

% Correct the weights.
weights = weights + (sum_of_weights_constraint - sum_of_weights)/(2*K+1);

% Weights are symmetric!
weights = [flipud(weights(2:K+1)); weights];

tmp = zeros(size(o.data));

% Filtering step.
for t = K+1:nobs(o)-K
    tmp(t,:)  = weights'*o.data(t-K:t+K,:);    
end

% Update dseries object.
o.data = tmp(K+1:end-K,:);
init = firstdate(o)+K;
o.dates = init:init+(nobs(o)-1);

%@test:1
%$ plot_flag = false;
%$
%$ % Create a dataset.
%$ e = .2*randn(200,1);
%$ u = randn(200,1);
%$ stochastic_trend = cumsum(e); 
%$ deterministic_trend = .1*transpose(1:200);
%$ x = zeros(200,1);
%$ for i=2:200
%$    x(i) = .75*x(i-1) + e(i);
%$ end
%$ y = x + stochastic_trend + deterministic_trend;
%$
%$ % Test the routine.
%$ try
%$     ts = dseries(y,'1950Q1');
%$     ts = ts.baxter_king_filter_();
%$     xx = dseries(x,'1950Q1');
%$     t(1) = 1;
%$ catch
%$     t(1) = 0;
%$ end
%$
%$ if t(1)
%$     t(2) = dassert(ts.freq,4);
%$     t(3) = dassert(ts.init.freq,4);
%$     t(4) = dassert(ts.init.time,[1953, 1]);
%$     t(5) = dassert(ts.vobs,1);
%$     t(6) = dassert(ts.nobs,176);
%$ end
%$
%$ % Show results
%$ if plot_flag
%$     plot(xx(ts.dates).data,'-k');
%$     hold on
%$     plot(ts.data,'--r');
%$     hold off
%$     axis tight
%$     id = get(gca,'XTick');
%$     set(gca,'XTickLabel',strings(ts.dates(id)));
%$     legend({'Stationary component of y', 'Filtered y'})
%$     print('-depsc2','../doc/dynare.plots/BaxterKingFilter.eps')
%$     system('convert -density 300 ../doc/dynare.plots/BaxterKingFilter.eps ../doc/dynare.plots/BaxterKingFilter.png');
%$     system('convert -density 300 ../doc/dynare.plots/BaxterKingFilter.eps ../doc/dynare.plots/BaxterKingFilter.pdf');
%$     system('convert -density 300 ../doc/dynare.plots/BaxterKingFilter.eps ../doc/dynare.plots/BaxterKingFilter.jpg');
%$ end
%$
%$ T = all(t);
%@eof:1