From e9f2be0c1a2d101bac0666e6b5f59f5e83b05ace Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sat, 1 Aug 2015 18:42:45 +0200 Subject: [PATCH 01/14] Cosmetic changes to UnivariateSpectralDensity.m --- matlab/UnivariateSpectralDensity.m | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m index 9aae99c4c..95e8dab69 100644 --- a/matlab/UnivariateSpectralDensity.m +++ b/matlab/UnivariateSpectralDensity.m @@ -104,19 +104,16 @@ bb = ghu(iky,:); if options_.hp_filter == 0 tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb'; - k = find(abs(tmp) < 1e-12); - tmp(k) = 0; + tmp(abs(tmp) < 1e-12) = 0; Gamma(:,1) = diag(tmp); vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb'); tmp = aa*vxy; - k = find(abs(tmp) < 1e-12); - tmp(k) = 0; + tmp(abs(tmp) < 1e-12) = 0; Gamma(:,2) = diag(tmp); for i=2:cutoff vxy = A*vxy; tmp = aa*vxy; - k = find(abs(tmp) < 1e-12); - tmp(k) = 0; + tmp(abs(tmp) < 1e-12) = 0; Gamma(:,i+1) = diag(tmp); end else @@ -129,28 +126,26 @@ else tpos = exp( sqrt(-1)*freqs); tneg = exp(-sqrt(-1)*freqs); hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2); - mathp_col = []; + mathp_col = NaN(ngrid,length(ivar)^2); IA = eye(size(A,1)); IE = eye(M_.exo_nbr); for ig = 1:ngrid - f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]... - *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) IE]); % state variables + f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\ghu1;IE]... + *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % state variables g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series - mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row + mathp_col(ig,:) = (f_hp(:))'; % store as matrix row % for ifft end; imathp_col = real(ifft(mathp_col))*(2*pi); tmp = reshape(imathp_col(1,:),nvar,nvar); - k = find(abs(tmp)<1e-12); - tmp(k) = 0; + tmp(abs(tmp)<1e-12) = 0; Gamma(:,1) = diag(tmp); sy = sqrt(Gamma(:,1)); sy = sy *sy'; for i=1:cutoff-1 tmp = reshape(imathp_col(i+1,:),nvar,nvar)./sy; - k = find(abs(tmp) < 1e-12); - tmp(k) = 0; + tmp(abs(tmp) < 1e-12) = 0; Gamma(:,i+1) = diag(tmp); end end From 25df325899a225351c182a6bac5be6ebe2fbece5 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sat, 1 Aug 2015 18:43:12 +0200 Subject: [PATCH 02/14] Add bandpass filter to th_autocovariances.m Also adds documentation --- matlab/th_autocovariances.m | 83 +++++++++++++++++++------------------ 1 file changed, 43 insertions(+), 40 deletions(-) diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m index 6852f9654..31ee5025c 100644 --- a/matlab/th_autocovariances.m +++ b/matlab/th_autocovariances.m @@ -1,8 +1,8 @@ function [Gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition) % Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process -% with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_ -% for a subset of variables ivar (indices in lgy_) -% Theoretical HP-filtering is available as an option +% with coefficients dr.ghx and dr.ghu and shock variances Sigma_e +% for a subset of variables ivar. +% Theoretical HP-filtering and band-pass filtering is available as an option % % INPUTS % dr: [structure] Reduced form solution of the DSGE model (decisions rules) @@ -156,7 +156,7 @@ if options_.order == 2 || options_.hp_filter == 0 end end end -if options_.hp_filter == 0 +if options_.hp_filter == 0 && ~options_.bandpass.indicator v = NaN*ones(nvar,nvar); v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb'; k = find(abs(v) < 1e-12); @@ -207,37 +207,44 @@ if options_.hp_filter == 0 end end end -else% ==> Theoretical HP filter. +else% ==> Theoretical filters. % By construction, all variables are stationary when HP filtered iky = inv_order_var(ivar); stationary_vars = (1:length(ivar))'; - aa = ghx(iky,:); - bb = ghu(iky,:); + aa = ghx(iky,:); %R in Uhlig (2001) + bb = ghu(iky,:); %S in Uhlig (2001) lambda = options_.hp_filter; ngrid = options_.hp_ngrid; - freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); - tpos = exp( sqrt(-1)*freqs); - tneg = exp(-sqrt(-1)*freqs); - hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2); - mathp_col = []; + freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); %[0,2*pi) + tpos = exp( sqrt(-1)*freqs); %positive frequencies + tneg = exp(-sqrt(-1)*freqs); %negative frequencies + if options_.bandpass.indicator + filter_gain = zeros(1,ngrid); + lowest_periodicity=options_.bandpass.passband(2); + highest_periodicity=options_.bandpass.passband(1); + highest_periodicity=max(2,highest_periodicity); % restrict to upper bound of pi + filter_gain(freqs>=2*pi/lowest_periodicity & freqs<=2*pi/highest_periodicity)=1; + filter_gain(freqs<=-2*pi/lowest_periodicity+2*pi & freqs>=-2*pi/highest_periodicity+2*pi)=1; + else + filter_gain = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2); %HP transfer function + end + mathp_col = NaN(ngrid,length(ivar)^2); IA = eye(size(A,1)); IE = eye(M_.exo_nbr); for ig = 1:ngrid - if hp1(ig)==0, + if filter_gain(ig)==0, f_hp = zeros(length(ivar),length(ivar)); else - f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]... - *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) ... - IE]); % state variables - g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables - f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series + f_omega =(1/(2*pi))*([(IA-A*tneg(ig))\ghu1;IE]... + *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % spectral density of state variables; top formula Uhlig (2001), p. 20 with N=0 + g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % spectral density of selected variables; middle formula Uhlig (2001), p. 20; only middle block, i.e. y_t' + f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series; top formula Uhlig (2001), p. 21; end - mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row - % for ifft + mathp_col(ig,:) = (f_hp(:))'; % store as matrix row for ifft end; % Covariance of filtered series - imathp_col = real(ifft(mathp_col))*(2*pi); + imathp_col = real(ifft(mathp_col))*(2*pi); % Inverse Fast Fourier Transformation; middle formula Uhlig (2001), p. 21; Gamma_y{1} = reshape(imathp_col(1,:),nvar,nvar); % Autocorrelations if nar > 0 @@ -253,44 +260,40 @@ else% ==> Theoretical HP filter. Gamma_y{nar+2} = ones(nvar,1); else Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr); - SS(exo_names_orig_ord,exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr); + SS(exo_names_orig_ord,exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr); %make sure Covariance matrix is positive definite cs = chol(SS)'; SS = cs*cs'; b1(:,exo_names_orig_ord) = ghu1; b2(:,exo_names_orig_ord) = ghu(iky,:); - mathp_col = []; + mathp_col = NaN(ngrid,length(ivar)^2); IA = eye(size(A,1)); IE = eye(M_.exo_nbr); for ig = 1:ngrid - if hp1(ig)==0, + if filter_gain(ig)==0, f_hp = zeros(length(ivar),length(ivar)); else - f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]... - *SS*[b1'*inv(IA-A'*tpos(ig)) ... - IE]); % state variables - g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables - f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series + f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\b1;IE]... + *SS*[b1'/(IA-A'*tpos(ig)) IE]); % spectral density of state variables; top formula Uhlig (2001), p. 20 with N=0 + g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % spectral density of selected variables; middle formula Uhlig (2001), p. 20; only middle block, i.e. y_t' + f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series; top formula Uhlig (2001), p. 21; end - mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row - % for ifft + mathp_col(ig,:) = (f_hp(:))'; % store as matrix row for ifft end; imathp_col = real(ifft(mathp_col))*(2*pi); vv = diag(reshape(imathp_col(1,:),nvar,nvar)); for i=1:M_.exo_nbr - mathp_col = []; + mathp_col = NaN(ngrid,length(ivar)^2); SSi = cs(:,i)*cs(:,i)'; for ig = 1:ngrid - if hp1(ig)==0, + if filter_gain(ig)==0, f_hp = zeros(length(ivar),length(ivar)); else - f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]... - *SSi*[b1'*inv(IA-A'*tpos(ig)) ... - IE]); % state variables - g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables - f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series + f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\b1;IE]... + *SSi*[b1'/(IA-A'*tpos(ig)) IE]); % spectral density of state variables; top formula Uhlig (2001), p. 20 with N=0 + g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % spectral density of selected variables; middle formula Uhlig (2001), p. 20; only middle block, i.e. y_t' + f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series; top formula Uhlig (2001), p. 21; end - mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row - % for ifft + mathp_col(ig,:) = (f_hp(:))'; % store as matrix row for ifft end; imathp_col = real(ifft(mathp_col))*(2*pi); Gamma_y{nar+2}(:,i) = abs(diag(reshape(imathp_col(1,:),nvar,nvar)))./vv; From a4b04ca9b4ca2383ca64b1038b57878f4af7db7c Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 4 Aug 2015 08:45:55 +0200 Subject: [PATCH 03/14] Fix UnivariateSpectralDensity.m for filtered variables - Adds bandpass filter - Writes results to oo_ --- matlab/UnivariateSpectralDensity.m | 137 +++++++++++++---------------- matlab/stoch_simul.m | 2 +- 2 files changed, 63 insertions(+), 76 deletions(-) diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m index 95e8dab69..22fee8739 100644 --- a/matlab/UnivariateSpectralDensity.m +++ b/matlab/UnivariateSpectralDensity.m @@ -1,11 +1,25 @@ -function [omega,f] = UnivariateSpectralDensity(dr,var_list) +function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list) % This function computes the theoretical spectral density of each % endogenous variable declared in var_list. Results are stored in -% oo_ and may be plotted. Plots are saved into the graphs-folder. +% oo_.SpectralDensity and may be plotted. Plots are saved into the +% graphs-folder. % +% INPUTS +% M_ [structure] Dynare's model structure +% oo_ [structure] Dynare's results structure +% options_ [structure] Dynare's options structure +% var_list [integer] Vector of indices for a subset of variables. +% +% OUTPUTS +% oo_ [structure] Dynare's results structure, +% containing the subfield +% SpectralDensity with fields freqs +% and density, which are of size nvar*ngrid. +% + % Adapted from th_autocovariances.m. -% Copyright (C) 2006-2013 Dynare Team +% Copyright (C) 2006-2015 Dynare Team % % This file is part of Dynare. % @@ -22,10 +36,6 @@ function [omega,f] = UnivariateSpectralDensity(dr,var_list) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ oo_ M_ - - -omega = []; f = []; if options_.order > 1 disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density') @@ -33,12 +43,7 @@ if options_.order > 1 disp('Please set order = 1. I abort') return end -pltinfo = options_.SpectralDensity.plot; -cutoff = options_.SpectralDensity.cutoff; -sdl = options_.SpectralDensity.sdl; -omega = (0:sdl:pi)'; -GridSize = length(omega); -exo_names_orig_ord = M_.exo_names_orig_ord; + if isoctave warning('off', 'Octave:divide-by-zero') else @@ -60,22 +65,19 @@ for i=1:nvar ivar(i) = i_tmp; end end -f = zeros(nvar,GridSize); -ghx = dr.ghx; -ghu = dr.ghu; + +ghx = oo_.dr.ghx; +ghu = oo_.dr.ghu; nspred = M_.nspred; nstatic = M_.nstatic; -kstate = dr.kstate; -order = dr.order_var; +kstate = oo_.dr.kstate; +order = oo_.dr.order_var; iv(order) = [1:length(order)]; nx = size(ghx,2); ikx = [nstatic+1:nstatic+nspred]; -A = zeros(nx,nx); k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:); i0 = find(k0(:,2) == M_.maximum_lag+1); i00 = i0; -n0 = length(i0); -A(i0,:) = ghx(ikx,:); AS = ghx(:,i0); ghu1 = zeros(nx,M_.exo_nbr); ghu1(i0,:) = ghu(ikx,:); @@ -91,87 +93,72 @@ for i=M_.maximum_lag:-1:2 AS(:,j1) = AS(:,j1)+ghx(:,i1); i0 = i1; end -Gamma = zeros(nvar,cutoff+1); -[A,B] = kalman_transition_matrix(dr,ikx',1:nx,M_.exo_nbr); + +[A,B] = kalman_transition_matrix(oo_.dr,ikx',1:nx,M_.exo_nbr); [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); iky = iv(ivar); if ~isempty(u) iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2))); - ivar = dr.order_var(iky); + ivar = oo_.dr.order_var(iky); end + +iky = iv(ivar); aa = ghx(iky,:); bb = ghu(iky,:); - -if options_.hp_filter == 0 - tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb'; - tmp(abs(tmp) < 1e-12) = 0; - Gamma(:,1) = diag(tmp); - vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb'); - tmp = aa*vxy; - tmp(abs(tmp) < 1e-12) = 0; - Gamma(:,2) = diag(tmp); - for i=2:cutoff - vxy = A*vxy; - tmp = aa*vxy; - tmp(abs(tmp) < 1e-12) = 0; - Gamma(:,i+1) = diag(tmp); - end -else - iky = iv(ivar); - aa = ghx(iky,:); - bb = ghu(iky,:); +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 + filter_gain=ones(ngrid,1); +elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && options_.bandpass.indicator %filter with bandpass + filter_gain = zeros(1,ngrid); + lowest_periodicity=options_.bandpass.passband(2); + highest_periodicity=options_.bandpass.passband(1); + highest_periodicity=max(2,highest_periodicity); % restrict to upper bound of pi + filter_gain(freqs>=2*pi/lowest_periodicity & freqs<=2*pi/highest_periodicity)=1; + filter_gain(freqs<=-2*pi/lowest_periodicity+2*pi & freqs>=-2*pi/highest_periodicity+2*pi)=1; +elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && ~options_.bandpass.indicator %filter with HP-filter lambda = options_.hp_filter; - ngrid = options_.hp_ngrid; - freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); - tpos = exp( sqrt(-1)*freqs); - tneg = exp(-sqrt(-1)*freqs); - hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2); - mathp_col = NaN(ngrid,length(ivar)^2); - IA = eye(size(A,1)); - IE = eye(M_.exo_nbr); - for ig = 1:ngrid - f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\ghu1;IE]... - *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % state variables - g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables - f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series - mathp_col(ig,:) = (f_hp(:))'; % store as matrix row - % for ifft - end; - imathp_col = real(ifft(mathp_col))*(2*pi); - tmp = reshape(imathp_col(1,:),nvar,nvar); - tmp(abs(tmp)<1e-12) = 0; - Gamma(:,1) = diag(tmp); - sy = sqrt(Gamma(:,1)); - sy = sy *sy'; - for i=1:cutoff-1 - tmp = reshape(imathp_col(i+1,:),nvar,nvar)./sy; - tmp(abs(tmp) < 1e-12) = 0; - Gamma(:,i+1) = diag(tmp); - end + filter_gain = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2); end -H = 1:cutoff; +mathp_col = NaN(ngrid,length(ivar)^2); +IA = eye(size(A,1)); +IE = eye(M_.exo_nbr); +for ig = 1:ngrid + f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\ghu1;IE]... + *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % state variables + g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables + f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series + mathp_col(ig,:) = (f_hp(:))'; % store as matrix row +end; + +f = zeros(nvar,ngrid); for i=1:nvar - f(i,:) = Gamma(i,1)/(2*pi) + Gamma(i,H+1)*cos(H'*omega')/pi; + f(i,:) = real(mathp_col(:,(i-1)*nvar+i)); %read out spectral density end +oo_.SpectralDensity.freqs=freqs; +oo_.SpectralDensity.density=f; + if isoctave warning('on', 'Octave:divide-by-zero') else warning on MATLAB:dividebyzero end -if pltinfo +if options_.nograph == 0 if ~exist(M_.fname, 'dir') mkdir('.',M_.fname); end - if ~exist([M_.fname '/graphs']) + if ~exist([M_.fname '/graphs'],'dir') mkdir(M_.fname,'graphs'); end for i= 1:nvar hh = dyn_figure(options_,'Name',['Spectral Density of ' deblank(M_.endo_names(ivar(i),:)) '.']); - plot(omega,f(i,:),'-k','linewidth',2) + plot(freqs,f(i,:),'-k','linewidth',2) xlabel('0 \leq \omega \leq \pi') ylabel('f(\omega)') box on diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index d6b6ceadf..ff8b53cfd 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -350,7 +350,7 @@ if options_.irf end if options_.SpectralDensity.trigger == 1 - [omega,f] = UnivariateSpectralDensity(oo_.dr,var_list); + [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list); end From 9459ff5d8f9ef1a1fffbc4339434e258ee5161db Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 4 Aug 2015 08:47:39 +0200 Subject: [PATCH 04/14] Add bandpass filtering to simulated moments Also takes precautions for future implementation of one-sided hp filter --- matlab/disp_moments.m | 12 ++++++++++-- matlab/global_initialization.m | 2 +- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index cd7ebf2b9..4e6de10d2 100644 --- a/matlab/disp_moments.m +++ b/matlab/disp_moments.m @@ -42,10 +42,18 @@ y = y(ivar,options_.drop+1:end)'; m = mean(y); -if options_.hp_filter +if options_.hp_filter && ~options.one_sided_hp_filter && ~options_.bandpass.indicator [hptrend,y] = sample_hp_filter(y,options_.hp_filter); -else +elseif ~options_.hp_filter && options_.one_sided_hp_filter && ~options_.bandpass.indicator + error('disp_moments:: The one-sided HP filter is not yet available') +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); + y=data_temp.data; +elseif ~options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator y = bsxfun(@minus, y, m); +else + error('disp_moments:: You cannot use more than one filter at the same time') end s2 = mean(y.*y); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 3f65232c9..f093cd2f1 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -153,7 +153,7 @@ options_.impulse_responses.plot_threshold=1e-10; options_.relative_irf = 0; options_.ar = 5; options_.hp_filter = 0; -options_.one_sided_hp_filter = 1600; +options_.one_sided_hp_filter = 0; options_.hp_ngrid = 512; options_.nodecomposition = 0; options_.nomoments = 0; From a524681e407c5722604ec06169b8324a2a9cc1e3 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 7 Aug 2015 16:57:28 +0200 Subject: [PATCH 05/14] Replace globals in disp_moments.m by function arguments --- matlab/disp_moments.m | 16 ++++++++++++---- matlab/stoch_simul.m | 2 +- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index 4e6de10d2..39966d0b4 100644 --- a/matlab/disp_moments.m +++ b/matlab/disp_moments.m @@ -1,7 +1,17 @@ -function disp_moments(y,var_list) +function oo_=disp_moments(y,var_list,M_,options_,oo_) +% function disp_moments(y,var_list,M_,options_,oo_) % Displays moments of simulated variables +% INPUTS +% y [double] nvar*nperiods vector of simulated variables. +% var_list [char] nvar character array with names of variables. +% M_ [structure] Dynare's model structure +% oo_ [structure] Dynare's results structure +% options_ [structure] Dynare's options structure +% +% OUTPUTS +% oo_ [structure] Dynare's results structure, -% Copyright (C) 2001-2012 Dynare Team +% Copyright (C) 2001-2015 Dynare Team % % This file is part of Dynare. % @@ -18,8 +28,6 @@ function disp_moments(y,var_list) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ options_ oo_ - warning_old_state = warning; warning off diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index ff8b53cfd..5f51a33e7 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -160,7 +160,7 @@ if options_.nomoments == 0 disp_th_moments(oo_.dr,var_list); end else - disp_moments(oo_.endo_simul,var_list); + oo_=disp_moments(oo_.endo_simul,var_list,M_,options_,oo_); end end From 281e87f807e2f948d33651cb2a77d3fa43b626cb Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 10 Aug 2015 17:05:56 +0200 Subject: [PATCH 06/14] Add simulated variance decomposition --- matlab/disp_moments.m | 83 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 70 insertions(+), 13 deletions(-) diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index 39966d0b4..06e8413b2 100644 --- a/matlab/disp_moments.m +++ b/matlab/disp_moments.m @@ -50,19 +50,8 @@ y = y(ivar,options_.drop+1:end)'; m = mean(y); -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') -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); - y=data_temp.data; -elseif ~options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator - y = bsxfun(@minus, y, m); -else - error('disp_moments:: You cannot use more than one filter at the same time') -end +% filter series +y=get_filtered_time_series(y,m,options_); s2 = mean(y.*y); s = sqrt(s2); @@ -121,4 +110,72 @@ if ar > 0 end end + +if ~options_.nodecomposition + if M_.exo_nbr == 1 + oo_.variance_decomposition = 100*ones(nvar,1); + else + oo_.variance_decomposition=zeros(nvar,M_.exo_nbr); + %get starting values + if isempty(M_.endo_histval) + y0 = oo_.dr.ys; + else + y0 = M_.endo_histval; + end + %back out shock matrix used for generating y + i_exo_var = setdiff([1:M_.exo_nbr],find(diag(M_.Sigma_e) == 0)); % find shocks with 0 variance + chol_S = chol(M_.Sigma_e(i_exo_var,i_exo_var)); %decompose rest + shock_mat=zeros(options_.periods,M_.exo_nbr); %initialize + shock_mat(:,i_exo_var)=oo_.exo_simul(:,i_exo_var)/chol_S; %invert construction of oo_.exo_simul from simult.m + + for shock_iter=1:length(i_exo_var) + temp_shock_mat=zeros(size(shock_mat)); + temp_shock_mat(:,i_exo_var(shock_iter))=shock_mat(:,i_exo_var(shock_iter)); + temp_shock_mat(:,i_exo_var) = temp_shock_mat(:,i_exo_var)*chol_S; + y_sim_one_shock = simult_(y0,oo_.dr,temp_shock_mat,options_.order); + y_sim_one_shock=y_sim_one_shock(ivar,1+options_.drop+1:end)'; + y_sim_one_shock=get_filtered_time_series(y_sim_one_shock,mean(y_sim_one_shock),options_); + oo_.variance_decomposition(:,i_exo_var(shock_iter))=var(y_sim_one_shock)./s2*100; + end + if ~options_.noprint %options_.nomoments == 0 + skipline() + title='VARIANCE DECOMPOSITION SIMULATING ONE SHOCK AT A TIME (in percent)'; + if options_.hp_filter + title = [title ' (HP filter, lambda = ' ... + num2str(options_.hp_filter) ')']; + end + headers = M_.exo_names; + headers(M_.exo_names_orig_ord,:) = headers; + headers = char(' ',headers); + lh = size(deblank(M_.endo_names(ivar,:)),2)+2; + dyntable(title,char(headers,'Total linear contrib.'),deblank(M_.endo_names(ivar,:)),[oo_.variance_decomposition sum(oo_.variance_decomposition,2)],lh,8,2); + if options_.order == 1 + fprintf('Note: numbers do not add up to 100 due to non-zero correlation of simulated shocks in small samples\n\n') + else + fprintf('Note: numbers do not add up to 100 due to i) non-zero correlation of simulated shocks in small samples and ii) nonlinearity\n\n') + end + end + + end +end + warning(warning_old_state); +end + +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') +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); + y=data_temp.data; +elseif ~options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator + y = bsxfun(@minus, y, m); +else + error('disp_moments:: You cannot use more than one filter at the same time') +end + +end \ No newline at end of file From 3902c8640757f46baf6d555f2600778115c189c1 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 10 Aug 2015 17:11:13 +0200 Subject: [PATCH 07/14] Add unit test for variance decomposition --- tests/Makefile.am | 1 + tests/moments/example1_var_decomp.mod | 89 +++++++++++++++++++++++++++ 2 files changed, 90 insertions(+) create mode 100644 tests/moments/example1_var_decomp.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index 07688dcf0..c8530c417 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -7,6 +7,7 @@ MODFILES = \ estimation/MH_recover/fs2000_recover.mod \ estimation/t_proposal/fs2000_student.mod \ estimation/TaRB/fs2000_tarb.mod \ + moments/example1_var_decomp \ gsa/ls2003.mod \ gsa/ls2003a.mod \ gsa/cod_ML_morris/cod_ML_morris.mod \ diff --git a/tests/moments/example1_var_decomp.mod b/tests/moments/example1_var_decomp.mod new file mode 100644 index 000000000..b02b59735 --- /dev/null +++ b/tests/moments/example1_var_decomp.mod @@ -0,0 +1,89 @@ +/* + * 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, junk; + +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; + +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; + +stoch_simul(relative_irf,order=1,periods=0); +oo_1_theoretic=oo_; +stoch_simul(relative_irf,order=1,periods=100000); +oo_1_simul=oo_; +stoch_simul(relative_irf,order=2,periods=0); +oo_2_theoretic=oo_; +set_dynare_seed('default'); +stoch_simul(relative_irf,order=2,periods=100000); +oo_2_simul=oo_; + +if max(max(abs(oo_1_theoretic.variance_decomposition-oo_1_simul.variance_decomposition)))>2 + error('Variance Decomposition wrong') +end + +if max(max(abs(oo_1_theoretic.variance_decomposition-oo_2_theoretic.variance_decomposition)))>1e-10 + error('Theoretical Variance Decomposition wrong') +end + +if max(max(abs(oo_2_theoretic.variance_decomposition-oo_2_simul.variance_decomposition)))>3 + error('Variance Decomposition wrong') +end + From f7ae5e4f60676dd57c56f3c32bedbf0e0f392c1b Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 10 Aug 2015 17:43:58 +0200 Subject: [PATCH 08/14] Adjust table titles to accommodate new filter options --- matlab/add_filter_subtitle.m | 15 +++++++++++++++ matlab/disp_moments.m | 28 +++++++++++----------------- matlab/disp_th_moments.m | 17 ++++------------- 3 files changed, 30 insertions(+), 30 deletions(-) create mode 100644 matlab/add_filter_subtitle.m diff --git a/matlab/add_filter_subtitle.m b/matlab/add_filter_subtitle.m new file mode 100644 index 000000000..652587517 --- /dev/null +++ b/matlab/add_filter_subtitle.m @@ -0,0 +1,15 @@ +function title=add_filter_subtitle(title,options_) + +if ~options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator %do not filter + %nothing to add here +elseif ~options_.hp_filter && ~options_.one_sided_hp_filter && options_.bandpass.indicator + title = [title ' (Bandpass filter, (' ... + num2str(options_.bandpass.passband(1)),' ',num2str(options_.bandpass.passband(2)), '))']; +elseif options_.hp_filter && ~options_.one_sided_hp_filter && ~options_.bandpass.indicator %filter with HP-filter + title = [title ' (HP filter, lambda = ' ... + num2str(options_.hp_filter) ')']; +elseif ~options_.hp_filter && options_.one_sided_hp_filter && ~options_.bandpass.indicator + title = [title ' (One-sided HP filter, lambda = ' ... + num2str(options_.one_sided_hp_filter) ')']; +end +end \ No newline at end of file diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index 06e8413b2..15f63f9eb 100644 --- a/matlab/disp_moments.m +++ b/matlab/disp_moments.m @@ -63,10 +63,9 @@ labels = deblank(M_.endo_names(ivar,:)); if options_.nomoments == 0 z = [ m' s' s2' (mean(y.^3)./s2.^1.5)' (mean(y.^4)./(s2.*s2)-3)' ]; title='MOMENTS OF SIMULATED VARIABLES'; - if options_.hp_filter - title = [title ' (HP filter, lambda = ' ... - num2str(options_.hp_filter) ')']; - end + + title=add_filter_subtitle(title,options_); + headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE','SKEWNESS', ... 'KURTOSIS'); dyntable(title,headers,labels,z,size(labels,2)+2,16,6); @@ -79,10 +78,9 @@ if options_.nocorr == 0 end if options_.noprint == 0 title = 'CORRELATION OF SIMULATED VARIABLES'; - if options_.hp_filter - title = [title ' (HP filter, lambda = ' ... - num2str(options_.hp_filter) ')']; - end + + title=add_filter_subtitle(title,options_); + headers = char('VARIABLE',M_.endo_names(ivar,:)); dyntable(title,headers,labels,corr,size(labels,2)+2,8,4); end @@ -101,10 +99,7 @@ if ar > 0 end if options_.noprint == 0 title = 'AUTOCORRELATION OF SIMULATED VARIABLES'; - if options_.hp_filter - title = [title ' (HP filter, lambda = ' ... - num2str(options_.hp_filter) ')']; - end + title=add_filter_subtitle(title,options_); headers = char('VARIABLE',int2str([1:ar]')); dyntable(title,headers,labels,autocorr,size(labels,2)+2,8,4); end @@ -140,10 +135,9 @@ if ~options_.nodecomposition if ~options_.noprint %options_.nomoments == 0 skipline() title='VARIANCE DECOMPOSITION SIMULATING ONE SHOCK AT A TIME (in percent)'; - if options_.hp_filter - title = [title ' (HP filter, lambda = ' ... - num2str(options_.hp_filter) ')']; - end + + title=add_filter_subtitle(title,options_); + headers = M_.exo_names; headers(M_.exo_names_orig_ord,:) = headers; headers = char(' ',headers); @@ -164,7 +158,7 @@ end function y=get_filtered_time_series(y,m,options_) -if options_.hp_filter && ~options.one_sided_hp_filter && ~options_.bandpass.indicator +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') diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index b6edede95..2b7cec250 100644 --- a/matlab/disp_th_moments.m +++ b/matlab/disp_th_moments.m @@ -62,9 +62,7 @@ if size(stationary_vars, 1) > 0 else title='THEORETICAL MOMENTS'; end - if options_.hp_filter - title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')']; - end + title=add_filter_subtitle(title,options_); headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE'); labels = deblank(M_.endo_names(ivar,:)); lh = size(labels,2)+2; @@ -77,10 +75,7 @@ if size(stationary_vars, 1) > 0 else title='VARIANCE DECOMPOSITION (in percent)'; end - if options_.hp_filter - title = [title ' (HP filter, lambda = ' ... - num2str(options_.hp_filter) ')']; - end + title=add_filter_subtitle(title,options_); headers = M_.exo_names; headers(M_.exo_names_orig_ord,:) = headers; headers = char(' ',headers); @@ -127,9 +122,7 @@ if options_.nocorr == 0 && size(stationary_vars, 1) > 0 else title='MATRIX OF CORRELATIONS'; end - if options_.hp_filter - title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')']; - end + title=add_filter_subtitle(title,options_); labels = deblank(M_.endo_names(ivar(i1),:)); headers = char('Variables',labels); lh = size(labels,2)+2; @@ -149,9 +142,7 @@ if options_.ar > 0 && size(stationary_vars, 1) > 0 else title='COEFFICIENTS OF AUTOCORRELATION'; end - if options_.hp_filter - title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')']; - end + title=add_filter_subtitle(title,options_); labels = deblank(M_.endo_names(ivar(i1),:)); headers = char('Order ',int2str([1:options_.ar]')); lh = size(labels,2)+2; From 38d293b86447df34da09e74875459eef8537620b Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 17 Aug 2015 08:57:07 +0200 Subject: [PATCH 09/14] Replace globals in disp_th_moments.m by function arguments --- matlab/cli/prior.m | 2 +- matlab/disp_th_moments.m | 6 ++---- matlab/stoch_simul.m | 2 +- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/matlab/cli/prior.m b/matlab/cli/prior.m index d00e0555b..9a4ead849 100644 --- a/matlab/cli/prior.m +++ b/matlab/cli/prior.m @@ -106,7 +106,7 @@ if ismember('moments', varargin) % Prior simulations (2nd order moments). % Solve model [dr,info, M_ ,options_ , oo_] = resol(0, M_ , options_ ,oo_); % Compute and display second order moments - disp_th_moments(oo_.dr,[]); + oo_=disp_th_moments(oo_.dr,[],M_,options_,oo_); skipline(2) donesomething = true; end diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index 2b7cec250..c8bd07b35 100644 --- a/matlab/disp_th_moments.m +++ b/matlab/disp_th_moments.m @@ -1,7 +1,7 @@ -function disp_th_moments(dr,var_list) +function oo_=disp_th_moments(dr,var_list,M_,options_,oo_) % Display theoretical moments of variables -% Copyright (C) 2001-2013 Dynare Team +% Copyright (C) 2001-2015 Dynare Team % % This file is part of Dynare. % @@ -18,8 +18,6 @@ function disp_th_moments(dr,var_list) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ oo_ options_ - nodecomposition = options_.nodecomposition; if size(var_list,1) == 0 diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index 5f51a33e7..6dd07ecc8 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -157,7 +157,7 @@ if options_.nomoments == 0 elseif options_.periods == 0 % There is no code for theoretical moments at 3rd order if options_.order <= 2 - disp_th_moments(oo_.dr,var_list); + oo_=disp_th_moments(oo_.dr,var_list,M_,options_,oo_); end else oo_=disp_moments(oo_.endo_simul,var_list,M_,options_,oo_); From 40877685f2e96f96c634de4750c44a30703df02c Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 17 Aug 2015 13:12:22 +0200 Subject: [PATCH 10/14] Add LaTeX-output option to table from stoch_simul.m --- matlab/disp_moments.m | 25 ++++++- matlab/disp_th_moments.m | 24 ++++++ ...splay_conditional_variance_decomposition.m | 16 +++- matlab/dyn_latex_table.m | 73 +++++++++++++++++++ matlab/stoch_simul.m | 6 ++ tests/TeX/fs2000_corr_ME.mod | 3 +- 6 files changed, 141 insertions(+), 6 deletions(-) create mode 100644 matlab/dyn_latex_table.m diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index 15f63f9eb..96ee22772 100644 --- a/matlab/disp_moments.m +++ b/matlab/disp_moments.m @@ -59,6 +59,7 @@ oo_.mean = transpose(m); oo_.var = y'*y/size(y,1); labels = deblank(M_.endo_names(ivar,:)); +labels_TeX = deblank(M_.endo_names_tex(ivar,:)); if options_.nomoments == 0 z = [ m' s' s2' (mean(y.^3)./s2.^1.5)' (mean(y.^4)./(s2.*s2)-3)' ]; @@ -69,6 +70,9 @@ if options_.nomoments == 0 headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE','SKEWNESS', ... 'KURTOSIS'); dyntable(title,headers,labels,z,size(labels,2)+2,16,6); + if options_.TeX + dyn_latex_table(M_,title,'sim_moments',headers,labels_TeX,z,size(labels,2)+2,16,6); + end end if options_.nocorr == 0 @@ -83,6 +87,11 @@ if options_.nocorr == 0 headers = char('VARIABLE',M_.endo_names(ivar,:)); dyntable(title,headers,labels,corr,size(labels,2)+2,8,4); + if options_.TeX + headers = char('VARIABLE',M_.endo_names_tex(ivar,:)); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'sim_corr_matrix',headers,labels_TeX,corr,size(labels,2)+2,8,4); + end end end @@ -102,7 +111,13 @@ if ar > 0 title=add_filter_subtitle(title,options_); headers = char('VARIABLE',int2str([1:ar]')); dyntable(title,headers,labels,autocorr,size(labels,2)+2,8,4); + if options_.TeX + headers = char('VARIABLE',int2str([1:ar]')); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'sim_autocorr_matrix',headers,labels_TeX,autocorr,size(labels_TeX,2)+2,8,4); + end end + end @@ -142,7 +157,15 @@ if ~options_.nodecomposition headers(M_.exo_names_orig_ord,:) = headers; headers = char(' ',headers); lh = size(deblank(M_.endo_names(ivar,:)),2)+2; - dyntable(title,char(headers,'Total linear contrib.'),deblank(M_.endo_names(ivar,:)),[oo_.variance_decomposition sum(oo_.variance_decomposition,2)],lh,8,2); + dyntable(title,char(headers,'Tot. lin. contr.'),deblank(M_.endo_names(ivar,:)),[oo_.variance_decomposition sum(oo_.variance_decomposition,2)],lh,8,2); + if options_.TeX + headers=M_.exo_names_tex; + headers = char(' ',headers); + labels = deblank(M_.endo_names_tex(ivar,:)); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'sim_var_decomp',char(headers,'Tot. lin. contr.'),labels_TeX,[oo_.variance_decomposition sum(oo_.variance_decomposition,2)],lh,8,2); + end + if options_.order == 1 fprintf('Note: numbers do not add up to 100 due to non-zero correlation of simulated shocks in small samples\n\n') else diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index c8bd07b35..d93a86948 100644 --- a/matlab/disp_th_moments.m +++ b/matlab/disp_th_moments.m @@ -65,6 +65,11 @@ if size(stationary_vars, 1) > 0 labels = deblank(M_.endo_names(ivar,:)); lh = size(labels,2)+2; dyntable(title,headers,labels,z,lh,11,4); + if options_.TeX + labels = deblank(M_.endo_names_tex(ivar,:)); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'th_moments',headers,labels,z,lh,11,4); + end if M_.exo_nbr > 1 && ~nodecomposition skipline() @@ -81,6 +86,13 @@ if size(stationary_vars, 1) > 0 dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ... :)),100* ... oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2); + if options_.TeX + headers=M_.exo_names_tex; + headers = char(' ',headers); + labels = deblank(M_.endo_names_tex(ivar(stationary_vars),:)); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'th_var_decomp_uncond',headers,labels,100*oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2); + end end end @@ -125,6 +137,12 @@ if options_.nocorr == 0 && size(stationary_vars, 1) > 0 headers = char('Variables',labels); lh = size(labels,2)+2; dyntable(title,headers,labels,corr,lh,8,4); + if options_.TeX + labels = deblank(M_.endo_names_tex(ivar(i1),:)); + headers=char('Variables',labels); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'th_corr_matrix',headers,labels,corr,lh,8,4); + end end end if options_.ar > 0 && size(stationary_vars, 1) > 0 @@ -145,5 +163,11 @@ if options_.ar > 0 && size(stationary_vars, 1) > 0 headers = char('Order ',int2str([1:options_.ar]')); lh = size(labels,2)+2; dyntable(title,headers,labels,z,lh,8,4); + if options_.TeX + labels = deblank(M_.endo_names_tex(ivar(i1),:)); + headers=char('Order ',int2str([1:options_.ar]')); + lh = size(labels,2)+2; + dyn_latex_table(M_,title,'th_autocorr_matrix',headers,labels,z,lh,8,4); + end end end diff --git a/matlab/display_conditional_variance_decomposition.m b/matlab/display_conditional_variance_decomposition.m index 2e15a1d9e..944469ef1 100644 --- a/matlab/display_conditional_variance_decomposition.m +++ b/matlab/display_conditional_variance_decomposition.m @@ -31,11 +31,13 @@ function display_conditional_variance_decomposition(conditional_decomposition_ar % along with Dynare. If not, see . if options_.order == 2 -skipline() -disp('APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)') + skipline() + title='APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)'; + disp(title) else -skipline() -disp('CONDITIONAL VARIANCE DECOMPOSITION (in percent)') + skipline() + title='CONDITIONAL VARIANCE DECOMPOSITION (in percent)'; + disp(title) end vardec_i = zeros(length(SubsetOfVariables),M_.exo_nbr); @@ -54,4 +56,10 @@ for i=1:length(Steps) dyntable('',headers,... deblank(M_.endo_names(SubsetOfVariables,:)),... vardec_i,lh,8,2); + if options_.TeX + labels_TeX = deblank(M_.endo_names_tex(SubsetOfVariables,:)); + headers_TeX=char('',deblank(M_.exo_names_tex)); + lh = size(labels_TeX,2)+2; + dyn_latex_table(M_,[title,'; Period ' int2str(Steps(i))],['th_var_decomp_cond_h',int2str(Steps(i))],headers_TeX,labels_TeX,vardec_i,lh,8,2); + end end \ No newline at end of file diff --git a/matlab/dyn_latex_table.m b/matlab/dyn_latex_table.m new file mode 100644 index 000000000..0251179e9 --- /dev/null +++ b/matlab/dyn_latex_table.m @@ -0,0 +1,73 @@ +function dyn_latex_table(M_,title,LaTeXtitle,headers,labels,values,label_width,val_width,val_precis) + +OutputDirectoryName = CheckPath('Output',M_.dname); + +%% get width of label column +if ~isempty(label_width) + label_width = max(size(deblank(char(headers(1,:),labels)),2)+2, ... + label_width); +else %use default length + label_width = max(size(deblank(char(headers(1,:),labels)),2))+2; +end +label_format_leftbound = sprintf('$%%-%ds$',label_width); + +%% get width of label column +if all(~isfinite(values)) + values_length = 4; +else + values_length = max(ceil(max(max(log10(abs(values(isfinite(values))))))),1)+val_precis+1; +end +if any(values) < 0 %add one character for minus sign + values_length = values_length+1; +end + +%% get width of header strings +headers_length = max(size(deblank(headers(2:end,:)),2)); +if ~isempty(val_width) + val_width = max(max(headers_length,values_length)+4,val_width); +else + val_width = max(headers_length,values_length)+4; +end +value_format = sprintf('%%%d.%df',val_width,val_precis); +header_string_format = sprintf('$%%%ds$',val_width); + +%Create and print header string +if length(headers) > 0 + header_string = sprintf(label_format_leftbound ,deblank(headers(1,:))); + header_code_string='l|'; + for i=2:size(headers,1) + header_string = [header_string '\t & \t ' sprintf(header_string_format,strrep(deblank(headers(i,:)),'\','\\'))]; + header_code_string= [header_code_string 'c']; + end +end +header_string=[header_string '\\\\\n']; + +filename = [OutputDirectoryName '/' M_.fname '_' LaTeXtitle '.TeX']; +fidTeX = fopen(filename,'w'); +fprintf(fidTeX,['%% ' datestr(now,0)]); +fprintf(fidTeX,' \n'); +fprintf(fidTeX,' \n'); +fprintf(fidTeX,'\\begin{center}\n'); +fprintf(fidTeX,['\\begin{longtable}{%s} \n'],header_code_string); + +fprintf(fidTeX,['\\caption{',title,'}\\\\\n ']); +fprintf(fidTeX,['\\label{Table:',LaTeXtitle,'}\\\\\n']); +fprintf(fidTeX,'\\hline\\hline \\\\ \n'); +fprintf(fidTeX,header_string); +fprintf(fidTeX,'\\hline \\endfirsthead \n'); +fprintf(fidTeX,'\\caption{(continued)}\\\\\n '); +fprintf(fidTeX,'\\hline\\hline \\\\ \n'); +fprintf(fidTeX,header_string); +fprintf(fidTeX,'\\hline \\endhead \n'); +fprintf(fidTeX,['\\hline \\multicolumn{',num2str(size(headers,1)),'}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n']); +fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n'); +for i=1:size(values,1) + fprintf(fidTeX,label_format_leftbound,deblank(labels(i,:))); + fprintf(fidTeX,['\t & \t' value_format],values(i,:)); + fprintf(fidTeX,' \\\\ \n'); +end + +fprintf(fidTeX,'\\end{longtable}\n '); +fprintf(fidTeX,'\\end{center}\n'); +fprintf(fidTeX,'%% End of TeX file.\n'); +fclose(fidTeX); \ No newline at end of file diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index 6dd07ecc8..750bd4929 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -107,6 +107,12 @@ if ~options_.noprint headers = char('Variables',labels); lh = size(labels,2)+2; dyntable(my_title,headers,labels,M_.Sigma_e,lh,10,6); + if options_.TeX + labels = deblank(M_.exo_names_tex); + headers = char('Variables',labels); + lh = size(labels,2)+2; + dyn_latex_table(M_,my_title,'covar_ex_shocks',headers,labels,M_.Sigma_e,lh,10,6); + end if options_.partial_information skipline() disp('SOLUTION UNDER PARTIAL INFORMATION') diff --git a/tests/TeX/fs2000_corr_ME.mod b/tests/TeX/fs2000_corr_ME.mod index bb19e5fef..dc7a6a277 100644 --- a/tests/TeX/fs2000_corr_ME.mod +++ b/tests/TeX/fs2000_corr_ME.mod @@ -123,7 +123,8 @@ end; steady; -stoch_simul(order=1,irf=20,graph_format=eps,contemporaneous_correlation); +stoch_simul(order=1,irf=20,graph_format=eps,periods=1000,contemporaneous_correlation,conditional_variance_decomposition=[1,3]); +stoch_simul(order=1,irf=20,graph_format=eps,periods=0,contemporaneous_correlation,conditional_variance_decomposition=[1,3]); write_latex_original_model; write_latex_static_model; From 40de49456891d2ba6de9ee90fd4b2baa7c4ae6b2 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 12 Oct 2015 17:13:50 +0200 Subject: [PATCH 11/14] Document new filter option --- doc/dynare.texi | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index eb7ceba22..97780b77b 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -3703,13 +3703,27 @@ Number of points (burnin) dropped at the beginning of simulation before computin @item hp_filter = @var{DOUBLE} Uses HP filter with @math{\lambda} = @var{DOUBLE} before computing -moments. Default: no filter. +moments. If theoretical moments are requested, the spectrum of the model solution is filtered +following the approach outlined in @cite{Uhlig (2001)}. +Default: no filter. @item hp_ngrid = @var{INTEGER} Number of points in the grid for the discrete Inverse Fast Fourier Transform used in the HP filter computation. It may be necessary to increase it for highly autocorrelated processes. Default: @code{512}. +@item bandpass_filter +Uses a bandpass filter with the default passband before computing moments. If theoretical moments are +requested, the spectrum of the model solution is filtered using an ideal bandpass +filter. If empirical moments are requested, the @cite{Baxter and King (1999)}-filter +is used. +Default: no filter. + +@item bandpass_filter = @var{[HIGHEST_PERIODICITY LOWEST_PERIODICITY]} +Uses a bandpass filter before computing moments. The passband is set to a periodicity of @code{HIGHEST_PERIODICITY} +to @code{LOWEST_PERIODICITY}, e.g. 6 to 32 quarters if the model frequency is quarterly. +Default: @code{[6,32]}. + @item irf = @var{INTEGER} @anchor{irf} Number of periods on which to compute the IRFs. Setting @code{irf=0}, @@ -10769,7 +10783,7 @@ ts1 is a dseries object: @deftypefn {dseries} {@var{B} = } baxter_king_filter (@var{A}, @var{hf}, @var{lf}, @var{K}) -Implementation of the Baxter and King (1999) band pass filter for @dseries objects. This filter isolates business cycle fluctuations with a period of length ranging between @var{hf} (high frequency) to @var{lf} (low frequency) using a symmetric moving average smoother with @math{2K+1} points, so that K observations at the beginning and at the end of the sample are lost in the computation of the filter. The default value for @var{hf} is @math{6}, for @var{lf} is @math{32}, and for @var{K} is 12. +Implementation of the @cite{Baxter and King (1999)} band pass filter for @dseries objects. This filter isolates business cycle fluctuations with a period of length ranging between @var{hf} (high frequency) to @var{lf} (low frequency) using a symmetric moving average smoother with @math{2K+1} points, so that K observations at the beginning and at the end of the sample are lost in the computation of the filter. The default value for @var{hf} is @math{6}, for @var{lf} is @math{32}, and for @var{K} is 12. @examplehead @example @@ -12933,6 +12947,11 @@ Backus, David K., Patrick J. Kehoe, and Finn E. Kydland (1992): ``International Real Business Cycles,'' @i{Journal of Political Economy}, 100(4), 745--775 +@item +Baxter, Marianne and Robert G. King (1999): +``Measuring Business Cycles: Approximate Band-pass Filters for Economic Time Series,'' +@i{Review of Economics and Statistics}, 81(4), 575--593 + @item Boucekkine, Raouf (1995): ``An alternative methodology for solving nonlinear forward-looking models,'' @i{Journal of Economic Dynamics @@ -13129,6 +13148,11 @@ Smets, Frank and Rafael Wouters (2003): ``An Estimated Dynamic Stochastic General Equilibrium Model of the Euro Area,'' @i{Journal of the European Economic Association}, 1(5), 1123--1175 +@item +Uhlig, Harald (2001): ``A Toolkit for Analysing Nonlinear Dynamic Stochastic Models Easily,'' +in @i{Computational Methods for the Study of Dynamic +Economies}, Eds. Ramon Marimon and Andrew Scott, Oxford University Press, 30--61 + @item Villemot, Sébastien (2011): ``Solving rational expectations models at first order: what Dynare does,'' @i{Dynare Working Papers}, 2, From 2c63ca88437b49b57a2f7fd9b75f2ac27ec8cdd2 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 12 Oct 2015 20:33:09 +0200 Subject: [PATCH 12/14] Add unit tests for HP filter, bandpass filter, and spectral density --- tests/Makefile.am | 5 +- tests/moments/example1_bp_test.mod | 116 +++++++++++++++ tests/moments/example1_hp_test.mod | 157 ++++++++++++++++++++ tests/moments/test_AR1_spectral_density.mod | 76 ++++++++++ 4 files changed, 353 insertions(+), 1 deletion(-) create mode 100644 tests/moments/example1_bp_test.mod create mode 100644 tests/moments/example1_hp_test.mod create mode 100644 tests/moments/test_AR1_spectral_density.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index c8530c417..ee0d77182 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -7,7 +7,10 @@ MODFILES = \ estimation/MH_recover/fs2000_recover.mod \ estimation/t_proposal/fs2000_student.mod \ estimation/TaRB/fs2000_tarb.mod \ - moments/example1_var_decomp \ + moments/example1_var_decomp.mod \ + moments/example1_hp_test.mod \ + moments/example1_bp_test.mod \ + moments/test_AR1_spectral_density.mod \ gsa/ls2003.mod \ gsa/ls2003a.mod \ gsa/cod_ML_morris/cod_ML_morris.mod \ diff --git a/tests/moments/example1_bp_test.mod b/tests/moments/example1_bp_test.mod new file mode 100644 index 000000000..5dfca6446 --- /dev/null +++ b/tests/moments/example1_bp_test.mod @@ -0,0 +1,116 @@ +/* + * 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); + +stoch_simul(order=1,nofunctions,irf=0,bandpass_filter=[6 32],hp_ngrid=8192); +oo_filtered_all_shocks_theoretical=oo_; + +stoch_simul(order=1,nofunctions,periods=1000000); +oo_filtered_all_shocks_simulated=oo_; + + +shocks; +var e; stderr 0; +var u; stderr 0.009; +var e, u = phi*0.009*0; +end; + +stoch_simul(order=1,nofunctions,irf=0,periods=0); + +oo_filtered_one_shock_theoretical=oo_; + +stoch_simul(order=1,nofunctions,periods=5000000); +oo_filtered_one_shock_simulated=oo_; + + +if max(abs((1-diag(oo_filtered_one_shock_simulated.var)./(diag(oo_filtered_all_shocks_simulated.var)))*100-oo_filtered_all_shocks_theoretical.variance_decomposition(:,1)))>2 + error('Variance Decomposition wrong') +end + +if max(max(abs(oo_filtered_all_shocks_simulated.var-oo_filtered_all_shocks_theoretical.var)))>2e-4; + error('Covariance wrong') +end + +if max(max(abs(oo_filtered_one_shock_simulated.var-oo_filtered_one_shock_theoretical.var)))>1e-4; + error('Covariance wrong') +end + +for ii=1:options_.ar + autocorr_model_all_shocks_simulated(:,ii)=diag(oo_filtered_all_shocks_simulated.autocorr{ii}); + autocorr_model_all_shocks_theoretical(:,ii)=diag(oo_filtered_all_shocks_theoretical.autocorr{ii}); + autocorr_model_one_shock_simulated(:,ii)=diag(oo_filtered_one_shock_simulated.autocorr{ii}); + autocorr_model_one_shock_theoretical(:,ii)=diag(oo_filtered_one_shock_theoretical.autocorr{ii}); +end + +if max(max(abs(autocorr_model_all_shocks_simulated-autocorr_model_all_shocks_theoretical)))>2e-2; + error('Correlation wrong') +end + +if max(max(abs(autocorr_model_one_shock_simulated-autocorr_model_one_shock_theoretical)))>2e-2; + error('Correlation wrong') +end diff --git a/tests/moments/example1_hp_test.mod b/tests/moments/example1_hp_test.mod new file mode 100644 index 000000000..c74adeab7 --- /dev/null +++ b/tests/moments/example1_hp_test.mod @@ -0,0 +1,157 @@ +/* + * 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; +options_.bandpass.indicator=0; +options_.bandpass.passband=[6 32]; + +stoch_simul(order=1,nofunctions,hp_filter=1600,irf=0); + +total_var_filtered=diag(oo_.var); +oo_filtered_all_shocks=oo_; + +stoch_simul(order=1,nofunctions,hp_filter=0,periods=5000000,nomoments); +options_.nomoments=0; +oo_unfiltered_all_shocks=oo_; + +[junk, y_filtered]=sample_hp_filter(y,1600); +[junk, c_filtered]=sample_hp_filter(c,1600); +[junk, k_filtered]=sample_hp_filter(k,1600); +[junk, a_filtered]=sample_hp_filter(a,1600); +[junk, h_filtered]=sample_hp_filter(h,1600); +[junk, b_filtered]=sample_hp_filter(b,1600); + +verbatim; +total_std_all_shocks_filtered_sim=std([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) +cov_filtered_all_shocks=cov([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) +acf(1,:)=autocorr([y_filtered ],5)'; +acf(2,:)=autocorr([c_filtered ],5)'; +acf(3,:)=autocorr([k_filtered ],5)'; +acf(4,:)=autocorr([a_filtered ],5)'; +acf(5,:)=autocorr([h_filtered ],5)'; +acf(6,:)=autocorr([b_filtered ],5)'; +autocorr_filtered_all_shocks=acf(:,2:end); +end; + +shocks; +var e; stderr 0; +var u; stderr 0.009; +var e, u = phi*0.009*0; +end; + +stoch_simul(order=1,nofunctions,hp_filter=1600,irf=0,periods=0); + +total_var_filtered_one_shock=diag(oo_.var); +oo_filtered_one_shock=oo_; + +stoch_simul(order=1,nofunctions,hp_filter=0,periods=5000000,nomoments); +oo_unfiltered_one_shock=oo_; + +[junk, y_filtered]=sample_hp_filter(y,1600); +[junk, c_filtered]=sample_hp_filter(c,1600); +[junk, k_filtered]=sample_hp_filter(k,1600); +[junk, a_filtered]=sample_hp_filter(a,1600); +[junk, h_filtered]=sample_hp_filter(h,1600); +[junk, b_filtered]=sample_hp_filter(b,1600); + +verbatim; +total_std_one_shock_filtered_sim=std([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) +cov_filtered_one_shock=cov([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) +acf(1,:)=autocorr([y_filtered ],5)'; +acf(2,:)=autocorr([c_filtered ],5)'; +acf(3,:)=autocorr([k_filtered ],5)'; +acf(4,:)=autocorr([a_filtered ],5)'; +acf(5,:)=autocorr([h_filtered ],5)'; +acf(6,:)=autocorr([b_filtered ],5)'; +autocorr_filtered_one_shock=acf(:,2:end); +end; + +if max(abs((1-(total_std_one_shock_filtered_sim.^2)./(total_std_all_shocks_filtered_sim.^2))*100-oo_filtered_all_shocks.variance_decomposition(:,1)'))>2 + error('Variance Decomposition wrong') +end + +if max(max(abs(oo_filtered_all_shocks.var-cov_filtered_all_shocks)))>1e-4; + error('Covariance wrong') +end + +if max(max(abs(oo_filtered_one_shock.var-cov_filtered_one_shock)))>5e-5; + error('Covariance wrong') +end + +for ii=1:options_.ar + autocorr_model_all_shocks(:,ii)=diag(oo_filtered_all_shocks.autocorr{ii}); + autocorr_model_one_shock(:,ii)=diag(oo_filtered_one_shock.autocorr{ii}); +end + +if max(max(abs(autocorr_model_all_shocks-autocorr_filtered_all_shocks)))>1e-2; + error('Covariance wrong') +end + +if max(max(abs(autocorr_model_one_shock-autocorr_filtered_one_shock)))>1e-2; + error('Covariance wrong') +end diff --git a/tests/moments/test_AR1_spectral_density.mod b/tests/moments/test_AR1_spectral_density.mod new file mode 100644 index 000000000..334776b98 --- /dev/null +++ b/tests/moments/test_AR1_spectral_density.mod @@ -0,0 +1,76 @@ +var white_noise ar1; +varexo e; + +parameters phi; + +phi=0.9; + +model; +white_noise=e; +ar1=phi*ar1(-1)+e; + +end; + +shocks; +var e = 1; +end; + +options_.SpectralDensity.trigger=1; + +options_.bandpass.indicator=0; +options_.hp_ngrid=2048; + +stoch_simul(order=1,nofunctions,hp_filter=0,irf=0,periods=1000000); + +white_noise_sample=white_noise; + +theoretical_spectrum_white_noise=1^2/(2*pi); %Hamilton (1994), 6.1.9 +if max(abs(oo_.SpectralDensity.density(strmatch('white_noise',M_.endo_names,'exact'),:)-theoretical_spectrum_white_noise))>1e-10 + error('Spectral Density is wrong') +end + +theoretical_spectrum_AR1=1/(2*pi)*(1^2./(1+phi^2-2*phi*cos(oo_.SpectralDensity.freqs))); %Hamilton (1994), 6.1.13 +if max(abs(oo_.SpectralDensity.density(strmatch('ar1',M_.endo_names,'exact'),:)-theoretical_spectrum_AR1'))>1e-10 + error('Spectral Density is wrong') +end + +stoch_simul(order=1,nofunctions,hp_filter=1600,irf=0,periods=0); +lambda=options_.hp_filter; +Kalman_gain=(4*lambda*(1 - cos(oo_.SpectralDensity.freqs)).^2 ./ (1 + 4*lambda*(1 - cos(oo_.SpectralDensity.freqs)).^2)); +theoretical_spectrum_white_noise_hp_filtered=1^2/(2*pi)*Kalman_gain.^2; %Hamilton (1994), 6.1.9 +if max(abs(oo_.SpectralDensity.density(strmatch('white_noise',M_.endo_names,'exact'),:)-theoretical_spectrum_white_noise_hp_filtered'))>1e-10 + error('Spectral Density is wrong') +end + +theoretical_spectrum_AR1_hp_filtered=1/(2*pi)*(1^2./(1+phi^2-2*phi*cos(oo_.SpectralDensity.freqs))).*Kalman_gain.^2; %Hamilton (1994), 6.1.13 +if max(abs(oo_.SpectralDensity.density(strmatch('ar1',M_.endo_names,'exact'),:)-theoretical_spectrum_AR1_hp_filtered'))>1e-10 + error('Spectral Density is wrong') +end + +options_.hp_filter=0; +stoch_simul(order=1,nofunctions,bandpass_filter=[6 32],irf=0); + +theoretical_spectrum_white_noise=repmat(theoretical_spectrum_white_noise,1,options_.hp_ngrid); +passband=oo_.SpectralDensity.freqs>=2*pi/options_.bandpass.passband(2) & oo_.SpectralDensity.freqs<=2*pi/options_.bandpass.passband(1); +if max(abs(oo_.SpectralDensity.density(strmatch('white_noise',M_.endo_names,'exact'),passband)-theoretical_spectrum_white_noise(passband)))>1e-10 + error('Spectral Density is wrong') +end +if max(abs(oo_.SpectralDensity.density(strmatch('white_noise',M_.endo_names,'exact'),~passband)-0))>1e-10 + error('Spectral Density is wrong') +end + +if max(abs(oo_.SpectralDensity.density(strmatch('ar1',M_.endo_names,'exact'),passband)-theoretical_spectrum_AR1(passband)'))>1e-10 + error('Spectral Density is wrong') +end +if max(abs(oo_.SpectralDensity.density(strmatch('ar1',M_.endo_names,'exact'),~passband)-0))>1e-10 + error('Spectral Density is wrong') +end + + +% [pow,f]=psd(a_sample,1024,1,[],512); +% figure +% plot(f,pow/(2*pi)) +% +% % figure +% % [pow,f]=psd(a_sample,1000,1,[],500); +% % plot(f(3:end)*2*pi,pow(3:end)/(2*pi)); \ No newline at end of file From f7175745acfb5d07e4ffd7b44f30f476348f92f5 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 13 Oct 2015 20:20:33 +0200 Subject: [PATCH 13/14] Add one-sided HP filter --- license.txt | 7 +- matlab/UnivariateSpectralDensity.m | 4 +- matlab/disp_moments.m | 2 +- matlab/disp_th_moments.m | 4 +- matlab/one_sided_hp_filter.m | 113 +++++++++++++++++++ tests/Makefile.am | 1 + tests/moments/example1_one_sided_hp_test.mod | 71 ++++++++++++ 7 files changed, 198 insertions(+), 4 deletions(-) create mode 100644 matlab/one_sided_hp_filter.m create mode 100644 tests/moments/example1_one_sided_hp_test.mod 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 From a13c852feb7af0a84ffbad0d01cb69349a2866cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 13 Oct 2015 23:50:38 +0200 Subject: [PATCH 14/14] Removed calls to autocorr routine in integration test. This routine is part of the matlab's statistical toolbox which is not required by Dynare (and not available on karaba). Use sample_autocovariance routine instead. --- tests/moments/example1_hp_test.mod | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/tests/moments/example1_hp_test.mod b/tests/moments/example1_hp_test.mod index c74adeab7..718db4417 100644 --- a/tests/moments/example1_hp_test.mod +++ b/tests/moments/example1_hp_test.mod @@ -89,13 +89,14 @@ oo_unfiltered_all_shocks=oo_; verbatim; total_std_all_shocks_filtered_sim=std([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) cov_filtered_all_shocks=cov([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) -acf(1,:)=autocorr([y_filtered ],5)'; -acf(2,:)=autocorr([c_filtered ],5)'; -acf(3,:)=autocorr([k_filtered ],5)'; -acf(4,:)=autocorr([a_filtered ],5)'; -acf(5,:)=autocorr([h_filtered ],5)'; -acf(6,:)=autocorr([b_filtered ],5)'; -autocorr_filtered_all_shocks=acf(:,2:end); +acf = zeros(6); +[junk, acf(:,1)] = sample_autocovariance([y_filtered ],5); +[junk, acf(:,2)] = sample_autocovariance([c_filtered ],5); +[junk, acf(:,3)] = sample_autocovariance([k_filtered ],5); +[junk, acf(:,4)] = sample_autocovariance([a_filtered ],5); +[junk, acf(:,5)] = sample_autocovariance([h_filtered ],5); +[junk, acf(:,6)] = sample_autocovariance([b_filtered ],5); +autocorr_filtered_all_shocks=acf(2:end,:)'; end; shocks; @@ -122,13 +123,14 @@ oo_unfiltered_one_shock=oo_; verbatim; total_std_one_shock_filtered_sim=std([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) cov_filtered_one_shock=cov([y_filtered c_filtered k_filtered a_filtered h_filtered b_filtered]) -acf(1,:)=autocorr([y_filtered ],5)'; -acf(2,:)=autocorr([c_filtered ],5)'; -acf(3,:)=autocorr([k_filtered ],5)'; -acf(4,:)=autocorr([a_filtered ],5)'; -acf(5,:)=autocorr([h_filtered ],5)'; -acf(6,:)=autocorr([b_filtered ],5)'; -autocorr_filtered_one_shock=acf(:,2:end); +acf = zeros(6); +[junk, acf(:,1)] = sample_autocovariance([y_filtered ],5); +[junk, acf(:,2)] = sample_autocovariance([c_filtered ],5); +[junk, acf(:,3)] = sample_autocovariance([k_filtered ],5); +[junk, acf(:,4)] = sample_autocovariance([a_filtered ],5); +[junk, acf(:,5)] = sample_autocovariance([h_filtered ],5); +[junk, acf(:,6)] = sample_autocovariance([b_filtered ],5); +autocorr_filtered_one_shock=acf(2:end,:)'; end; if max(abs((1-(total_std_one_shock_filtered_sim.^2)./(total_std_all_shocks_filtered_sim.^2))*100-oo_filtered_all_shocks.variance_decomposition(:,1)'))>2