diff --git a/license.txt b/license.txt index 0d947f8ae..034eb86ee 100644 --- a/license.txt +++ b/license.txt @@ -1,4 +1,4 @@ -Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/ +Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/ Upstream-Name: Dynare Upstream-Contact: Dynare Team, whose members in 2015 are: Stéphane Adjemian @@ -76,6 +76,11 @@ Copyright: 2008-2012 VZLU Prague, a.s. 2014 Dynare Team License: GPL-3+ +Files: matlab/one_sided_hp_filter.m +Copyright: 2010-2015 Alexander Meyer-Gohde + 2015 Dynare Team +License: GPL-3+ + Files: matlab/optimization/simpsa.m matlab/optimization/simpsaget.m matlab/optimization/simpsaset.m Copyright: 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m index 22fee8739..7a69585cb 100644 --- a/matlab/UnivariateSpectralDensity.m +++ b/matlab/UnivariateSpectralDensity.m @@ -109,7 +109,9 @@ ngrid = options_.hp_ngrid; %number of grid points freqs = (0 : pi/(ngrid-1):pi)'; % grid on which to compute tpos = exp( sqrt(-1)*freqs); %positive frequencies tneg = exp(-sqrt(-1)*freqs); %negative frequencies -if options_.hp_filter == 0 && ~options_.bandpass.indicator %do not filter +if options_.one_sided_hp_filter + error('UnivariateSpectralDensity:: spectral density estimate not available with one-sided HP filter') +elseif options_.hp_filter == 0 && ~options_.bandpass.indicator %do not filter filter_gain=ones(ngrid,1); elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && options_.bandpass.indicator %filter with bandpass filter_gain = zeros(1,ngrid); diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index 96ee22772..e6a9f3bcc 100644 --- a/matlab/disp_moments.m +++ b/matlab/disp_moments.m @@ -184,7 +184,7 @@ function y=get_filtered_time_series(y,m,options_) if options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator [hptrend,y] = sample_hp_filter(y,options_.hp_filter); elseif ~options_.hp_filter && options_.one_sided_hp_filter && ~options_.bandpass.indicator - error('disp_moments:: The one-sided HP filter is not yet available') + [hptrend,y] = one_sided_hp_filter(y,options_.one_sided_hp_filter); elseif ~options_.hp_filter && ~options_.one_sided_hp_filter && options_.bandpass.indicator data_temp=dseries(y,'0q1'); data_temp=baxter_king_filter(data_temp,options_.bandpass.passband(1),options_.bandpass.passband(2),200); diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index d93a86948..35406a01f 100644 --- a/matlab/disp_th_moments.m +++ b/matlab/disp_th_moments.m @@ -19,7 +19,9 @@ function oo_=disp_th_moments(dr,var_list,M_,options_,oo_) % along with Dynare. If not, see . nodecomposition = options_.nodecomposition; - +if options_.one_sided_hp_filter + error(['disp_th_moments:: theoretical moments incompatible with one-sided HP filter. Use simulated moments instead']) +end if size(var_list,1) == 0 var_list = M_.endo_names(1:M_.orig_endo_nbr, :); end diff --git a/matlab/one_sided_hp_filter.m b/matlab/one_sided_hp_filter.m new file mode 100644 index 000000000..38af32ea2 --- /dev/null +++ b/matlab/one_sided_hp_filter.m @@ -0,0 +1,113 @@ +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 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 . + +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 \ No newline at end of file diff --git a/tests/Makefile.am b/tests/Makefile.am index ee0d77182..18d54b0c1 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -11,6 +11,7 @@ MODFILES = \ moments/example1_hp_test.mod \ moments/example1_bp_test.mod \ moments/test_AR1_spectral_density.mod \ + moments/example1_one_sided_hp_test.mod \ gsa/ls2003.mod \ gsa/ls2003a.mod \ gsa/cod_ML_morris/cod_ML_morris.mod \ diff --git a/tests/moments/example1_one_sided_hp_test.mod b/tests/moments/example1_one_sided_hp_test.mod new file mode 100644 index 000000000..9da0798b5 --- /dev/null +++ b/tests/moments/example1_one_sided_hp_test.mod @@ -0,0 +1,71 @@ +/* + * Example 1 from F. Collard (2001): "Stochastic simulations with DYNARE: + * A practical guide" (see "guide.pdf" in the documentation directory). + */ + +/* + * Copyright (C) 2001-2010 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 . + */ + + +var y, c, k, a, h, b; +varexo e, u; + +parameters beta, rho, alpha, delta, theta, psi, tau; + +alpha = 0.36; +rho = 0.95; +tau = 0.025; +beta = 0.99; +delta = 0.025; +psi = 0; +theta = 2.95; + +phi = 0.1; + +model; +c*theta*h^(1+psi)=(1-alpha)*y; +k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1))) + *(exp(b(+1))*alpha*y(+1)+(1-delta)*k)); +y = exp(a)*(k(-1)^alpha)*(h^(1-alpha)); +k = exp(b)*(y-c)+(1-delta)*k(-1); +a = rho*a(-1)+tau*b(-1) + e; +b = tau*a(-1)+rho*b(-1) + u; +end; + +initval; +y = 1.08068253095672; +c = 0.80359242014163; +h = 0.29175631001732; +k = 11.08360443260358; +a = 0; +b = 0; +e = 0; +u = 0; +end; + +shocks; +var e; stderr 0.009; +var u; stderr 0.009; +var e, u = phi*0.009*0.009; +end; + +steady(solve_algo=4,maxit=1000); +options_.hp_ngrid=2048*4; + +stoch_simul(order=1,nofunctions,one_sided_hp_filter=1600,irf=0,periods=5000); \ No newline at end of file