Added one sided hp trend and cycle methods.

parent 2b6074a7
function o = onesidedhpcycle(o, lambda) % --*-- Unitary tests --*--
% Extracts the cycle component from a dseries object using the one sided HP filter.
%
% INPUTS
% - o [dseries] Original time series.
% - lambda [double] scalar, trend smoothness parameter.
%
% OUTPUTS
% - o [dseries] Cycle component of the original time series.
% Copyright (C) 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>1
if lambda<=0
error(['dseries::onesidedhptrend: Lambda must be a positive integer!'])
end
else
lambda = [];
end
o = copy(o);
o.onesidedhpcycle_(lambda);
return
%@test:1
d = randn(100,1);
o = dseries(d);
try
p = o.onesidedhpcycle();
t(1) = 1;
catch
t(1) = 0;
end
if t(1)
t(2) = dassert(o.data,d);
end
T = all(t);
%@eof:1
\ No newline at end of file
function o = onesidedhpcycle_(o, lambda) % --*-- Unitary tests --*--
% Extracts the cycle component from a dseries object using a one sided HP filter.
%
% INPUTS
% - o [dseries] Original time series.
% - lambda [double] scalar, trend smoothness parameter.
%
% OUTPUTS
% - o [dseries] Cycle component of the original time series.
% Copyright (C) 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>1
if lambda<=0
error(['dseries::onesidedhptrend: Lambda must be a positive integer!'])
end
else
lambda = [];
end
for i=1:vobs(o)
o.name(i) = {['onesidedhpcycle(' o.name{i} ')']};
o.tex(i) = {['\text{onesidedhpcycle}(' o.tex{i} ')']};
end
[~, o.data] = one_sided_hp_filter(o.data, lambda);
return
%@test:1
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;
try
ts0 = dseries(y,'1950Q1');
ts0.onesidedhpcycle_(1600);
t(1) = 1;
catch
t(1) = 0;
end
T = all(t);
%@eof:1
%@test:2
e = .2*randn(200,2);
u = randn(200,2);
stochastic_trends = cumsum(e);
deterministic_trends = transpose(1:200)*[.1, -.2];
x = zeros(200,2);
for i=2:200
x(i, 1) = .75*x(i-1, 1) + e(i, 1);
x(i, 2) = -.10*x(i-1, 2) + e(i, 2);
end
y = x + stochastic_trends + deterministic_trends;
try
ts0 = dseries(y,'1950Q1');
ts0.onesidedhpcycle_(1600);
t(1) = 1;
catch
t(1) = 0;
end
T = all(t);
%@eof:2
\ No newline at end of file
function o = onesidedhptrend(o, lambda) % --*-- Unitary tests --*--
% Extracts the trend component from a dseries object using the one sided HP filter.
%
% INPUTS
% - o [dseries] Original time series.
% - lambda [double] scalar, trend smoothness parameter.
%
% OUTPUTS
% - o [dseries] Trend component of the original time series.
% Copyright (C) 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>1
if lambda<=0
error(['dseries::onesidedhptrend: Lambda must be a positive integer!'])
end
else
lambda = [];
end
o = copy(o);
o.onesidedhptrend_(lambda);
return
%@test:1
d = randn(100,1);
o = dseries(d);
try
p = o.onesidedhptrend();
t(1) = 1;
catch
t(1) = 0;
end
if t(1)
t(2) = dassert(o.data,d);
end
T = all(t);
%@eof:1
\ No newline at end of file
function o = onesidedhptrend_(o, lambda) % --*-- Unitary tests --*--
% Extracts the trend component from a dseries object using a one sided HP filter.
%
% INPUTS
% - o [dseries] Original time series.
% - lambda [double] scalar, trend smoothness parameter.
%
% OUTPUTS
% - o [dseries] Trend component of the original time series.
% Copyright (C) 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>1
if lambda<=0
error(['dseries::onesidedhptrend: Lambda must be a positive integer!'])
end
else
lambda = [];
end
for i=1:vobs(o)
o.name(i) = {['onesidedhptrend(' o.name{i} ')']};
o.tex(i) = {['\text{onesidedhptrend}(' o.tex{i} ')']};
end
o.data = one_sided_hp_filter(o.data, lambda);
return
%@test:1
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;
try
ts0 = dseries(y,'1950Q1');
ts0.onesidedhptrend_(1600);
t(1) = 1;
catch
t(1) = 0;
end
T = all(t);
%@eof:1
%@test:2
e = .2*randn(200,2);
u = randn(200,2);
stochastic_trends = cumsum(e);
deterministic_trends = transpose(1:200)*[.1, -.2];
x = zeros(200,2);
for i=2:200
x(i, 1) = .75*x(i-1, 1) + e(i, 1);
x(i, 2) = -.10*x(i-1, 2) + e(i, 2);
end
y = x + stochastic_trends + deterministic_trends;
try
ts0 = dseries(y,'1950Q1');
ts0.onesidedhptrend_(1600);
t(1) = 1;
catch
t(1) = 0;
end
T = all(t);
%@eof:2
\ No newline at end of file
......@@ -90,7 +90,7 @@ switch S(1).type
case 'freq'
% Returns an integer characterizing the data frequency (1, 4, 12 or 52)
B = A.dates.freq;
case {'lag','lag_','lead','lead_','hptrend','hptrend_','hpcycle','hpcycle_','chain','chain_','detrend','detrend_','exist','mean','std','center','center_'} % Methods with less than two arguments.
case {'lag','lag_','lead','lead_','hptrend','hptrend_','hpcycle','hpcycle_','onesidedhptrend','onesidedhptrend_','onesidedhpcycle','onesidedhpcycle_','chain','chain_','detrend','detrend_','exist','mean','std','center','center_'} % Methods with less than two arguments.
if length(S)>1 && isequal(S(2).type,'()')
if isempty(S(2).subs)
B = feval(S(1).subs,A);
......
......@@ -71,6 +71,10 @@ if ~exist('randomstring','file')
p{end+1} = '/utilities/missing/randomstring';
end
if ~exist('one_sided_hp_filter','file')
p{end+1} = '/utilities/missing/one_sided_hp_filter';
end
% Install X13 binaries
opath = pwd();
cd([dseries_src_root '/../externals/x13'])
......
function [ytrend,ycycle]=one_sided_hp_filter(y,lambda,x_user,P_user,discard)
% function [ytrend,ycycle]=one_sided_hp_filter(y,lambda,x_user,P_user,discard)
% Conducts one-sided HP-filtering, derived using the Kalman filter
%
% Inputs:
% y [T*n] double data matrix in column format
% lambda [scalar] Smoothing parameter. Default value of 1600 will be used.
% x_user [2*n] double matrix with initial values of the state
% estimate for each variable in y. The underlying
% state vector is 2x1 for each variable in y.
% Default: use backwards extrapolations
% based on the first two observations
% P_user [n*1] struct structural array with n elements, each a two
% 2x2 matrix of intial MSE estimates for each
% variable in y.
% Default: matrix with large variances
% discard [scalar] number of initial periods to be discarded
% Default: 0
%
% Output:
% ytrend [(T-discard)*n] matrix of extracted trends
% ycycle [(T-discard)*n] matrix of extracted deviations from the extracted trends
%
% Algorithms:
%
% Implements the procedure described on p. 301 of Stock, J.H. and M.W. Watson (1999):
% "Forecasting inflation," Journal of Monetary Economics, vol. 44(2), pages 293-335, October.
% that states on page 301:
%
% "The one-sided HP trend estimate is constructed as the Kalman
% filter estimate of tau_t in the model:
%
% y_t=tau_t+epsilon_t
% (1-L)^2 tau_t=eta_t"
%
% The Kalman filter notation follows Chapter 13 of Hamilton, J.D. (1994).
% Time Series Analysis, with the exception of H, which is equivalent to his H'.
% Copyright (C) 200?-2015 Alexander Meyer-Gohde
% Copyright (C) 2015-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 < 2 || isempty(lambda)
lambda = 1600; %If the user didn't provide a value for lambda, set it to the default value 1600
end
[T,n] = size (y);% Calculate the number of periods and the number of variables in the series
%Set up state space
q=1/lambda; % the signal-to-noise ration: i.e. var eta_t / var epsilon_t
F=[2,-1;
1,0]; % state transition matrix
H=[1,0]; % observation matrix
Q=[q,0;
0,0]; % covariance matrix state equation errors
R=1; % variance observation equation error
for k=1:n %Run the Kalman filter for each variable
if nargin < 4 || isempty(x_user) %no intial value for state, extrapolate back two periods from the observations
x=[2*y(1,k)-y(2,k);
3*y(1,k)-2*y(2,k)];
else
x=x_user(:,k);
end
if nargin < 4 || isempty(P_user) %no initial value for the MSE, set a rather high one
P= [1e5 0;
0 1e5];
else
P=P_user{k};
end
for j=1:T %Get the estimates for each period
[x,P]=kalman_update(F,H,Q,R,y(j,k),x,P); %get new state estimate and update recursion
ytrend(j,k)=x(2);%second state is trend estimate
end
end
if nargout==2
ycycle=y-ytrend;
end
if nargin==5 %user provided a discard parameter
ytrend=ytrend(discard+1:end,:);%Remove the first "discard" periods from the trend series
if nargout==2
ycycle=ycycle(discard+1:end,:);
end
end
end
function [x,P]=kalman_update(F,H,Q,R,obs,x,P)
% Updates the Kalman filter estimation of the state and MSE
S=H*P*H'+R;
K=F*P*H';
K=K/S;
x=F*x+K*(obs -H*x); %State estimate
Temp=F-K*H;
P=Temp*P*Temp';
P=P+Q+K*R*K';%MSE estimate
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment