From 05189497a580e96f3522ce0b14dfcbc9ce6cf9eb Mon Sep 17 00:00:00 2001 From: michel Date: Sat, 4 Jul 2009 14:11:51 +0000 Subject: [PATCH] changed handling of nonstationary variables: - oo_.mean, oo_.var, oo_.gamma_y contains all selected variables - moments of non-stationary variables are set to NaN options_.Schur_vec_tolerance lowered to 10^-11 git-svn-id: https://www.dynare.org/svn/dynare/trunk@2810 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/disp_th_moments.m | 24 ++++++++++++-------- matlab/global_initialization.m | 6 ++--- matlab/th_autocovariances.m | 41 +++++++++++++++++++++------------- 3 files changed, 43 insertions(+), 28 deletions(-) diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index 52d423bde..ae123f832 100644 --- a/matlab/disp_th_moments.m +++ b/matlab/disp_th_moments.m @@ -35,8 +35,10 @@ function disp_th_moments(dr,var_list) end end - [oo_.gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_); + [oo_.gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_); m = dr.ys(ivar); + non_stationary_vars = setdiff(1:length(ivar),stationary_vars); + m(ivar(non_stationary_vars)) = NaN; i1 = find(abs(diag(oo_.gamma_y{1})) > 1e-12); @@ -50,14 +52,15 @@ function disp_th_moments(dr,var_list) oo_.mean = m; oo_.var = oo_.gamma_y{1}; - lh = size(deblank(M_.endo_names(ivar,:)),2)+2; if options_.nomoments == 0 title='THEORETICAL MOMENTS'; if options_.hp_filter title = [title ' (HP filter, lambda = ' int2str(options_.hp_filter) ')']; end headers=strvcat('VARIABLE','MEAN','STD. DEV.','VARIANCE'); - dyntable(title,headers,deblank(M_.endo_names(ivar,:)),z,lh,11,4); + labels = deblank(M_.endo_names(ivar,:)); + lh = size(labels,2)+2; + dyntable(title,headers,labels,z,lh,11,4); if M_.exo_nbr > 1 disp(' ') title='VARIANCE DECOMPOSITION (in percent)'; @@ -68,8 +71,9 @@ function disp_th_moments(dr,var_list) headers = M_.exo_names; headers(M_.exo_names_orig_ord,:) = headers; headers = strvcat(' ',headers); - dyntable(title,headers,deblank(M_.endo_names(ivar(i1),:)),100*oo_.gamma_y{options_.ar+2}(i1,:), ... - lh,8,2); + lh = size(deblank(M_.endo_names(ivar(stationary_vars),:)),2)+2; + dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ... + :)),100*oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2); end end @@ -79,10 +83,11 @@ function disp_th_moments(dr,var_list) if options_.hp_filter title = [title ' (HP filter, lambda = ' int2str(options_.hp_filter) ')']; end - labels = deblank(M_.endo_names(ivar,:)); - headers = strvcat('Variables',labels(i1,:)); + labels = deblank(M_.endo_names(ivar(i1),:)); + headers = strvcat('Variables',labels); corr = oo_.gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)'); - dyntable(title,headers,labels(i1,:),corr,lh,8,4); + lh = size(labels,2)+2; + dyntable(title,headers,labels,corr,lh,8,4); end if options_.ar > 0 @@ -98,7 +103,8 @@ function disp_th_moments(dr,var_list) oo_.autocorr{i} = oo_.gamma_y{i+1}; z(:,i) = diag(oo_.gamma_y{i+1}(i1,i1)); end - dyntable(title,headers,labels,z,0,8,4); + lh = size(labels,2)+2; + dyntable(title,headers,labels,z,lh,8,4); end % 10/09/02 MJ diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 5a916c18d..6f06b19f7 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -41,9 +41,9 @@ function global_initialization() options_.gstep = 1e-2; options_.debug = 0; options_.initval_file = 0; - options_.Schur_vec_tol = 1e-8; % used to find nonstationary variables - % in Schur decomposition of the - % transition matrix + options_.Schur_vec_tol = 1e-11; % used to find nonstationary variables + % in Schur decomposition of the + % transition matrix options_.qz_criterium = 1.000001; options_.lyapunov_complex_threshold = 1e-15; options_.solve_tolf = eps^(1/3); diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m index b60a692cb..12f4fa2cd 100644 --- a/matlab/th_autocovariances.m +++ b/matlab/th_autocovariances.m @@ -1,4 +1,4 @@ -function [Gamma_y,ivar] = 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 % with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_ % for a subset of variables ivar (indices in lgy_) @@ -19,7 +19,8 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition % Gamma_y{nar+2} [double] Variance decomposition. % Gamma_y{nar+3} [double] Expectation of the endogenous variables associated with a second % order approximation. -% ivar [integer] Vector of indices for a subset of variables. +% stationary_vars [integer] Vector of indices of stationary +% variables in declaration order % % SPECIAL REQUIREMENTS % @@ -45,6 +46,7 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition nodecomposition = 0; end + endo_nbr = M_.endo_nbr; exo_names_orig_ord = M_.exo_names_orig_ord; if exist('OCTAVE_VERSION') warning('off', 'Octave:divide-by-zero') @@ -54,7 +56,7 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition nar = options_.ar; Gamma_y = cell(nar+1,1); if isempty(ivar) - ivar = [1:M_.endo_nbr]'; + ivar = [1:endo_nbr]'; end nvar = size(ivar,1); @@ -63,8 +65,8 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition npred = dr.npred; nstatic = dr.nstatic; kstate = dr.kstate; - order = dr.order_var; - iv(order) = [1:length(order)]; + order_var = dr.order_var; + inv_order_var = dr.inv_order_var; nx = size(ghx,2); ikx = [nstatic+1:nstatic+npred]; @@ -94,10 +96,12 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition [A,B] = kalman_transition_matrix(dr,ipred,1:nx,dr.transition_auxiliary_variables,M_.exo_nbr); if options_.order == 2 | options_.hp_filter == 0 [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold); - iky = iv(ivar); + iky = inv_order_var(ivar); + stationary_vars = (1:length(ivar))'; if ~isempty(u) - iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2))); - ivar = dr.order_var(iky); + x = abs(ghx*u); + iky = iky(find(all(x(iky,:) < options_.Schur_vec_tol,2))); + stationary_vars = find(all(x(inv_order_var(ivar(stationary_vars)),:) < options_.Schur_vec_tol,2)); end aa = ghx(iky,:); bb = ghu(iky,:); @@ -109,23 +113,28 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition end end if options_.hp_filter == 0 - Gamma_y{1} = aa*vx*aa'+ bb*M_.Sigma_e*bb'; - k = find(abs(Gamma_y{1}) < 1e-12); - Gamma_y{1}(k) = 0; + v = NaN*ones(nvar,nvar); + v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb'; + k = find(abs(v) < 1e-12); + v(k) = 0; + Gamma_y{1} = v; % autocorrelations if nar > 0 vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb'); sy = sqrt(diag(Gamma_y{1})); + sy = sy(stationary_vars); sy = sy *sy'; - Gamma_y{2} = aa*vxy./sy; + Gamma_y{2} = NaN*ones(nvar,nvar); + Gamma_y{2}(stationary_vars,stationary_vars) = aa*vxy./sy; for i=2:nar vxy = A*vxy; - Gamma_y{i+1} = aa*vxy./sy; + Gamma_y{i+1} = NaN*ones(nvar,nvar); + Gamma_y{i+1}(stationary_vars,stationary_vars) = aa*vxy./sy; end end % variance decomposition if ~nodecomposition && M_.exo_nbr > 1 - Gamma_y{nar+2} = zeros(length(ivar),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); cs = chol(SS)'; b1(:,exo_names_orig_ord) = ghu1; @@ -136,12 +145,12 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition vv = diag(aa*vx*aa'+b2*b2'); for i=1:M_.exo_nbr vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,options_.lyapunov_complex_threshold,2); - Gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv; + Gamma_y{nar+2}(stationary_vars,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv; end end else% ==> Theoretical HP filter. if options_.order < 2 - iky = iv(ivar); + iky = inv_order_var(ivar); aa = ghx(iky,:); bb = ghu(iky,:); end