diff --git a/doc/dynare.texi b/doc/dynare.texi index 48f3f3ed2..a7a68c01a 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}, @@ -10774,7 +10788,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 @@ -12938,6 +12952,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 @@ -13138,6 +13157,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, 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 9aae99c4c..7a69585cb 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,92 +93,74 @@ 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'; - k = find(abs(tmp) < 1e-12); - tmp(k) = 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; - Gamma(:,2) = diag(tmp); - for i=2:cutoff - vxy = A*vxy; - tmp = aa*vxy; - k = find(abs(tmp) < 1e-12); - tmp(k) = 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_.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); + 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 = []; - 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 - 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 - % 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; - 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; - 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/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/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_moments.m b/matlab/disp_moments.m index cd7ebf2b9..e6a9f3bcc 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 @@ -42,11 +50,8 @@ y = y(ivar,options_.drop+1:end)'; m = mean(y); -if options_.hp_filter - [hptrend,y] = sample_hp_filter(y,options_.hp_filter); -else - y = bsxfun(@minus, y, m); -end +% filter series +y=get_filtered_time_series(y,m,options_); s2 = mean(y.*y); s = sqrt(s2); @@ -54,17 +59,20 @@ 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)' ]; 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); + 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 @@ -74,12 +82,16 @@ 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); + 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 @@ -96,13 +108,91 @@ 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); + 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 + +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)'; + + title=add_filter_subtitle(title,options_); + + 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,'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 + 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 + [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); + 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 diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index b6edede95..35406a01f 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,10 +18,10 @@ 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 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 @@ -62,13 +62,16 @@ 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; 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() @@ -77,10 +80,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); @@ -88,6 +88,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 @@ -127,13 +134,17 @@ 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; 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 @@ -149,12 +160,16 @@ 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; 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/global_initialization.m b/matlab/global_initialization.m index 40c3ff858..f4ac9dbaa 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; 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/matlab/stoch_simul.m b/matlab/stoch_simul.m index d6b6ceadf..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') @@ -157,10 +163,10 @@ 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 - disp_moments(oo_.endo_simul,var_list); + oo_=disp_moments(oo_.endo_simul,var_list,M_,options_,oo_); end end @@ -350,7 +356,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 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; diff --git a/tests/Makefile.am b/tests/Makefile.am index 706ec8907..ef4df9c81 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -7,6 +7,11 @@ MODFILES = \ estimation/MH_recover/fs2000_recover.mod \ estimation/t_proposal/fs2000_student.mod \ estimation/TaRB/fs2000_tarb.mod \ + moments/example1_var_decomp.mod \ + 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/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; 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..718db4417 --- /dev/null +++ b/tests/moments/example1_hp_test.mod @@ -0,0 +1,159 @@ +/* + * 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 = 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; +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 = 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 + 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/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 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 + 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