stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'johannes_bandpass'

time-shift
Stéphane Adjemian (Charybdis) 2015-10-13 23:56:22 +02:00
parent 8c44936846 a13c852feb
commit 6ea5bdde34
20 changed files with 1025 additions and 172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
doc/dynare.texi
View File

@ -3703,13 +3703,27 @@ Number of points (burnin) dropped at the beginning of simulation before computin
@item hp_filter = @var{DOUBLE} @item hp_filter = @var{DOUBLE}
Uses HP filter with @math{\lambda} = @var{DOUBLE} before computing 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} @item hp_ngrid = @var{INTEGER}
Number of points in the grid for the discrete Inverse Fast Fourier Number of points in the grid for the discrete Inverse Fast Fourier
Transform used in the HP filter computation. It may be necessary to Transform used in the HP filter computation. It may be necessary to
increase it for highly autocorrelated processes. Default: @code{512}. 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} @item irf = @var{INTEGER}
@anchor{irf} @anchor{irf}
Number of periods on which to compute the IRFs. Setting @code{irf=0}, 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}) @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 @examplehead
@example @example
@ -12938,6 +12952,11 @@ Backus, David K., Patrick J. Kehoe, and Finn E. Kydland (1992):
``International Real Business Cycles,'' @i{Journal of Political ``International Real Business Cycles,'' @i{Journal of Political
Economy}, 100(4), 745--775 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 @item
Boucekkine, Raouf (1995): ``An alternative methodology for solving Boucekkine, Raouf (1995): ``An alternative methodology for solving
nonlinear forward-looking models,'' @i{Journal of Economic Dynamics 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 Stochastic General Equilibrium Model of the Euro Area,'' @i{Journal of
the European Economic Association}, 1(5), 1123--1175 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 @item
Villemot, Sébastien (2011): ``Solving rational expectations models at Villemot, Sébastien (2011): ``Solving rational expectations models at
first order: what Dynare does,'' @i{Dynare Working Papers}, 2, first order: what Dynare does,'' @i{Dynare Working Papers}, 2,

7
license.txt
View File

@ -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-Name: Dynare
Upstream-Contact: Dynare Team, whose members in 2015 are: Upstream-Contact: Dynare Team, whose members in 2015 are:
Stéphane Adjemian <stephane.adjemian@univ-lemans.fr> Stéphane Adjemian <stephane.adjemian@univ-lemans.fr>
@ -76,6 +76,11 @@ Copyright: 2008-2012 VZLU Prague, a.s.
2014 Dynare Team 2014 Dynare Team
License: GPL-3+ 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 Files: matlab/optimization/simpsa.m matlab/optimization/simpsaget.m matlab/optimization/simpsaset.m
Copyright: 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se Copyright: 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se
2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be

144
matlab/UnivariateSpectralDensity.m
View File

@ -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 % This function computes the theoretical spectral density of each
% endogenous variable declared in var_list. Results are stored in % 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. % Adapted from th_autocovariances.m.
% Copyright (C) 2006-2013 Dynare Team % Copyright (C) 2006-2015 Dynare Team
% %
% This file is part of Dynare. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ oo_ M_
omega = []; f = [];
if options_.order > 1 if options_.order > 1
disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density') disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density')
@ -33,12 +43,7 @@ if options_.order > 1
disp('Please set order = 1. I abort') disp('Please set order = 1. I abort')
return return
end 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 if isoctave
warning('off', 'Octave:divide-by-zero') warning('off', 'Octave:divide-by-zero')
else else
@ -60,22 +65,19 @@ for i=1:nvar
ivar(i) = i_tmp; ivar(i) = i_tmp;
end end
end end
f = zeros(nvar,GridSize);
ghx = dr.ghx; ghx = oo_.dr.ghx;
ghu = dr.ghu; ghu = oo_.dr.ghu;
nspred = M_.nspred; nspred = M_.nspred;
nstatic = M_.nstatic; nstatic = M_.nstatic;
kstate = dr.kstate; kstate = oo_.dr.kstate;
order = dr.order_var; order = oo_.dr.order_var;
iv(order) = [1:length(order)]; iv(order) = [1:length(order)];
nx = size(ghx,2); nx = size(ghx,2);
ikx = [nstatic+1:nstatic+nspred]; ikx = [nstatic+1:nstatic+nspred];
A = zeros(nx,nx);
k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:); k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
i0 = find(k0(:,2) == M_.maximum_lag+1); i0 = find(k0(:,2) == M_.maximum_lag+1);
i00 = i0; i00 = i0;
n0 = length(i0);
A(i0,:) = ghx(ikx,:);
AS = ghx(:,i0); AS = ghx(:,i0);
ghu1 = zeros(nx,M_.exo_nbr); ghu1 = zeros(nx,M_.exo_nbr);
ghu1(i0,:) = ghu(ikx,:); ghu1(i0,:) = ghu(ikx,:);
@ -91,92 +93,74 @@ for i=M_.maximum_lag:-1:2
AS(:,j1) = AS(:,j1)+ghx(:,i1); AS(:,j1) = AS(:,j1)+ghx(:,i1);
i0 = i1; i0 = i1;
end 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); [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); iky = iv(ivar);
if ~isempty(u) if ~isempty(u)
iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2))); 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 end
iky = iv(ivar);
aa = ghx(iky,:); aa = ghx(iky,:);
bb = ghu(iky,:); bb = ghu(iky,:);
ngrid = options_.hp_ngrid; %number of grid points
if options_.hp_filter == 0 freqs = (0 : pi/(ngrid-1):pi)'; % grid on which to compute
tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb'; tpos = exp( sqrt(-1)*freqs); %positive frequencies
k = find(abs(tmp) < 1e-12); tneg = exp(-sqrt(-1)*freqs); %negative frequencies
tmp(k) = 0; if options_.one_sided_hp_filter
Gamma(:,1) = diag(tmp); error('UnivariateSpectralDensity:: spectral density estimate not available with one-sided HP filter')
vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb'); elseif options_.hp_filter == 0 && ~options_.bandpass.indicator %do not filter
tmp = aa*vxy; filter_gain=ones(ngrid,1);
k = find(abs(tmp) < 1e-12); elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && options_.bandpass.indicator %filter with bandpass
tmp(k) = 0; filter_gain = zeros(1,ngrid);
Gamma(:,2) = diag(tmp); lowest_periodicity=options_.bandpass.passband(2);
for i=2:cutoff highest_periodicity=options_.bandpass.passband(1);
vxy = A*vxy; highest_periodicity=max(2,highest_periodicity); % restrict to upper bound of pi
tmp = aa*vxy; filter_gain(freqs>=2*pi/lowest_periodicity & freqs<=2*pi/highest_periodicity)=1;
k = find(abs(tmp) < 1e-12); filter_gain(freqs<=-2*pi/lowest_periodicity+2*pi & freqs>=-2*pi/highest_periodicity+2*pi)=1;
tmp(k) = 0; elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && ~options_.bandpass.indicator %filter with HP-filter
Gamma(:,i+1) = diag(tmp);
end
else
iky = iv(ivar);
aa = ghx(iky,:);
bb = ghu(iky,:);
lambda = options_.hp_filter; lambda = options_.hp_filter;
ngrid = options_.hp_ngrid; filter_gain = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
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
end 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 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 end
oo_.SpectralDensity.freqs=freqs;
oo_.SpectralDensity.density=f;
if isoctave if isoctave
warning('on', 'Octave:divide-by-zero') warning('on', 'Octave:divide-by-zero')
else else
warning on MATLAB:dividebyzero warning on MATLAB:dividebyzero
end end
if pltinfo if options_.nograph == 0
if ~exist(M_.fname, 'dir') if ~exist(M_.fname, 'dir')
mkdir('.',M_.fname); mkdir('.',M_.fname);
end end
if ~exist([M_.fname '/graphs']) if ~exist([M_.fname '/graphs'],'dir')
mkdir(M_.fname,'graphs'); mkdir(M_.fname,'graphs');
end end
for i= 1:nvar for i= 1:nvar
hh = dyn_figure(options_,'Name',['Spectral Density of ' deblank(M_.endo_names(ivar(i),:)) '.']); 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') xlabel('0 \leq \omega \leq \pi')
ylabel('f(\omega)') ylabel('f(\omega)')
box on box on

15
matlab/add_filter_subtitle.m Normal file
View File

@ -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

2
matlab/cli/prior.m
View File

@ -106,7 +106,7 @@ if ismember('moments', varargin) % Prior simulations (2nd order moments).
% Solve model % Solve model
[dr,info, M_ ,options_ , oo_] = resol(0, M_ , options_ ,oo_); [dr,info, M_ ,options_ , oo_] = resol(0, M_ , options_ ,oo_);
% Compute and display second order moments % Compute and display second order moments
disp_th_moments(oo_.dr,[]); oo_=disp_th_moments(oo_.dr,[],M_,options_,oo_);
skipline(2) skipline(2)
donesomething = true; donesomething = true;
end end

132
matlab/disp_moments.m
View File

@ -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 % 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. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ oo_
warning_old_state = warning; warning_old_state = warning;
warning off warning off
@ -42,11 +50,8 @@ y = y(ivar,options_.drop+1:end)';
m = mean(y); m = mean(y);
if options_.hp_filter % filter series
[hptrend,y] = sample_hp_filter(y,options_.hp_filter); y=get_filtered_time_series(y,m,options_);
else
y = bsxfun(@minus, y, m);
end
s2 = mean(y.*y); s2 = mean(y.*y);
s = sqrt(s2); s = sqrt(s2);
@ -54,17 +59,20 @@ oo_.mean = transpose(m);
oo_.var = y'*y/size(y,1); oo_.var = y'*y/size(y,1);
labels = deblank(M_.endo_names(ivar,:)); labels = deblank(M_.endo_names(ivar,:));
labels_TeX = deblank(M_.endo_names_tex(ivar,:));
if options_.nomoments == 0 if options_.nomoments == 0
z = [ m' s' s2' (mean(y.^3)./s2.^1.5)' (mean(y.^4)./(s2.*s2)-3)' ]; z = [ m' s' s2' (mean(y.^3)./s2.^1.5)' (mean(y.^4)./(s2.*s2)-3)' ];
title='MOMENTS OF SIMULATED VARIABLES'; title='MOMENTS OF SIMULATED VARIABLES';
if options_.hp_filter
title = [title ' (HP filter, lambda = ' ... title=add_filter_subtitle(title,options_);
num2str(options_.hp_filter) ')'];
end
headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE','SKEWNESS', ... headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE','SKEWNESS', ...
'KURTOSIS'); 'KURTOSIS');
dyntable(title,headers,labels,z,size(labels,2)+2,16,6); 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 end
if options_.nocorr == 0 if options_.nocorr == 0
@ -74,12 +82,16 @@ if options_.nocorr == 0
end end
if options_.noprint == 0 if options_.noprint == 0
title = 'CORRELATION OF SIMULATED VARIABLES'; title = 'CORRELATION OF SIMULATED VARIABLES';
if options_.hp_filter
title = [title ' (HP filter, lambda = ' ... title=add_filter_subtitle(title,options_);
num2str(options_.hp_filter) ')'];
end
headers = char('VARIABLE',M_.endo_names(ivar,:)); headers = char('VARIABLE',M_.endo_names(ivar,:));
dyntable(title,headers,labels,corr,size(labels,2)+2,8,4); 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
end end
@ -96,13 +108,91 @@ if ar > 0
end end
if options_.noprint == 0 if options_.noprint == 0
title = 'AUTOCORRELATION OF SIMULATED VARIABLES'; title = 'AUTOCORRELATION OF SIMULATED VARIABLES';
if options_.hp_filter title=add_filter_subtitle(title,options_);
title = [title ' (HP filter, lambda = ' ...
num2str(options_.hp_filter) ')'];
end
headers = char('VARIABLE',int2str([1:ar]')); headers = char('VARIABLE',int2str([1:ar]'));
dyntable(title,headers,labels,autocorr,size(labels,2)+2,8,4); 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
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); 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

51
matlab/disp_th_moments.m
View File

@ -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 % Display theoretical moments of variables
% Copyright (C) 2001-2013 Dynare Team % Copyright (C) 2001-2015 Dynare Team
% %
% This file is part of Dynare. % 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 % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_ options_
nodecomposition = options_.nodecomposition; 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 if size(var_list,1) == 0
var_list = M_.endo_names(1:M_.orig_endo_nbr, :); var_list = M_.endo_names(1:M_.orig_endo_nbr, :);
end end
@ -62,13 +62,16 @@ if size(stationary_vars, 1) > 0
else else
title='THEORETICAL MOMENTS'; title='THEORETICAL MOMENTS';
end end
if options_.hp_filter title=add_filter_subtitle(title,options_);
title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];
end
headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE'); headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE');
labels = deblank(M_.endo_names(ivar,:)); labels = deblank(M_.endo_names(ivar,:));
lh = size(labels,2)+2; lh = size(labels,2)+2;
dyntable(title,headers,labels,z,lh,11,4); 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 if M_.exo_nbr > 1 && ~nodecomposition
skipline() skipline()
@ -77,10 +80,7 @@ if size(stationary_vars, 1) > 0
else else
title='VARIANCE DECOMPOSITION (in percent)'; title='VARIANCE DECOMPOSITION (in percent)';
end end
if options_.hp_filter title=add_filter_subtitle(title,options_);
title = [title ' (HP filter, lambda = ' ...
num2str(options_.hp_filter) ')'];
end
headers = M_.exo_names; headers = M_.exo_names;
headers(M_.exo_names_orig_ord,:) = headers; headers(M_.exo_names_orig_ord,:) = headers;
headers = char(' ',headers); headers = char(' ',headers);
@ -88,6 +88,13 @@ if size(stationary_vars, 1) > 0
dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ... dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ...
:)),100* ... :)),100* ...
oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2); 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
end end
@ -127,13 +134,17 @@ if options_.nocorr == 0 && size(stationary_vars, 1) > 0
else else
title='MATRIX OF CORRELATIONS'; title='MATRIX OF CORRELATIONS';
end end
if options_.hp_filter title=add_filter_subtitle(title,options_);
title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];
end
labels = deblank(M_.endo_names(ivar(i1),:)); labels = deblank(M_.endo_names(ivar(i1),:));
headers = char('Variables',labels); headers = char('Variables',labels);
lh = size(labels,2)+2; lh = size(labels,2)+2;
dyntable(title,headers,labels,corr,lh,8,4); 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
end end
if options_.ar > 0 && size(stationary_vars, 1) > 0 if options_.ar > 0 && size(stationary_vars, 1) > 0
@ -149,12 +160,16 @@ if options_.ar > 0 && size(stationary_vars, 1) > 0
else else
title='COEFFICIENTS OF AUTOCORRELATION'; title='COEFFICIENTS OF AUTOCORRELATION';
end end
if options_.hp_filter title=add_filter_subtitle(title,options_);
title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];
end
labels = deblank(M_.endo_names(ivar(i1),:)); labels = deblank(M_.endo_names(ivar(i1),:));
headers = char('Order ',int2str([1:options_.ar]')); headers = char('Order ',int2str([1:options_.ar]'));
lh = size(labels,2)+2; lh = size(labels,2)+2;
dyntable(title,headers,labels,z,lh,8,4); 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
end end

