1function [ytrend,ycycle]=one_sided_hp_filter(y,lambda,x_user,P_user,discard)
2% function [ytrend,ycycle]=one_sided_hp_filter(y,lambda,x_user,P_user,discard)
3% Conducts one-sided HP-filtering, derived using the Kalman filter
4%
5% Inputs:
6%   y           [T*n] double    data matrix in column format
7%   lambda      [scalar]        Smoothing parameter. Default value of 1600 will be used.
8%   x_user      [2*n] double    matrix with initial values of the state
9%                               estimate for each variable in y. The underlying
10%                               state vector is 2x1 for each variable in y.
11%                               Default: use backwards extrapolations
12%                               based on the first two observations
13%   P_user      [n*1] struct    structural array with n elements, each a two
14%                               2x2 matrix of intial MSE estimates for each
15%                               variable in y.
16%                               Default: matrix with large variances
17%   discard     [scalar]        number of initial periods to be discarded
18%                               Default: 0
19%
20% Output:
21%   ytrend      [(T-discard)*n] matrix of extracted trends
22%   ycycle      [(T-discard)*n] matrix of extracted deviations from the extracted trends
23%
24% Algorithms:
25%
26%   Implements the procedure described on p. 301 of Stock, J.H. and M.W. Watson (1999):
27%   "Forecasting inflation," Journal of Monetary Economics,  vol. 44(2), pages 293-335, October.
28%   that states on page 301:
29%
30%       "The one-sided HP trend estimate is constructed as the Kalman
31%       filter estimate of tau_t in the model:
32%
33%       y_t=tau_t+epsilon_t
34%       (1-L)^2 tau_t=eta_t"
35%
36%  The Kalman filter notation follows Chapter 13 of Hamilton, J.D. (1994).
37%   Time Series Analysis, with the exception of H, which is equivalent to his H'.
38
39
40% Copyright (C) 2010-2015 Alexander Meyer-Gohde
41% Copyright (C) 2015-2017 Dynare Team
42%
43% This file is part of Dynare.
44%
45% Dynare is free software: you can redistribute it and/or modify
46% it under the terms of the GNU General Public License as published by
47% the Free Software Foundation, either version 3 of the License, or
48% (at your option) any later version.
49%
50% Dynare is distributed in the hope that it will be useful,
51% but WITHOUT ANY WARRANTY; without even the implied warranty of
52% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
53% GNU General Public License for more details.
54%
55% You should have received a copy of the GNU General Public License
56% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
57
58if nargin < 2 || isempty(lambda)
59    lambda = 1600; %If the user didn't provide a value for lambda, set it to the default value 1600
60end
61[T,n] = size (y);% Calculate the number of periods and the number of variables in the series
62
63%Set up state space
64q=1/lambda;     % the signal-to-noise ration: i.e. var eta_t / var epsilon_t
65F=[2,-1;
66   1,0];       % state transition matrix
67H=[1,0];        % observation matrix
68Q=[q,0;
69   0,0];        % covariance matrix state equation errors
70R=1;            % variance observation equation error
71
72for k=1:n %Run the Kalman filter for each variable
73    if nargin < 3 || isempty(x_user) %no intial value for state, extrapolate back two periods from the observations
74        x=[2*y(1,k)-y(2,k);
75           3*y(1,k)-2*y(2,k)];
76    else
77        x=x_user(:,k);
78    end
79    if nargin < 4 || isempty(P_user) %no initial value for the MSE, set a rather high one
80        P= [1e5 0;
81            0 1e5];
82    else
83        P=P_user{k};
84    end
85
86    for j=1:T %Get the estimates for each period
87        [x,P]=kalman_update(F,H,Q,R,y(j,k),x,P); %get new state estimate and update recursion
88        ytrend(j,k)=x(2);%second state is trend estimate
89    end
90end
91
92if nargout==2
93    ycycle=y-ytrend;
94end
95
96if nargin==5 %user provided a discard parameter
97    ytrend=ytrend(discard+1:end,:);%Remove the first "discard" periods from the trend series
98    if nargout==2
99        ycycle=ycycle(discard+1:end,:);
100    end
101end
102end
103
104function   [x,P]=kalman_update(F,H,Q,R,obs,x,P)
105% Updates the Kalman filter estimation of the state and MSE
106S=H*P*H'+R;
107K=F*P*H';
108K=K/S;
109x=F*x+K*(obs -H*x); %State estimate
110Temp=F-K*H;
111P=Temp*P*Temp';
112P=P+Q+K*R*K';%MSE estimate
113end
114