From d4c07260266484b38eeb751ab883a92200c89ef8 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 11 Mar 2015 18:15:28 +0100 Subject: [PATCH 01/13] Fix size of filtered variables matrix Fixes crash occurring due to non-conformable matrices --- matlab/prior_posterior_statistics.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index 35a0b481f..4908d6b64 100644 --- a/matlab/prior_posterior_statistics.m +++ b/matlab/prior_posterior_statistics.m @@ -98,13 +98,13 @@ MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8)); MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8)); if naK - MAX_naK = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)* ... - length(options_.filter_step_ahead)*gend)/8)); + MAX_naK = min(B,ceil(MaxNumberOfBytes/(endo_nbr* ... + length(options_.filter_step_ahead)*(gend+max(options_.filter_step_ahead)))/8)); end if horizon - MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*horizon)/8)); - MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*horizon)/ ... + MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8)); + MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ... 8)); IdObs = bayestopt_.mfys; From 214dc747237b8a3ddce4267039ada25d0ab2f0e6 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 1 Apr 2015 09:00:51 +0200 Subject: [PATCH 02/13] - Fixed bugs around analytic derivation. - Fixed test routine, eliminating diffuse filter. - Trapped incompatibility of diffuse filter with analytic derivation. --- matlab/dsge_likelihood.m | 2 +- matlab/dynare_estimation_init.m | 3 +++ matlab/optimization/dynare_minimize_objective.m | 3 --- tests/analytic_derivatives/fs2000_analytic_derivation.mod | 8 +++++--- 4 files changed, 9 insertions(+), 7 deletions(-) diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index 7ac811416..aae1224e4 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -631,7 +631,7 @@ if analytic_derivation, analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess}; else analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P}; - clear DT DYss DOm DH DP D2T D2Yss D2Om D2H D2P, + clear DT DYss DOm DP D2T D2Yss D2Om D2H D2P, end else analytic_deriv_info={0}; diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index e54618a7b..eab7b8b7e 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -442,6 +442,9 @@ else end; if options_.analytic_derivation, + if options_.lik_init == 3, + error('analytic derivation is incompatible with diffuse filter') + end options_.analytic_derivation = 1; if ~(exist('sylvester3','file')==2), dynareroot = strrep(which('dynare'),'dynare.m',''); diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m index b345b8d55..4113bf119 100644 --- a/matlab/optimization/dynare_minimize_objective.m +++ b/matlab/optimization/dynare_minimize_objective.m @@ -178,9 +178,6 @@ switch minimizer_algorithm end [opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,analytic_grad,crit,nit,0,varargin{:}); %hessian_mat is the plain outer product gradient Hessian - if options_.analytic_derivation %Hessian is already analytic one, reset option - options_.analytic_derivation = ana_deriv; - end case 6 [opt_par_values, hessian_mat, Scale, fval] = gmhmaxlik(objective_function, start_par_value, ... Initial_Hessian, options_.mh_jscale, bounds, prior_information.p2, options_.gmhmaxlik, options_.optim_opt, varargin{:}); diff --git a/tests/analytic_derivatives/fs2000_analytic_derivation.mod b/tests/analytic_derivatives/fs2000_analytic_derivation.mod index c4f2a7479..3dc8b7ad0 100644 --- a/tests/analytic_derivatives/fs2000_analytic_derivation.mod +++ b/tests/analytic_derivatives/fs2000_analytic_derivation.mod @@ -61,7 +61,7 @@ alp, beta_pdf, 0.356, 0.02; bet, beta_pdf, 0.993, 0.002; gam, normal_pdf, 0.0085, 0.003; mst, normal_pdf, 1.0002, 0.007; -rho, beta_pdf, 0.129, 0.223; +rho, beta_pdf, 0.129, 0.1; psi, beta_pdf, 0.65, 0.05; del, beta_pdf, 0.01, 0.005; stderr e_a, inv_gamma_pdf, 0.035449, inf; @@ -74,5 +74,7 @@ options_.solve_tolf = 1e-12; estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=1,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=2,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); -estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=3,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); -estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=4,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); \ No newline at end of file +estimation(order=1,mode_compute=4,analytic_derivation,kalman_algo=1,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); +estimation(order=1,mode_compute=4,analytic_derivation,kalman_algo=2,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); +//estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=3,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); +//estimation(order=1,mode_compute=5,analytic_derivation,kalman_algo=4,datafile=fsdat_simul,nobs=192,mh_replic=0,mh_nblocks=2,mh_jscale=0.8); \ No newline at end of file From 9c3c0d727f7fe336c02e580c57b13f5ab9f29314 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 17 Sep 2014 20:06:57 +0200 Subject: [PATCH 03/13] Add check for point priors and force user to fix this parameter. Traps Inf-prior --- matlab/set_prior.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/matlab/set_prior.m b/matlab/set_prior.m index 7bc3cc83d..a72a6204d 100644 --- a/matlab/set_prior.m +++ b/matlab/set_prior.m @@ -155,6 +155,12 @@ end bayestopt_.p6 = NaN(size(bayestopt_.p1)) ; bayestopt_.p7 = bayestopt_.p6 ; +%% check for point priors and disallow them as they do not work with MCMC +if any(bayestopt_.p2 ==0) + error(sprintf(['Error in prior for %s: you cannot use a point prior in estimation. Either increase the prior standard deviation',... + ' or fix the parameter completely.'], bayestopt_.name{bayestopt_.p2 ==0})) +end + % generalized location parameters by default for beta distribution k = find(bayestopt_.pshape == 1); k1 = find(isnan(bayestopt_.p3(k))); From 27403fb2dd495e578325c8e57d58fe448e31ac25 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 17 Sep 2014 20:09:33 +0200 Subject: [PATCH 04/13] Add check whether initial prior is Inf Requires adding fourth output argument of priordens.m that stores position of problematic parameter. --- matlab/initial_estimation_checks.m | 8 +++++++- matlab/priordens.m | 20 +++++++++++++++++++- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m index 54af16f7b..7b1bf2380 100644 --- a/matlab/initial_estimation_checks.m +++ b/matlab/initial_estimation_checks.m @@ -89,8 +89,14 @@ if any(BayesInfo.pshape) % if Bayesian estimation end end -%% display warning if some parameters are still NaN +% display warning if some parameters are still NaN test_for_deep_parameters_calibration(Model); + +[lnprior, junk1,junk2,info]= priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4); +if info + fprintf('The prior density evaluated at the initial values is Inf for the following parameters: %s\n',BayesInfo.name{info,1}) + error('The initial value of the prior is -Inf') +end % Evaluate the likelihood. ana_deriv = DynareOptions.analytic_derivation; diff --git a/matlab/priordens.m b/matlab/priordens.m index 02ba5c698..f3e5fd619 100644 --- a/matlab/priordens.m +++ b/matlab/priordens.m @@ -1,4 +1,4 @@ -function [logged_prior_density, dlprior, d2lprior] = priordens(x, pshape, p6, p7, p3, p4,initialization) +function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4,initialization) % Computes a prior density for the structural parameters of DSGE models % % INPUTS @@ -12,6 +12,7 @@ function [logged_prior_density, dlprior, d2lprior] = priordens(x, pshape, p6, p7 % % OUTPUTS % logged_prior_density [double] scalar, log of the prior density evaluated at x. +% info [double] error code for index of Inf-prior parameter % % Copyright (C) 2003-2012 Dynare Team @@ -34,6 +35,8 @@ function [logged_prior_density, dlprior, d2lprior] = priordens(x, pshape, p6, p7 persistent id1 id2 id3 id4 id5 id6 persistent tt1 tt2 tt3 tt4 tt5 tt6 +info=0; + if nargin > 6 && initialization == 1 % Beta indices. tt1 = 1; @@ -81,6 +84,9 @@ d2lprior = 0.0; if tt1 logged_prior_density = logged_prior_density + sum(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1))) ; if isinf(logged_prior_density) + if nargout ==4 + info=id1(isinf(lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1)))); + end return end if nargout == 2, @@ -93,6 +99,9 @@ end if tt2 logged_prior_density = logged_prior_density + sum(lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2))) ; if isinf(logged_prior_density) + if nargout ==4 + info=id2(isinf(lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2)))); + end return end if nargout == 2, @@ -114,6 +123,9 @@ end if tt4 logged_prior_density = logged_prior_density + sum(lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4))) ; if isinf(logged_prior_density) + if nargout ==4 + info=id4(isinf(lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4)))); + end return end if nargout == 2, @@ -126,6 +138,9 @@ end if tt5 if any(x(id5)-p3(id5)<0) || any(x(id5)-p4(id5)>0) logged_prior_density = -Inf ; + if nargout ==4 + info=id5((x(id5)-p3(id5)<0) || (x(id5)-p4(id5)>0)); + end return end logged_prior_density = logged_prior_density + sum(log(1./(p4(id5)-p3(id5)))) ; @@ -140,6 +155,9 @@ end if tt6 logged_prior_density = logged_prior_density + sum(lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6))) ; if isinf(logged_prior_density) + if nargout ==4 + info=id6(isinf(lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6)))); + end return end if nargout == 2, From 8cc9fbd69a78290da3edd3b802d840ca2f8ae7ad Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 2 Apr 2015 20:52:49 +0200 Subject: [PATCH 05/13] Remove hard-coded coloring of legend in mode_check.m Now displays correct colors with new Matlab plot interface (>2014a) --- matlab/mode_check.m | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/matlab/mode_check.m b/matlab/mode_check.m index fd80782d7..180ba3592 100644 --- a/matlab/mode_check.m +++ b/matlab/mode_check.m @@ -149,7 +149,7 @@ for plt = 1:nbplt, y(i,2) = (y(i,1)+lnprior-dy); end end - plot(z,-y); + fighandle=plot(z,-y); hold on yl=get(gca,'ylim'); plot( [x(kk) x(kk)], yl, 'c', 'LineWidth', 1) @@ -173,8 +173,9 @@ for plt = 1:nbplt, else axes('position',[0.3 0.01 0.42 0.05],'box','on'), end - plot([0.48 0.68],[0.5 0.5],'color',[0 0.5 0]) - hold on, plot([0.04 0.24],[0.5 0.5],'b') + line_color=get(fighandle,'color'); + plot([0.48 0.68],[0.5 0.5],'color',line_color{2}) + hold on, plot([0.04 0.24],[0.5 0.5],'color',line_color{1}) set(gca,'xlim',[0 1],'ylim',[0 1],'xtick',[],'ytick',[]) text(0.25,0.5,'log-post') text(0.69,0.5,'log-lik kernel') From 7c838b133f11ce58101bf611b770de74feec48c5 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 1 Apr 2015 14:55:47 +0200 Subject: [PATCH 06/13] Use local variables for restrict var list, to avoid overwriting the global and messing likelihood evaluations later on. --- matlab/DsgeSmoother.m | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index f59236c69..ef38cf231 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -69,9 +69,10 @@ M_ = set_all_parameters(xparam1,estim_params_,M_); %------------------------------------------------------------------------------ % 2. call model setup & reduction program %------------------------------------------------------------------------------ -oo_.dr.restrict_var_list = bayestopt_.smoother_var_list; -oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns; -[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_); +DynareResults = oo_; +DynareResults.dr.restrict_var_list = bayestopt_.smoother_var_list; +DynareResults.dr.restrict_columns = bayestopt_.smoother_restrict_columns; +[T,R,SteadyState,info,M_,options_,DynareResults] = dynare_resolve(M_,options_,DynareResults); bayestopt_.mf = bayestopt_.smoother_mf; if options_.noconstant constant = zeros(vobs,1); From 52978365772b67a7e4ff78e59b00b89438dbefd4 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 1 Apr 2015 14:58:10 +0200 Subject: [PATCH 07/13] Harmonize criteria for exiting diffuse steps in likelihood with the smoother. Since initial Pinf is well scaled to unity, crit1= 1.e-6 is used for smoother and should also apply to likelihood evaluations. --- matlab/kalman/likelihood/kalman_filter_d.m | 3 ++- .../likelihood/missing_observations_kalman_filter_d.m | 3 ++- matlab/kalman/likelihood/univariate_kalman_filter_d.m | 7 ++++--- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/matlab/kalman/likelihood/kalman_filter_d.m b/matlab/kalman/likelihood/kalman_filter_d.m index ced4d1305..bbec1fe7c 100644 --- a/matlab/kalman/likelihood/kalman_filter_d.m +++ b/matlab/kalman/likelihood/kalman_filter_d.m @@ -59,8 +59,9 @@ dlik = zeros(smpl,1); % Initialization of the vector gathering the densitie dLIK = Inf; % Default value of the log likelihood. oldK = Inf; s = 0; +crit1=1.e-6; -while rank(Pinf,kalman_tol) && (t<=last) +while rank(Pinf,crit1) && (t<=last) s = t-start+1; v = Y(:,t)-Z*a; Finf = Z*Pinf*Z'; diff --git a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m index 608a07b7d..c09194820 100644 --- a/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m +++ b/matlab/kalman/likelihood/missing_observations_kalman_filter_d.m @@ -66,13 +66,14 @@ t = start; % Initialization of the time index. dlik = zeros(smpl,1); % Initialization of the vector gathering the densities. dLIK = Inf; % Default value of the log likelihood. oldK = Inf; +crit1=1.e-6; if isequal(H,0) H = zeros(pp,pp); end s = 0; -while rank(Pinf,kalman_tol) && (t<=last) +while rank(Pinf,crit1) && (t<=last) s = t-start+1; d_index = data_index{t}; if isempty(d_index) diff --git a/matlab/kalman/likelihood/univariate_kalman_filter_d.m b/matlab/kalman/likelihood/univariate_kalman_filter_d.m index ba047a291..e3368a6c1 100644 --- a/matlab/kalman/likelihood/univariate_kalman_filter_d.m +++ b/matlab/kalman/likelihood/univariate_kalman_filter_d.m @@ -110,7 +110,8 @@ dLIK = Inf; % Default value of the log likelihood. oldK = Inf; llik = zeros(smpl,pp); -newRank = rank(Pinf,kalman_tol); +crit1 = 1.e-6; +newRank = rank(Pinf,crit1); l2pi = log(2*pi); while newRank && (t<=last) @@ -138,7 +139,7 @@ while newRank && (t<=last) end end if newRank - oldRank = rank(Pinf,kalman_tol); + oldRank = rank(Pinf,crit1); else oldRank = 0; end @@ -146,7 +147,7 @@ while newRank && (t<=last) Pstar = T*Pstar*T'+QQ; Pinf = T*Pinf*T'; if newRank - newRank = rank(Pinf,kalman_tol); + newRank = rank(Pinf,crit1); end if oldRank ~= newRank disp('univariate_diffuse_kalman_filter:: T does influence the rank of Pinf!') From e97e5c340757098a2904b24aaeb5f8011812cdf3 Mon Sep 17 00:00:00 2001 From: Marco Ratto Date: Wed, 1 Apr 2015 15:10:26 +0200 Subject: [PATCH 08/13] Exclude zero columns of T from Kitagawa transformation, since ordschur is extremely noisy for multiple zero eigenvalues. This can make a lot of difference for large models that have hundreds of definitions. --- matlab/schur_statespace_transformation.m | 61 +++++++++++++++++++++--- 1 file changed, 55 insertions(+), 6 deletions(-) diff --git a/matlab/schur_statespace_transformation.m b/matlab/schur_statespace_transformation.m index cef99ed78..c9bdec8a9 100644 --- a/matlab/schur_statespace_transformation.m +++ b/matlab/schur_statespace_transformation.m @@ -1,8 +1,15 @@ -function [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_criterium) -% function [Z,ST,QT,R1,Pstar,Pinf] = schur_statespace(mf,T,R,Q,qz_criterium) +function [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_criterium, restrict_columns) +% function [Z,ST,QT,R1,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_criterium, restrict_columns) % Kitagawa transformation of state space system with a quasi-triangular -% transition matrix with unit roots at the top. Computation of Pstar and -% Pinf for Durbin and Koopman Diffuse filter +% transition matrix with unit roots at the top, but excluding zero columns of the transition matrix. +% Computation of Pstar and Pinf for Durbin and Koopman Diffuse filter +% +% The transformed state space is +% y = [ss; z; x]; +% s = static variables (zero columns of T) +% z = unit roots +% x = stable roots +% ss = s - z = stationarized static variables % % INPUTS % mf [integer] vector of indices of observed variables in @@ -44,6 +51,20 @@ function [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_c % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +np = size(T,1); +if nargin == 6, + indx = restrict_columns; + indx0=find(~ismember([1:np],indx)); +else + indx=(find(max(abs(T))>0)); + indx0=(find(max(abs(T))==0)); +end +T0=T(indx0,indx); %static variables vs. dynamic ones +R0=R(indx0,:); % matrix of shocks for static variables + +% perform Kitagawa transformation only for non-zero columns of T +T=T(indx,indx); +R=R(indx,:); np = size(T,1); [QT,ST] = schur(T); e1 = abs(ordeig(ST)) > 2-qz_criterium; @@ -93,14 +114,40 @@ if i == nk+1 Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1)); end -Z = QT(mf,:); R1 = QT'*R; +ST1=ST; + +% now I recover stationarized static variables +% using +% ss = s-z and +% z-z(-1) (growth rates of unit roots) only depends on stationary variables +np0=length(indx0); +Pstar = blkdiag(zeros(np0),Pstar); +ST = [zeros(length(Pstar),length(indx0)) [T0*QT ;ST]]; +R1 = [R0; R1]; +ST0=ST; +ST0(:,1:np0+nk)=0; +ST0(np0+1:np0+nk,:)=0; +ST0(1:np0,np0+nk+1:end) = ST(1:np0,np0+nk+1:end)-ST(1:np0,np0+1:np0+nk)*ST(np0+1:np0+nk,np0+nk+1:end); +R10 = R1; +R10(np0+1:np0+nk,:)=0; +R10(1:np0,:) = R1(1:np0,:)-ST(1:np0,np0+1:np0+nk)*R1(np0+1:np0+nk,:); +Pstar = ... + ST0*Pstar*ST0'+ ... + R10*Q*R10'; +QT = blkdiag(eye(np0),QT); +QT(1:np0,np0+1:np0+nk) = QT(1:np0,np0+1:np0+nk)+ST(1:np0,np0+1:np0+nk); +% reorder QT entries +QT([indx0(:); indx(:)],:) = QT; +Z = QT(mf,:); +ST(1:np0,:) = ST0(1:np0,:); +R1(1:np0,:) = R10(1:np0,:); % stochastic trends with no influence on observed variables are % arbitrarily initialized to zero Pinf = zeros(np,np); Pinf(1:nk,1:nk) = eye(nk); -[QQ,RR,EE] = qr(Z*ST(:,1:nk),0); +[QQ,RR,EE] = qr(Z*ST(:,1+np0:nk+np0),0); k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8); if length(k) > 0 k1 = EE(:,k); @@ -108,3 +155,5 @@ if length(k) > 0 dd(k1) = zeros(length(k1),1); Pinf(1:nk,1:nk) = diag(dd); end +Pinf = blkdiag(zeros(np0),Pinf); + From 090c4fedbd28dcd184fdcdaccc750d922ec4bf59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 3 Apr 2015 17:48:25 +0200 Subject: [PATCH 09/13] Added new option (diffuse_kalman_tol) and fixed tolerance paremeters in diffuse smoother routines. --- doc/dynare.texi | 2 ++ matlab/DsgeSmoother.m | 5 +++-- matlab/global_initialization.m | 1 + matlab/missing_DiffuseKalmanSmootherH1_Z.m | 16 ++++++++-------- matlab/missing_DiffuseKalmanSmootherH3_Z.m | 14 +++++++------- preprocessor/DynareBison.yy | 4 +++- preprocessor/DynareFlex.ll | 1 + 7 files changed, 25 insertions(+), 18 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 0e60bdf3e..094089cb4 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5221,6 +5221,8 @@ singularity is encountered, Dynare by default automatically switches to the univ @item kalman_tol = @var{DOUBLE} @anchor{kalman_tol} Numerical tolerance for determining the singularity of the covariance matrix of the prediction errors during the Kalman filter (minimum allowed reciprocal of the matrix condition number). Default value is @code{1e-10} +@item diffuse_kalman_tol = @var{DOUBLE} +@anchor{diffuse_kalman_tol} Numerical tolerance for determining the singularity of the covariance matrix of the prediction errors (@math{F_{\infty}}) and the rank of the covariance matrix of the state variables (@math{P_{\infty}}) during the Diffuse Kalman filter. Default value is @code{1e-8} @item filter_covariance @anchor{filter_covariance} Saves the series of one step ahead error of diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index ef38cf231..2dba85016 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -175,6 +175,7 @@ elseif options_.lik_init == 5 % Old diffuse Kalman filter only for th Pinf = []; end kalman_tol = options_.kalman_tol; +diffuse_kalman_tol = options_.diffuse_kalman_tol; riccati_tol = options_.riccati_tol; data1 = Y-trend; % ----------------------------------------------------------------------------- @@ -200,7 +201,7 @@ if kalman_algo == 1 || kalman_algo == 3 [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(ST, ... Z,R1,Q,H,Pinf,Pstar, ... data1,vobs,np,smpl,data_index, ... - options_.nk,kalman_tol,options_.filter_decomposition); + options_.nk,kalman_tol,diffuse_kalman_tol,options_.filter_decomposition); if isinf(alphahat) if kalman_algo == 1 kalman_algo = 2; @@ -227,7 +228,7 @@ if kalman_algo == 2 || kalman_algo == 4 [alphahat,epsilonhat,etahat,ahat,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH3_Z(ST, ... Z,R1,Q,diag(H), ... Pinf,Pstar,data1,vobs,np,smpl,data_index, ... - options_.nk,kalman_tol,... + options_.nk,kalman_tol,diffuse_kalman_tol, ... options_.filter_decomposition); end diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 1470d533f..f169bfe6a 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -387,6 +387,7 @@ options_.nobs = NaN; options_.kalman_algo = 0; options_.fast_kalman = 0; options_.kalman_tol = 1e-10; +options_.diffuse_kalman_tol = 1e-8; options_.use_univariate_filters_if_singularity_is_detected = 1; options_.riccati_tol = 1e-6; options_.lik_algo = 1; diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m index 4b5014771..9fc09eb30 100644 --- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m +++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m @@ -1,6 +1,6 @@ -function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag) +function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKalmanSmootherH1_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag) -% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag) +% function [alphahat,epsilonhat,etahat,a,aK,PK,decomp] = DiffuseKalmanSmoother1(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag) % Computes the diffuse kalman smoother without measurement error, in the case of a non-singular var-cov matrix % % INPUTS @@ -18,6 +18,7 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % data_index [cell] 1*smpl cell of column vectors of indices. % nk number of forecasting periods % kalman_tol tolerance for reciprocal condition number +% diffuse_kalman_tol tolerance for reciprocal condition number (for Finf) and the rank of Pinf % decomp_flag if true, compute filter decomposition % % OUTPUTS @@ -38,7 +39,7 @@ function [alphahat,epsilonhat,etahat,atilde,P,aK,PK,decomp] = missing_DiffuseKal % Models", S.J. Koopman and J. Durbin (2003, in Journal of Time Series % Analysis, vol. 24(1), pp. 85-98). -% Copyright (C) 2004-2011 Dynare Team +% Copyright (C) 2004-2015 Dynare Team % % This file is part of Dynare. % @@ -80,7 +81,6 @@ Kstar = zeros(mm,pp,smpl); P = zeros(mm,mm,smpl+1); Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1; Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1; -crit1 = 1.e-8; steady = smpl; rr = size(Q,1); QQ = R*Q*transpose(R); @@ -91,7 +91,7 @@ epsilonhat = zeros(rr,smpl); r = zeros(mm,smpl+1); t = 0; -while rank(Pinf(:,:,t+1),crit1) && t INT_NUMBER %token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS IRF_PLOT_THRESHOLD IRF_CALIBRATION -%token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS TOLF +%token KALMAN_ALGO KALMAN_TOL DIFFUSE_KALMAN_TOL SUBSAMPLES OPTIONS TOLF %token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LOGDATA LYAPUNOV %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT %token MFS MH_CONF_SIG MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS @@ -1664,6 +1664,7 @@ estimation_options : o_datafile | o_filtered_vars | o_kalman_algo | o_kalman_tol + | o_diffuse_kalman_tol | o_xls_sheet | o_xls_range | o_filter_step_ahead @@ -2643,6 +2644,7 @@ o_filtered_vars : FILTERED_VARS { driver.option_num("filtered_vars", "1"); }; o_relative_irf : RELATIVE_IRF { driver.option_num("relative_irf", "1"); }; o_kalman_algo : KALMAN_ALGO EQUAL INT_NUMBER { driver.option_num("kalman_algo", $3); }; o_kalman_tol : KALMAN_TOL EQUAL non_negative_number { driver.option_num("kalman_tol", $3); }; +o_diffuse_kalman_tol : DIFFUSE_KALMAN_TOL EQUAL non_negative_number { driver.option_num("diffuse_kalman_tol", $3); }; o_marginal_density : MARGINAL_DENSITY EQUAL LAPLACE { driver.option_str("mc_marginal_density", "laplace"); } | MARGINAL_DENSITY EQUAL MODIFIEDHARMONICMEAN diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 994aaa421..0e7a31c16 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -288,6 +288,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 nodiagnostic {return token::NODIAGNOSTIC;} kalman_algo {return token::KALMAN_ALGO;} kalman_tol {return token::KALMAN_TOL;} +diffuse_kalman_tol {return token::DIFFUSE_KALMAN_TOL;} forecast {return token::FORECAST;} smoother {return token::SMOOTHER;} bayesian_irf {return token::BAYESIAN_IRF;} From 6aa5ceabb2de37ded2315fa9b143203170fbf403 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 3 Apr 2015 17:59:51 +0200 Subject: [PATCH 10/13] Fixed bug introduced in 363c4a61222103fda22adff9bdb1115a6073428a. --- matlab/DsgeSmoother.m | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 2dba85016..a692965f1 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -69,10 +69,16 @@ M_ = set_all_parameters(xparam1,estim_params_,M_); %------------------------------------------------------------------------------ % 2. call model setup & reduction program %------------------------------------------------------------------------------ -DynareResults = oo_; -DynareResults.dr.restrict_var_list = bayestopt_.smoother_var_list; -DynareResults.dr.restrict_columns = bayestopt_.smoother_restrict_columns; -[T,R,SteadyState,info,M_,options_,DynareResults] = dynare_resolve(M_,options_,DynareResults); +oldoo.restrict_var_list = oo_.dr.restrict_var_list; +oldoo.restrict_columns = oo_.dr.restrict_columns; +oo_.dr.restrict_var_list = bayestopt_.smoother_var_list; +oo_.dr.restrict_columns = bayestopt_.smoother_restrict_columns; + +[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_); + +oo_.dr.restrict_var_list = oldoo.restrict_var_list; +oo_.dr.restrict_columns = oldoo.restrict_columns; + bayestopt_.mf = bayestopt_.smoother_mf; if options_.noconstant constant = zeros(vobs,1); From 104e4767a473d7298c76741720668c85150cb1ce Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 7 Apr 2015 12:01:23 +0200 Subject: [PATCH 11/13] update license file, #874 --- license.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/license.txt b/license.txt index a354f2a03..52fce2880 100644 --- a/license.txt +++ b/license.txt @@ -39,16 +39,16 @@ License: public-domain Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3, pages 472-489 -Files: matlab/bfgsi1.m matlab/csolve.m matlab/csminit1.m matlab/numgrad2.m - matlab/numgrad3.m matlab/numgrad3_.m matlab/numgrad5.m - matlab/numgrad5_.m matlab/csminwel1.m matlab/bvar_density.m +Files: matlab/bfgsi1.m matlab/csolve.m matlab/optimization/csminit1.m matlab/optimization/numgrad2.m + matlab/optimization/numgrad3.m matlab/optimization/numgrad3_.m matlab/optimization/numgrad5.m + matlab/optimization/numgrad5_.m matlab/optimization/csminwel1.m matlab/bvar_density.m matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m matlab/qzswitch.m matlab/qzdiv.m Copyright: 1993-2009 Christopher Sims 2006-2013 Dynare Team License: GPL-3+ -Files: matlab/cmaes.m +Files: matlab/optimization/cmaes.m Copyright: 2001-2012 Nikolaus Hansen 2012 Dynare Team License: GPL-3+ @@ -63,7 +63,7 @@ Copyright: 2008-2012 VZLU Prague, a.s. 2014 Dynare Team License: GPL-3+ -Files: matlab/simpsa.m matlab/simpsaget.m matlab/simpsaset.m +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 2013 Dynare Team From d8b132e39a701a2d8dcda305e0c02497bf6c7f75 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Apr 2015 12:36:07 +0200 Subject: [PATCH 12/13] preprocessor: remove commented code --- preprocessor/ComputingTasks.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 31545e5e3..a342d041b 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -1940,8 +1940,6 @@ void EstimationDataStatement::writeOutput(ostream &output, const string &basename) const { options_list.writeOutput(output, "options_.dataset"); - //if (options_list.date_options.find("first_obs") == options_list.date_options.end()) - // output << "options_.dataset.first_obs = options_.initial_period;" << endl; } SubsamplesStatement::SubsamplesStatement(const string &name1_arg, From 44bb0672d27241752137bf298a40c7c6d3261fca Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 14 Apr 2015 11:33:02 +0200 Subject: [PATCH 13/13] preprocessor: method out of order --- preprocessor/ComputingTasks.cc | 42 +++++++++++++++++----------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index a342d041b..889694726 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -2501,27 +2501,6 @@ CorrPriorStatement::writeOutput(ostream &output, const string &basename) const writePriorOutput(output, lhs_field, name1); } -PriorEqualStatement::PriorEqualStatement(const string &to_declaration_type_arg, - const string &to_name1_arg, - const string &to_name2_arg, - const string &to_subsample_name_arg, - const string &from_declaration_type_arg, - const string &from_name1_arg, - const string &from_name2_arg, - const string &from_subsample_name_arg, - const SymbolTable &symbol_table_arg) : - to_declaration_type(to_declaration_type_arg), - to_name1(to_name1_arg), - to_name2(to_name2_arg), - to_subsample_name(to_subsample_name_arg), - from_declaration_type(from_declaration_type_arg), - from_name1(from_name1_arg), - from_name2(from_name2_arg), - from_subsample_name(from_subsample_name_arg), - symbol_table(symbol_table_arg) -{ -} - void CorrPriorStatement::writeCOutput(ostream &output, const string &basename) { @@ -2554,6 +2533,27 @@ CorrPriorStatement::writeCOutput(ostream &output, const string &basename) output << endl <<" index, index1, shape, mean, mode, stdev, variance, domain));" << endl; } +PriorEqualStatement::PriorEqualStatement(const string &to_declaration_type_arg, + const string &to_name1_arg, + const string &to_name2_arg, + const string &to_subsample_name_arg, + const string &from_declaration_type_arg, + const string &from_name1_arg, + const string &from_name2_arg, + const string &from_subsample_name_arg, + const SymbolTable &symbol_table_arg) : + to_declaration_type(to_declaration_type_arg), + to_name1(to_name1_arg), + to_name2(to_name2_arg), + to_subsample_name(to_subsample_name_arg), + from_declaration_type(from_declaration_type_arg), + from_name1(from_name1_arg), + from_name2(from_name2_arg), + from_subsample_name(from_subsample_name_arg), + symbol_table(symbol_table_arg) +{ +} + void PriorEqualStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) {