16
matlab/display_conditional_variance_decomposition.m
View File

@ -31,11 +31,13 @@ function display_conditional_variance_decomposition(conditional_decomposition_ar
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if options_.order == 2 if options_.order == 2
skipline() skipline()
disp('APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)') title='APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)';
disp(title)
else else
skipline() skipline()
disp('CONDITIONAL VARIANCE DECOMPOSITION (in percent)') title='CONDITIONAL VARIANCE DECOMPOSITION (in percent)';
disp(title)
end end
vardec_i = zeros(length(SubsetOfVariables),M_.exo_nbr); vardec_i = zeros(length(SubsetOfVariables),M_.exo_nbr);
@ -54,4 +56,10 @@ for i=1:length(Steps)
dyntable('',headers,... dyntable('',headers,...
deblank(M_.endo_names(SubsetOfVariables,:)),... deblank(M_.endo_names(SubsetOfVariables,:)),...
vardec_i,lh,8,2); 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 end

73
matlab/dyn_latex_table.m Normal file
View File

@ -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);

2
matlab/global_initialization.m
View File

@ -153,7 +153,7 @@ options_.impulse_responses.plot_threshold=1e-10;
options_.relative_irf = 0; options_.relative_irf = 0;
options_.ar = 5; options_.ar = 5;
options_.hp_filter = 0; options_.hp_filter = 0;
options_.one_sided_hp_filter = 1600; options_.one_sided_hp_filter = 0;
options_.hp_ngrid = 512; options_.hp_ngrid = 512;
options_.nodecomposition = 0; options_.nodecomposition = 0;
options_.nomoments = 0; options_.nomoments = 0;

113
matlab/one_sided_hp_filter.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
if nargin < 2 || isempty(lambda)
lambda = 1600; %If the user didn't provide a value for lambda, set it to the default value 1600
end
[T,n] = size (y);% Calculate the number of periods and the number of variables in the series
%Set up state space
q=1/lambda; % the signal-to-noise ration: i.e. var eta_t / var epsilon_t
F=[2,-1;
1,0]; % state transition matrix
H=[1,0]; % observation matrix
Q=[q,0;
0,0]; % covariance matrix state equation errors
R=1; % variance observation equation error
for k=1:n %Run the Kalman filter for each variable
if nargin < 4 || isempty(x_user) %no intial value for state, extrapolate back two periods from the observations
x=[2*y(1,k)-y(2,k);
3*y(1,k)-2*y(2,k)];
else
x=x_user(:,k);
end
if nargin < 4 || isempty(P_user) %no initial value for the MSE, set a rather high one
P= [1e5 0;
0 1e5];
else
P=P_user{k};
end
for j=1:T %Get the estimates for each period
[x,P]=kalman_update(F,H,Q,R,y(j,k),x,P); %get new state estimate and update recursion
ytrend(j,k)=x(2);%second state is trend estimate
end
end
if nargout==2
ycycle=y-ytrend;
end
if nargin==5 %user provided a discard parameter
ytrend=ytrend(discard+1:end,:);%Remove the first "discard" periods from the trend series
if nargout==2
ycycle=ycycle(discard+1:end,:);
end
end
end
function [x,P]=kalman_update(F,H,Q,R,obs,x,P)
% Updates the Kalman filter estimation of the state and MSE
S=H*P*H'+R;
K=F*P*H';
K=K/S;
x=F*x+K*(obs -H*x); %State estimate
Temp=F-K*H;
P=Temp*P*Temp';
P=P+Q+K*R*K';%MSE estimate
end

12
matlab/stoch_simul.m
View File

@ -107,6 +107,12 @@ if ~options_.noprint
headers = char('Variables',labels); headers = char('Variables',labels);
lh = size(labels,2)+2; lh = size(labels,2)+2;
dyntable(my_title,headers,labels,M_.Sigma_e,lh,10,6); 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 if options_.partial_information
skipline() skipline()
disp('SOLUTION UNDER PARTIAL INFORMATION') disp('SOLUTION UNDER PARTIAL INFORMATION')
@ -157,10 +163,10 @@ if options_.nomoments == 0
elseif options_.periods == 0 elseif options_.periods == 0
% There is no code for theoretical moments at 3rd order % There is no code for theoretical moments at 3rd order
if options_.order <= 2 if options_.order <= 2
disp_th_moments(oo_.dr,var_list); oo_=disp_th_moments(oo_.dr,var_list,M_,options_,oo_);
end end
else else
disp_moments(oo_.endo_simul,var_list); oo_=disp_moments(oo_.endo_simul,var_list,M_,options_,oo_);
end end
end end
@ -350,7 +356,7 @@ if options_.irf
end end
if options_.SpectralDensity.trigger == 1 if options_.SpectralDensity.trigger == 1
[omega,f] = UnivariateSpectralDensity(oo_.dr,var_list); [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list);
end end

83
matlab/th_autocovariances.m
View File

@ -1,8 +1,8 @@
function [Gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition) function [Gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition)
% Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process % Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process
% with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_ % with coefficients dr.ghx and dr.ghu and shock variances Sigma_e
% for a subset of variables ivar (indices in lgy_) % for a subset of variables ivar.
% Theoretical HP-filtering is available as an option % Theoretical HP-filtering and band-pass filtering is available as an option
% %
% INPUTS % INPUTS
% dr: [structure] Reduced form solution of the DSGE model (decisions rules) % 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 end
end end
if options_.hp_filter == 0 if options_.hp_filter == 0 && ~options_.bandpass.indicator
v = NaN*ones(nvar,nvar); v = NaN*ones(nvar,nvar);
v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb'; v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb';
k = find(abs(v) < 1e-12); k = find(abs(v) < 1e-12);
@ -207,37 +207,44 @@ if options_.hp_filter == 0
end end
end end
end end
else% ==> Theoretical HP filter. else% ==> Theoretical filters.
% By construction, all variables are stationary when HP filtered % By construction, all variables are stationary when HP filtered
iky = inv_order_var(ivar); iky = inv_order_var(ivar);
stationary_vars = (1:length(ivar))'; stationary_vars = (1:length(ivar))';
aa = ghx(iky,:); aa = ghx(iky,:); %R in Uhlig (2001)
bb = ghu(iky,:); bb = ghu(iky,:); %S in Uhlig (2001)
lambda = options_.hp_filter; lambda = options_.hp_filter;
ngrid = options_.hp_ngrid; ngrid = options_.hp_ngrid;
freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); %[0,2*pi)
tpos = exp( sqrt(-1)*freqs); tpos = exp( sqrt(-1)*freqs); %positive frequencies
tneg = exp(-sqrt(-1)*freqs); tneg = exp(-sqrt(-1)*freqs); %negative frequencies
hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2); if options_.bandpass.indicator
mathp_col = []; 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)); IA = eye(size(A,1));
IE = eye(M_.exo_nbr); IE = eye(M_.exo_nbr);
for ig = 1:ngrid for ig = 1:ngrid
if hp1(ig)==0, if filter_gain(ig)==0,
f_hp = zeros(length(ivar),length(ivar)); f_hp = zeros(length(ivar),length(ivar));
else else
f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]... f_omega =(1/(2*pi))*([(IA-A*tneg(ig))\ghu1;IE]...
*M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) ... *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % spectral density of state variables; top formula Uhlig (2001), p. 20 with N=0
IE]); % state variables 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'
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; top formula Uhlig (2001), p. 21;
f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
end end
mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row mathp_col(ig,:) = (f_hp(:))'; % store as matrix row for ifft
% for ifft
end; end;
% Covariance of filtered series % 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); Gamma_y{1} = reshape(imathp_col(1,:),nvar,nvar);
% Autocorrelations % Autocorrelations
if nar > 0 if nar > 0
@ -253,44 +260,40 @@ else% ==> Theoretical HP filter.
Gamma_y{nar+2} = ones(nvar,1); Gamma_y{nar+2} = ones(nvar,1);
else else
Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr); 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)'; cs = chol(SS)';
SS = cs*cs'; SS = cs*cs';
b1(:,exo_names_orig_ord) = ghu1; b1(:,exo_names_orig_ord) = ghu1;
b2(:,exo_names_orig_ord) = ghu(iky,:); b2(:,exo_names_orig_ord) = ghu(iky,:);
mathp_col = []; mathp_col = NaN(ngrid,length(ivar)^2);
IA = eye(size(A,1)); IA = eye(size(A,1));
IE = eye(M_.exo_nbr); IE = eye(M_.exo_nbr);
for ig = 1:ngrid for ig = 1:ngrid
if hp1(ig)==0, if filter_gain(ig)==0,
f_hp = zeros(length(ivar),length(ivar)); f_hp = zeros(length(ivar),length(ivar));
else else
f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]... f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\b1;IE]...
*SS*[b1'*inv(IA-A'*tpos(ig)) ... *SS*[b1'/(IA-A'*tpos(ig)) IE]); % spectral density of state variables; top formula Uhlig (2001), p. 20 with N=0
IE]); % state variables 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'
g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series; top formula Uhlig (2001), p. 21;
f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
end end
mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row mathp_col(ig,:) = (f_hp(:))'; % store as matrix row for ifft
% for ifft
end; end;
imathp_col = real(ifft(mathp_col))*(2*pi); imathp_col = real(ifft(mathp_col))*(2*pi);
vv = diag(reshape(imathp_col(1,:),nvar,nvar)); vv = diag(reshape(imathp_col(1,:),nvar,nvar));
for i=1:M_.exo_nbr for i=1:M_.exo_nbr
mathp_col = []; mathp_col = NaN(ngrid,length(ivar)^2);
SSi = cs(:,i)*cs(:,i)'; SSi = cs(:,i)*cs(:,i)';
for ig = 1:ngrid for ig = 1:ngrid
if hp1(ig)==0, if filter_gain(ig)==0,
f_hp = zeros(length(ivar),length(ivar)); f_hp = zeros(length(ivar),length(ivar));
else else
f_omega =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]... f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\b1;IE]...
*SSi*[b1'*inv(IA-A'*tpos(ig)) ... *SSi*[b1'/(IA-A'*tpos(ig)) IE]); % spectral density of state variables; top formula Uhlig (2001), p. 20 with N=0
IE]); % state variables 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'
g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series; top formula Uhlig (2001), p. 21;
f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
end end
mathp_col = [mathp_col ; (f_hp(:))']; % store as matrix row mathp_col(ig,:) = (f_hp(:))'; % store as matrix row for ifft
% for ifft
end; end;
imathp_col = real(ifft(mathp_col))*(2*pi); imathp_col = real(ifft(mathp_col))*(2*pi);
Gamma_y{nar+2}(:,i) = abs(diag(reshape(imathp_col(1,:),nvar,nvar)))./vv; Gamma_y{nar+2}(:,i) = abs(diag(reshape(imathp_col(1,:),nvar,nvar)))./vv;

5
tests/Makefile.am
View File

@ -7,6 +7,11 @@ MODFILES = \
estimation/MH_recover/fs2000_recover.mod \ estimation/MH_recover/fs2000_recover.mod \
estimation/t_proposal/fs2000_student.mod \ estimation/t_proposal/fs2000_student.mod \
estimation/TaRB/fs2000_tarb.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/ls2003.mod \
gsa/ls2003a.mod \ gsa/ls2003a.mod \
gsa/cod_ML_morris/cod_ML_morris.mod \ gsa/cod_ML_morris/cod_ML_morris.mod \

3
tests/TeX/fs2000_corr_ME.mod
View File

@ -123,7 +123,8 @@ end;
steady; 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_original_model;
write_latex_static_model; write_latex_static_model;

116
tests/moments/example1_bp_test.mod Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
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

159
tests/moments/example1_hp_test.mod Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
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

71
tests/moments/example1_one_sided_hp_test.mod Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
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);

89
tests/moments/example1_var_decomp.mod Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
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

76
tests/moments/test_AR1_spectral_density.mod Normal file
View File

@ -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));