Fix UnivariateSpectralDensity.m for filtered variables
- Adds bandpass filter - Writes results to oo_time-shift
parent
25df325899
commit
a4b04ca9b4
|
@ -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,87 +93,72 @@ 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
|
||||||
tmp(abs(tmp) < 1e-12) = 0;
|
tneg = exp(-sqrt(-1)*freqs); %negative frequencies
|
||||||
Gamma(:,1) = diag(tmp);
|
if options_.hp_filter == 0 && ~options_.bandpass.indicator %do not filter
|
||||||
vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
|
filter_gain=ones(ngrid,1);
|
||||||
tmp = aa*vxy;
|
elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && options_.bandpass.indicator %filter with bandpass
|
||||||
tmp(abs(tmp) < 1e-12) = 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;
|
||||||
tmp(abs(tmp) < 1e-12) = 0;
|
filter_gain(freqs<=-2*pi/lowest_periodicity+2*pi & freqs>=-2*pi/highest_periodicity+2*pi)=1;
|
||||||
Gamma(:,i+1) = diag(tmp);
|
elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && ~options_.bandpass.indicator %filter with HP-filter
|
||||||
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 = NaN(ngrid,length(ivar)^2);
|
|
||||||
IA = eye(size(A,1));
|
|
||||||
IE = eye(M_.exo_nbr);
|
|
||||||
for ig = 1:ngrid
|
|
||||||
f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\ghu1;IE]...
|
|
||||||
*M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % state variables
|
|
||||||
g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
|
|
||||||
f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
|
|
||||||
mathp_col(ig,:) = (f_hp(:))'; % store as matrix row
|
|
||||||
% for ifft
|
|
||||||
end;
|
|
||||||
imathp_col = real(ifft(mathp_col))*(2*pi);
|
|
||||||
tmp = reshape(imathp_col(1,:),nvar,nvar);
|
|
||||||
tmp(abs(tmp)<1e-12) = 0;
|
|
||||||
Gamma(:,1) = diag(tmp);
|
|
||||||
sy = sqrt(Gamma(:,1));
|
|
||||||
sy = sy *sy';
|
|
||||||
for i=1:cutoff-1
|
|
||||||
tmp = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
|
|
||||||
tmp(abs(tmp) < 1e-12) = 0;
|
|
||||||
Gamma(:,i+1) = diag(tmp);
|
|
||||||
end
|
|
||||||
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
|
||||||
|
|
|
@ -350,7 +350,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
|
||||||
|
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue