Merge branch 'master' into remove-dynDate-class

Conflicts:
	preprocessor/DynareBison.yy
time-shift
Stéphane Adjemian (Penelope) 2013-11-14 16:41:08 +01:00
commit 0c00151092
27 changed files with 663 additions and 60 deletions

View File

@ -827,12 +827,14 @@ Allows Dynare to issue a warning and continue processing when
@outputhead
Depending on the computing tasks requested in the @file{.mod} file,
executing command @code{dynare} will leave in the workspace variables
containing results available for further processing. More details are
given under the relevant computing tasks.
executing the @code{dynare} command will leave variables containing
results in the workspace available for further processing. More
details are given under the relevant computing tasks.
The @code{M_}, @code{oo_} and @code{options_} structures are also saved
in a file called @file{@var{FILENAME}_results.mat}.
The @code{M_}, @code{oo_}, and @code{options_} structures are saved in
a file called @file{@var{FILENAME}_results.mat}. If they exist,
@code{estim_params_}, @code{bayestopt_}, @code{dataset_}, and
@code{estimation_info} are saved in the same file.
@examplehead
@ -3358,6 +3360,10 @@ The exogenous variables for which to compute IRFs. Default: all.
Requests the computation of normalized IRFs in percentage of the
standard error of each shock.
@item irf_plot_threshold = @var{DOUBLE}
@anchor{irf_plot_threshold}
Threshold size for plotting IRFs. All IRFs for a particular variable with a maximum absolute deviation from the steady state smaller than this value are not displayed. Default: @code{1e-10}.
@item nocorr
Don't print the correlation matrix (printing them is the default).
@ -4201,9 +4207,10 @@ end;
@end deffn
@deffn Block estimated_params_init ;
@deffnx Block estimated_params_init (@var{OPTIONS}@dots{});
This block declares numerical initial values for the optimizer when
these ones are different from the prior mean.
these ones are different from the prior mean. It should be specified after the @code{estimated_params}-block as otherwise the specified starting values are overwritten by the latter.
Each line has the following syntax:
@ -4212,6 +4219,14 @@ stderr VARIABLE_NAME | corr VARIABLE_NAME_1, VARIABLE_NAME_2 | PARAMETER_NAME
, INITIAL_VALUE;
@end example
@optionshead
@table @code
@item use_calibration
For not specifically intialized parameters, use the deep parameters and the elements of the covariance matrix specified in the @code{shocks}-block from calibration as starting values for estimation. For components of the @code{shocks}-block that were not explicitly specified during calibration or which violate the prior, the prior mean is used.
@end table
@xref{estimated_params}, for the meaning and syntax of the various components.
@end deffn
@ -4896,6 +4911,9 @@ using a standard Kalman filter.
@xref{irf_shocks}. Only used if @ref{bayesian_irf} is passed. Cannot be used
with @ref{dsge_var}.
@item irf_plot_threshold = @var{DOUBLE}
@xref{irf_plot_threshold}. Only used if @ref{bayesian_irf} is passed.
@item aim_solver
@xref{aim_solver}.
@ -7058,11 +7076,14 @@ The total number of draws is equal to
Tuning period for Metropolis-Hasting draws. Default: @code{30,000}
@item save_draws
Save all elements of @math{A^0}, @math{A^+}, @math{\xi}, and the
transition matrix to a file named @code{draws_<<file_tag>>.out} with
each draw on a separate line. A file that describes how these matrices
are laid out is contained in
@code{draws_header_<<file_tag>>.out}. Default: @code{off}
Save all elements of @math{A^0}, @math{A^+}, @math{Q}, and
@math{\zeta}, to a file named @code{draws_<<file_tag>>.out} with each
draw on a separate line. A file that describes how these matrices are
laid out is contained in @code{draws_header_<<file_tag>>.out}. A file
called @code{load_flat_file.m} is provided to simplify loading the
saved files into the corresponding variables @code{A0}, @code{Aplus},
@code{Q}, and @code{Zeta} in your MATLAB/Octave workspace. Default:
@code{off}
@end table
@end deffn

View File

@ -404,7 +404,7 @@ if options_.TeX
TEXNAMES = [];
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
name = deblank(varlist(j,:));
texname = deblank(varlist_TeX(j,:));

View File

@ -1,6 +1,7 @@
function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% Generates and stores Posterior IRFs
% PARALLEL CONTEXT
% This function perfom in parallel a portion of PosteriorIRF.m code.
% This function perfoms in parallel execution a portion of the PosteriorIRF.m code.
% This is a special kind of parallel function. Unlike of other parallel functions,
% that running in parallel a 'for' cycle, this function run in parallel a
% 'while' loop! The parallelization of 'while' loop (when possible) is a more

View File

@ -1,6 +1,8 @@
function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Generates the Posterior IRFs plot from the IRFs generated in
% PosteriorIRF_core1
% PARALLEL CONTEXT
% Perfome in parallel a portion of PosteriorIRF.m code.
% Perform in parallel execution a portion of the PosteriorIRF.m code.
% See also the comment in random_walk_metropolis_hastings_core.m funtion.
%
% INPUTS
@ -17,7 +19,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% SPECIAL REQUIREMENTS.
% None.
%
% Copyright (C) 2006-2012 Dynare Team
% Copyright (C) 2006-2013 Dynare Team
%
% This file is part of Dynare.
%
@ -87,7 +89,7 @@ for i=fpar:npar,
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > 10^(-6)
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
subplotnum = subplotnum+1;
if subplotnum == 1 && options_.relative_irf
hh = dyn_figure(options_,'Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
@ -131,6 +133,10 @@ for i=fpar:npar,
end
name = deblank(varlist(j,:));
title(name,'Interpreter','none')
else
if options_.debug
fprintf('POSTERIOR_IRF: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n',deblank(varlist(j,:)),tit(i,:),options_.impulse_responses.plot_threshold)
end
end
if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0)

View File

@ -0,0 +1,51 @@
function check_prior_bounds(xparam1,bounds,M_,estim_params_,options_,bayestopt_)
% function check_prior_bounds(xparam1,bounds,M_,estim_params_,options_)
% checks the parameter vector of violations of the prior bounds
% Inputs:
% -xparam1 [double] vector of parameters to be estimated (initial values)
% -bounds [vector] vector containing the lower and upper
% bounds
% -M_ [structure] characterizing the model.
% -estim_params_ [structure] characterizing parameters to be estimated
% -options_ [structure] characterizing the options
% -bayestopt_ [structure] characterizing priors
% Copyright (C) 2013 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/>.
outside_bound_pars=find(xparam1 < bounds(:,1) | xparam1 > bounds(:,2));
if ~isempty(outside_bound_pars)
for ii=1:length(outside_bound_pars)
outside_bound_par_names{ii,1}=get_the_name(outside_bound_pars(ii),0,M_,estim_params_,options_);
end
disp_string=[outside_bound_par_names{1,:}];
for ii=2:size(outside_bound_par_names,1)
disp_string=[disp_string,', ',outside_bound_par_names{ii,:}];
end
error(['Initial value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0. If you used the mode_file-option, check whether your mode-file is consistent with the priors.'])
end
inadmissible_inverse_gamma_values=find(bayestopt_.pshape==4 & xparam1 == 0);
if ~isempty(inadmissible_inverse_gamma_values)
for ii=1:length(inadmissible_inverse_gamma_values)
inadmissible_inverse_gamma_par_names{ii,1}=get_the_name(inadmissible_inverse_gamma_values(ii),0,M_,estim_params_,options_);
end
disp_string=[inadmissible_inverse_gamma_par_names{1,:}];
for ii=2:size(inadmissible_inverse_gamma_par_names,1)
disp_string=[disp_string,', ',inadmissible_inverse_gamma_par_names{ii,:}];
end
error(['Initial value(s) of ', disp_string ,' is zero. This is not allowed when using an inverse gamma prior.\n'])
end

View File

@ -0,0 +1,143 @@
function [xparam1,estim_params_,xparam1_explicitly_initialized,xparam1_properly_calibrated]=do_parameter_initialization(estim_params_,xparam1_calib,xparam1_NaN_set_to_prior_mean)
% function [xparam1,estim_params_]=get_initialized_parameters(estim_params_,xparam1_calib)
% gets explicitly initialized variables and properly calibrated parameters
%
% INPUTS
% o estim_params_ [structure] characterizing parameters to be estimated.
% o xparam1_calib [double] vector of parameters to be estimated, with parameters
% initialized from calibration using get_all_parameters
%
% o xparam1_NaN_set_to_prior_mean [double] vector of parameters to be estimated, with parameters
% initialized using dynare_estimation_init; not explicitly initialized
% parameters are at prior mean
% OUTPUTS
% o xparam1 [double] vector of initialized parameters; uses the hierarchy: 1) explicitly initialized parameters,
% 2) calibrated parameters, 3) prior mean
% o estim_params_ [structure] characterizing parameters to be estimated; it is
% updated here to reflect calibrated parameters
% o xparam1_explicitly_initialized [double] vector of parameters to be estimated that
% were explicitly initialized
% o xparam1_properly_calibrated [double] vector of parameters to be estimated that
% were properly calibrated
%
% SPECIAL REQUIREMENTS
% None
% Copyright (C) 2013 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/>.
nvx = size(estim_params_.var_exo,1);
nvn = size(estim_params_.var_endo,1);
ncx = size(estim_params_.corrx,1);
ncn = size(estim_params_.corrn,1);
np = size(estim_params_.param_vals,1);
estim_params_.nvx = nvx; %exogenous shock variances
estim_params_.nvn = nvn; %endogenous variances, i.e. measurement error
estim_params_.ncx = ncx; %exogenous shock correlations
estim_params_.ncn = ncn; % correlation between endogenous variables, i.e. measurement error.
estim_params_.np = np; % other parameters of the model
xparam1_explicitly_initialized = NaN(nvx+nvn+ncx+ncn+np,1);
xparam1_properly_calibrated = NaN(nvx+nvn+ncx+ncn+np,1);
offset=0;
if nvx
initialized_par_index=find(~isnan(estim_params_.var_exo(:,2)));
calibrated_par_index=find(isnan(estim_params_.var_exo(:,2)) & ~isnan(xparam1_calib(offset+1:offset+nvx,1)));
uninitialized_par_index=find(isnan(estim_params_.var_exo(:,2)) & isnan(xparam1_calib(offset+1:offset+nvx,1)));
xparam1_explicitly_initialized(offset+initialized_par_index,1) = estim_params_.var_exo(initialized_par_index,2);
%update estim_params_ with calibrated starting values
estim_params_.var_exo(calibrated_par_index,2)=xparam1_calib(offset+calibrated_par_index,1);
%find parameters that are calibrated and do not violate inverse gamma prior
xparam1_properly_calibrated(offset+calibrated_par_index,1) = xparam1_calib(offset+calibrated_par_index,1);
inv_gamma_violation=find(estim_params_.var_exo(calibrated_par_index,2)==0 & estim_params_.var_exo(calibrated_par_index,4)==0);
if inv_gamma_violation
estim_params_.var_exo(calibrated_par_index(inv_gamma_violation),2)=NaN;
xparam1_properly_calibrated(offset+calibrated_par_index(inv_gamma_violation),1)=NaN;
fprintf('PARAMETER INITIALIZATION: Some standard deviations of shocks of the calibrated model are 0 and\n')
fprintf('PARAMETER INITIALIZATION: violate the inverse gamma prior. They will instead be initialized with the prior mean.\n')
end
if uninitialized_par_index
fprintf('PARAMETER INITIALIZATION: Warning, some estimated standard deviations of shocks are not\n')
fprintf('PARAMETER INITIALIZATION: initialized. They will be initialized with the prior mean.\n')
end
end
offset=offset+nvx;
if nvn
initialized_par_index=find(~isnan(estim_params_.var_endo(:,2)));
calibrated_par_index=find(isnan(estim_params_.var_endo(:,2)) & ~isnan(xparam1_calib(offset+1:offset+nvn,1)));
uninitialized_par_index=find(isnan(estim_params_.var_endo(:,2)) & isnan(xparam1_calib(offset+1:offset+nvn,1)));
xparam1_explicitly_initialized(offset+initialized_par_index,1) = estim_params_.var_endo(initialized_par_index,2);
estim_params_.var_endo(calibrated_par_index,2)=xparam1_calib(offset+calibrated_par_index,1);
%find parameters that are calibrated and do not violate inverse gamma prior
xparam1_properly_calibrated(offset+calibrated_par_index,1) = xparam1_calib(offset+calibrated_par_index,1);
inv_gamma_violation=find(estim_params_.var_endo(calibrated_par_index,2)==0 & estim_params_.var_endo(calibrated_par_index,5)==4);
if inv_gamma_violation
estim_params_.var_endo(calibrated_par_index(inv_gamma_violation),2)=NaN;
xparam1_properly_calibrated(offset+calibrated_par_index(inv_gamma_violation),1)=NaN;
fprintf('PARAMETER INITIALIZATION: Some measurement errors of the calibrated model are 0 and violate the\n')
fprintf('PARAMETER INITIALIZATION: inverse gamma prior. They will instead be initialized with the prior mean.\n')
end
if uninitialized_par_index
fprintf('PARAMETER INITIALIZATION: Warning, some measurement errors are not initialized. They will be initialized\n')
fprintf('PARAMETER INITIALIZATION: with the prior mean.\n')
end
end
offset=offset+nvn;
if ncx
initialized_par_index=find(~isnan(estim_params_.corrx(:,3)));
calibrated_par_index=find(isnan(estim_params_.corrx(:,3)) & ~isnan(xparam1_calib(offset+1:offset+ncx,1)));
uninitialized_par_index=find(isnan(estim_params_.corrx(:,3)) & isnan(xparam1_calib(offset+1:offset+ncx,1)));
xparam1_explicitly_initialized(offset+initialized_par_index,1) = estim_params_.corrx(initialized_par_index,3);
estim_params_.corrx(calibrated_par_index,3)=xparam1_calib(offset+calibrated_par_index,1);
xparam1_properly_calibrated(offset+calibrated_par_index,1) = xparam1_calib(offset+calibrated_par_index,1);
if uninitialized_par_index
fprintf('PARAMETER INITIALIZATION: Warning, some correlations between structural shocks are not initialized.\n')
fprintf('PARAMETER INITIALIZATION: They will be initialized with the prior mean.\n')
end
end
offset=offset+ncx;
if ncn
initialized_par_index=find(~isnan(estim_params_.corrn(:,3)));
calibrated_par_index=find(isnan(estim_params_.corrn(:,3)) & ~isnan(xparam1_calib(offset+1:offset+ncn,1)));
uninitialized_par_index=find(isnan(estim_params_.corrn(:,3)) & isnan(xparam1_calib(offset+1:offset+ncn,1)));
xparam1_explicitly_initialized(offset+initialized_par_index,1) = estim_params_.corrn(initialized_par_index,3);
estim_params_.corrn(calibrated_par_index,3)=xparam1_calib(offset+calibrated_par_index,1);
xparam1_properly_calibrated(offset+calibrated_par_index,1) = xparam1_calib(offset+calibrated_par_index,1);
if uninitialized_par_index
fprintf('PARAMETER INITIALIZATION: Warning, some correlations between measurement errors are not initialized.\n')
fprintf('PARAMETER INITIALIZATION: They will be initialized with the prior mean.\n')
end
end
offset=offset+ncn;
if np
initialized_par_index=find(~isnan(estim_params_.param_vals(:,2)));
calibrated_par_index=find(isnan(estim_params_.param_vals(:,2)) & ~isnan(xparam1_calib(offset+1:offset+np,1)));
uninitialized_par_index=find(isnan(estim_params_.param_vals(:,2)) & isnan(xparam1_calib(offset+1:offset+np,1)));
xparam1_explicitly_initialized(offset+initialized_par_index,1) = estim_params_.param_vals(initialized_par_index,2);
estim_params_.param_vals(calibrated_par_index,2)=xparam1_calib(offset+calibrated_par_index,1);
xparam1_properly_calibrated(offset+calibrated_par_index,1) = xparam1_calib(offset+calibrated_par_index,1);
if uninitialized_par_index
fprintf('PARAMETER INITIALIZATION: Warning, some deep parameters are not initialized. They will be\n')
fprintf('PARAMETER INITIALIZATION: initialized with the prior mean.\n')
end
end
xparam1=xparam1_explicitly_initialized;
xparam1(isnan(xparam1))=xparam1_properly_calibrated(isnan(xparam1)); %set not explicitly initialized parameters that do not obviously violate prior distribution to calibrated parameter values
xparam1(isnan(xparam1))=xparam1_NaN_set_to_prior_mean(isnan(xparam1)); %set not yet initialized parameters to prior mean coming from dynare_estimation_init

View File

@ -84,6 +84,25 @@ end
if ~isempty(estim_params_)
estim_params_=check_for_calibrated_covariances(xparam1,estim_params_,M_);
end
%%read out calibration that was set in mod-file and can be used for initialization
xparam1_calib=get_all_parameters(estim_params_,M_); %get calibrated parameters
if ~any(isnan(xparam1_calib)) %all estimated parameters are calibrated
full_calibration_detected=1;
else
full_calibration_detected=0;
end
if options_.use_calibration_initialization %set calibration as starting values
[xparam1,estim_params_]=do_parameter_initialization(estim_params_,xparam1_calib,xparam1); %get explicitly initialized parameters that have precedence to calibrated values
try
check_prior_bounds(xparam1,[bayestopt_.lb bayestopt_.ub],M_,estim_params_,options_,bayestopt_); %check whether calibration satisfies prior bounds
catch
e = lasterror();
fprintf('Cannot use parameter values from calibration as they violate the prior bounds.')
rethrow(e);
end
end
% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file).
M_.sigma_e_is_diagonal = 1;
if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
@ -152,7 +171,21 @@ if options_.dsge_var
end
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,M_,estim_params_,options_,bayestopt_,oo_);
%% perform initial estimation checks;
try
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,M_,estim_params_,options_,bayestopt_,oo_);
catch % if check fails, provide info on using calibration if present
e = lasterror();
if full_calibration_detected %calibrated model present and no explicit starting values
skipline(1);
fprintf('ESTIMATION_CHECKS: There was an error in computing the likelihood for initial parameter values.\n')
fprintf('ESTIMATION_CHECKS: You should try using the calibrated version of the model as starting values. To do\n')
fprintf('ESTIMATION_CHECKS: this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n')
fprintf('ESTIMATION_CHECKS: command (and after the estimated_params-block so that it does not get overwritten):\n');
skipline(2);
end
rethrow(e);
end
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
if options_.smoother == 1

View File

@ -242,28 +242,7 @@ if ~isempty(estim_params_)
bounds(:,2) = ub;
end
% Test if initial values of the estimated parameters are all between the prior lower and upper bounds.
outside_bound_pars=find(xparam1 < bounds(:,1) | xparam1 > bounds(:,2));
if ~isempty(outside_bound_pars)
for ii=1:length(outside_bound_pars)
outside_bound_par_names{ii,1}=get_the_name(outside_bound_pars(ii),0,M_,estim_params_,options_);
end
disp_string=[outside_bound_par_names{1,:}];
for ii=2:size(outside_bound_par_names,1)
disp_string=[disp_string,', ',outside_bound_par_names{ii,:}];
end
error(['Initial value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0. If you used the mode_file-option, check whether your mode-file is consistent with the priors.'])
end
inadmissible_inverse_gamma_values=find(bayestopt_.pshape==4 & xparam1 == 0);
if ~isempty(inadmissible_inverse_gamma_values)
for ii=1:length(inadmissible_inverse_gamma_values)
inadmissible_inverse_gamma_par_names{ii,1}=get_the_name(inadmissible_inverse_gamma_values(ii),0,M_,estim_params_,options_);
end
disp_string=[inadmissible_inverse_gamma_par_names{1,:}];
for ii=2:size(inadmissible_inverse_gamma_par_names,1)
disp_string=[disp_string,', ',inadmissible_inverse_gamma_par_names{ii,:}];
end
error(['Initial value(s) of ', disp_string ,' is zero. This is not allowed when using an inverse gamma prior.\n'])
end
check_prior_bounds(xparam1,bounds,M_,estim_params_,options_,bayestopt_)
lb = bounds(:,1);
ub = bounds(:,2);

View File

@ -0,0 +1,96 @@
function xparam1=get_all_parameters(estim_params_,M_)
% function xparam1=get_parameters
% gets parameters values from M_.params into xparam1 (inverse mapping to set_all_parameters)
% This is called if a model was calibrated before estimation to back out
% parameter values
%
% INPUTS
% estim_params_: Dynare structure describing the estimated parameters.
% M_: Dynare structure describing the model.
%
% OUTPUTS
% xparam1: N*1 double vector of parameters from calibrated model that are to be estimated
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013 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/>.
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
np = estim_params_.np;
Sigma_e = M_.Sigma_e;
Correlation_matrix = M_.Correlation_matrix;
H = M_.H;
Correlation_matrix_ME = M_.Correlation_matrix_ME;
xparam1=NaN(nvx+ncx+nvn+ncn+np,1);
% stderrs of the exogenous shocks
if nvx
var_exo = estim_params_.var_exo;
for i=1:nvx
k = var_exo(i,1);
xparam1(i)=sqrt(Sigma_e(k,k));
end
end
% update offset
offset = nvx;
% setting measument error variance
if nvn
for i=1:nvn
k = estim_params_.nvn_observable_correspondence(i,1);
xparam1(offset+i)=sqrt(H(k,k));
end
end
% update offset
offset = nvx+nvn;
% correlations among shocks (ncx)
if ncx
corrx = estim_params_.corrx;
for i=1:ncx
k1 = corrx(i,1);
k2 = corrx(i,2);
xparam1(i+offset)=Correlation_matrix(k1,k2);
end
end
% update offset
offset = nvx+nvn+ncx;
if ncn
corrn_observable_correspondence = estim_params_.corrn_observable_correspondence;
for i=1:ncn
k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2);
xparam1(i+offset)=Correlation_matrix_ME(k1,k2);
end
end
% update offset
offset = nvx+ncx+nvn+ncn;
% structural parameters
if np
xparam1(offset+1:end)=M_.params(estim_params_.param_vals(:,1));
end

View File

@ -140,6 +140,7 @@ end
% IRFs & other stoch_simul output
options_.irf = 40;
options_.impulse_responses.plot_threshold=1e-10;
options_.relative_irf = 0;
options_.ar = 5;
options_.hp_filter = 0;
@ -382,6 +383,7 @@ options_.mh_recover = 0;
options_.mh_replic = 20000;
options_.recursive_estimation_restart = 0;
options_.MCMC_jumping_covariance='hessian';
options_.use_calibration_initialization = 0;
options_.mode_compute = 4;
options_.mode_file = '';

View File

@ -263,7 +263,8 @@ if fload==0,
if prepSA,
try
T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
catch ME
catch
ME = lasterror();
if strcmp('MATLAB:nomem',ME.identifier),
prepSA=0;
disp('The model is too large for storing state space matrices ...')

View File

@ -146,8 +146,8 @@ for b=1:nb
end
end
if singularity_problem
fprint('The presence of a singularity problem typically indicates that there is one\n')
fprint('redundant equation entered in the model block, while another non-redundant equation\n')
fprint('is missing. The problem often derives from Walras Law.\n')
fprintf('The presence of a singularity problem typically indicates that there is one\n')
fprintf('redundant equation entered in the model block, while another non-redundant equation\n')
fprintf('is missing. The problem often derives from Walras Law.\n')
end

View File

@ -0,0 +1,71 @@
function [Q, A0, Aplus, Zeta] = load_flat_file(file_tag)
%function [Q, A0, Aplus, Zeta] = load_flat_file(file_tag)
% Loads file saved by save_draws option of ms_simulation.
%
% INPUTS
% file_tag: the file tag used with ms_simulation
%
% OUTPUTS
% Q: Q matrix
% A0: A0 matrix
% Aplus: A+ matrix
% Zeta: Zeta matrx
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013 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/>.
Q = [];
A0 = [];
Aplus = [];
Zeta = [];
try
headerfid = fopen(['est_flat_header_' file_tag '.out'], 'r');
catch
error(['Can''t find est_flat_header_' file_tag '.out'])
end
headerfile=textscan(headerfid, '%s');
fclose(headerfid);
flatfile = dlmread(['est_flat_' file_tag '.out'], ' ');
headerfile = headerfile{:};
for i=1:length(headerfile)
line = char(headerfile{i});
indob = strfind(line, '[');
indcb = strfind(line, ']');
indop = strfind(line, '(');
indcp = strfind(line, ')');
indc = strfind(line, ',');
if isempty(indob)
name = line(1:indop-1);
dim3 = [];
else
name = line(1:indob-1);
dim3 = line(indob+1:indcb-1);
end
row = line(indop+1:indc-1);
col = line(indc+1:indcp-1);
if isempty(dim3)
eval([name '(' row ',' col ') = ' num2str(flatfile(i)) ';']);
else
eval([name '(' row ',' col ',' dim3 ') = ' num2str(flatfile(i)) ';']);
end
end

View File

@ -67,8 +67,8 @@ end
%build covariance matrix from correlation matrix and variances already on
%diagonal
Sigma_e = diag(sqrt(diag(Sigma_e)))*Correlation_matrix*diag(sqrt(diag(Sigma_e)));
if isfield(estim_params,'calibrated_covariances')
Sigma_e(estim_params.calibrated_covariances.position)=estim_params.calibrated_covariances.cov_value;
if isfield(estim_params_,'calibrated_covariances')
Sigma_e(estim_params_.calibrated_covariances.position)=estim_params_.calibrated_covariances.cov_value;
end
% and update offset

View File

@ -34,7 +34,7 @@ function dr=set_state_space(dr,DynareModel,DynareOptions)
%! @end deftypefn
%@eod:
% Copyright (C) 1996-2012 Dynare Team
% Copyright (C) 1996-2013 Dynare Team
%
% This file is part of Dynare.
%
@ -84,6 +84,9 @@ if max_lag > 0
end
kmask = [kmask; lead_lag_incidence(1,order_var)] ;
else
if max_lead==0 %%in this case lead_lag_incidence has no entry max_lag+2
error('Dynare currently does not allow to solve purely static models in a stochastic context.')
end
kmask = lead_lag_incidence(max_lag+2,order_var) ;
end

View File

@ -195,7 +195,7 @@ if options_.irf
y(i_var(j),:)');
eval(['oo_.irfs.' deblank(M_.endo_names(i_var(j),:)) '_' ...
deblank(M_.exo_names(i,:)) ' = y(i_var(j),:);']);
if max(y(i_var(j),:)) - min(y(i_var(j),:)) > 1e-10
if max(y(i_var(j),:)) - min(y(i_var(j),:)) > options_.impulse_responses.plot_threshold
irfs = cat(1,irfs,y(i_var(j),:));
if isempty(mylist)
mylist = deblank(var_list(j,:));
@ -209,6 +209,10 @@ if options_.irf
mylistTeX = char(mylistTeX,deblank(var_listTeX(j,:)));
end
end
else
if options_.debug
fprintf('STOCH_SIMUL: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n',deblank(M_.endo_names(i_var(j),:)),deblank(M_.exo_names(i,:)),options_.impulse_responses.plot_threshold)
end
end
end
if options_.nograph == 0

View File

@ -581,15 +581,20 @@ EstimatedParamsStatement::writeOutput(ostream &output, const string &basename) c
}
EstimatedParamsInitStatement::EstimatedParamsInitStatement(const vector<EstimationParams> &estim_params_list_arg,
const SymbolTable &symbol_table_arg) :
const SymbolTable &symbol_table_arg,
const bool use_calibration_arg) :
estim_params_list(estim_params_list_arg),
symbol_table(symbol_table_arg)
symbol_table(symbol_table_arg),
use_calibration(use_calibration_arg)
{
}
void
EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basename) const
{
if (use_calibration)
output << "options_.use_calibration_initialization = 1;" << endl;
vector<EstimationParams>::const_iterator it;
for (it = estim_params_list.begin(); it != estim_params_list.end(); it++)

View File

@ -290,9 +290,11 @@ class EstimatedParamsInitStatement : public Statement
private:
const vector<EstimationParams> estim_params_list;
const SymbolTable &symbol_table;
const bool use_calibration;
public:
EstimatedParamsInitStatement(const vector<EstimationParams> &estim_params_list_arg,
const SymbolTable &symbol_table_arg);
const SymbolTable &symbol_table_arg,
const bool use_calibration_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};

View File

@ -90,7 +90,7 @@ class ParsingDriver;
%}
%token AIM_SOLVER ANALYTIC_DERIVATION AR AUTOCORR
%token BAYESIAN_IRF BETA_PDF BLOCK
%token BAYESIAN_IRF BETA_PDF BLOCK USE_CALIBRATION
%token BVAR_DENSITY BVAR_FORECAST
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
@ -106,7 +106,7 @@ class ParsingDriver;
%token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HOMOTOPY_FORCE_CONTINUE HP_FILTER HP_NGRID HYBRID
%token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE BOUNDS JSCALE INIT
%token <string_val> INT_NUMBER
%token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS
%token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS IRF_PLOT_THRESHOLD
%token KALMAN_ALGO KALMAN_TOL SUBSAMPLES OPTIONS TOLF
%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR 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
@ -990,6 +990,7 @@ stoch_simul_primary_options : o_dr_algo
| o_dr_cycle_reduction_tol
| o_dr_logarithmic_reduction_tol
| o_dr_logarithmic_reduction_maxiter
| o_irf_plot_threshold
;
stoch_simul_options : stoch_simul_primary_options
@ -1136,7 +1137,10 @@ estimated_elem3 : expression_or_empty COMMA expression_or_empty
;
estimated_params_init : ESTIMATED_PARAMS_INIT ';' estimated_init_list END ';'
{ driver.estimated_params_init(); };
{ driver.estimated_params_init(); }
| ESTIMATED_PARAMS_INIT '(' USE_CALIBRATION ')' ';' estimated_init_list END ';'
{ driver.estimated_params_init(true); }
;
estimated_init_list : estimated_init_list estimated_init_elem
{ driver.add_estimated_params_element(); }
@ -1562,6 +1566,7 @@ estimation_options : o_datafile
| o_taper_steps
| o_geweke_interval
| o_mcmc_jumping_covariance
| o_irf_plot_threshold
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@ -2689,6 +2694,7 @@ o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN
| MCMC_JUMPING_COVARIANCE EQUAL filename
{ driver.option_str("MCMC_jumping_covariance", $3); }
;
o_irf_plot_threshold : IRF_PLOT_THRESHOLD EQUAL non_negative_number { driver.option_num("impulse_responses.plot_threshold", $3); };
range : symbol ':' symbol
{

View File

@ -397,7 +397,7 @@ string eofbuff;
return token::CNUM;
}
<DYNARE_STATEMENT>banact {return token::BANACT;}
<DYNARE_BLOCK>use_calibration {return token::USE_CALIBRATION;}
<DYNARE_STATEMENT>output_file_tag {return token::OUTPUT_FILE_TAG;}
<DYNARE_STATEMENT>file_tag {return token::FILE_TAG;};
<DYNARE_STATEMENT>no_create_init {return token::NO_CREATE_INIT;};
@ -569,6 +569,7 @@ string eofbuff;
<DYNARE_STATEMENT>mh_recover {return token::MH_RECOVER;}
<DYNARE_STATEMENT>planner_discount {return token::PLANNER_DISCOUNT;}
<DYNARE_STATEMENT>calibration {return token::CALIBRATION;}
<DYNARE_STATEMENT>irf_plot_threshold {return token::IRF_PLOT_THRESHOLD;}
<DYNARE_BLOCK>equation {return token::EQUATION;}
<DYNARE_BLOCK>exclusion {return token::EXCLUSION;}

View File

@ -703,7 +703,15 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool no_log, b
if (block && !byte_code)
mOutputFile << "rmpath " << basename << ";" << endl;
mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl;
mOutputFile << "save('" << basename << "_results.mat', 'oo_', 'M_', 'options_');" << endl
<< "if exist('estim_params_', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'estim_params_', '-append');" << endl << "end" << endl
<< "if exist('bayestopt_', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'bayestopt_', '-append');" << endl << "end" << endl
<< "if exist('dataset_', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'dataset_', '-append');" << endl << "end" << endl
<< "if exist('estimation_info', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'estimation_info', '-append');" << endl << "end" << endl;
config_file.writeEndParallel(mOutputFile);

View File

@ -1209,9 +1209,9 @@ ParsingDriver::estimated_params()
}
void
ParsingDriver::estimated_params_init()
ParsingDriver::estimated_params_init(bool use_calibration)
{
mod_file->addStatement(new EstimatedParamsInitStatement(estim_params_list, mod_file->symbol_table));
mod_file->addStatement(new EstimatedParamsInitStatement(estim_params_list, mod_file->symbol_table, use_calibration));
estim_params_list.clear();
}

View File

@ -394,7 +394,7 @@ public:
//! Writes estimated params command
void estimated_params();
//! Writes estimated params init command
void estimated_params_init();
void estimated_params_init(bool use_calibration = false);
//! Writes estimated params bound command
void estimated_params_bounds();
//! Adds a declaration for a user-defined external function

View File

@ -4,6 +4,8 @@ MODFILES = \
estimation/fs2000_mc6.mod \
estimation/fs2000_mc6_mf.mod \
estimation/fs2000_MCMC_jumping_covariance.mod \
estimation/fs2000_initialize_from_calib.mod \
estimation/fs2000_calibrated_covariance.mod \
gsa/ls2003.mod \
ramst.mod \
ramst_a.mod \
@ -311,7 +313,8 @@ EXTRA_DIST = \
third_order/policyfunctions.mat \
shock_decomposition/example1_calib_shock_decomp_data.mat \
shock_decomposition/fsdat_simul.m \
estimation/fs2000_MCMC_jumping_covariance_steadystate.m
estimation/fs2000_MCMC_jumping_covariance_steadystate.m \
estimation/fs2000_initialize_from_calib_steadystate.m
TARGETS =

View File

@ -0,0 +1,94 @@
//uses a prior mean for bet that does not allow solving the model. Then uses the calibrated value to initialize.
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m e_f;
parameters alp bet gam mst rho psi del;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
model;
dA = exp(gam+e_a);
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
W = l/n;
-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
P*c = m;
m-1+d = l;
e = exp(e_a);
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
initval;
k = 6;
m = mst;
P = 2.25;
c = 0.45;
e = 1;
W = 4;
R = 1.02;
d = 0.85;
n = 0.19;
l = 0.86;
y = 0.6;
gy_obs = exp(gam);
gp_obs = exp(-gam);
dA = exp(gam);
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
var e_f; stderr 0.005;
corr gy_obs, gp_obs = 0.1;
var gy_obs; stderr 0.005;
end;
steady;
//check;
estimated_params;
alp, beta_pdf, 0.356, 0.02;
bet, normal_pdf, 2, 1;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
rho, normal_pdf, 0.129, 0.223;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
stderr e_m, inv_gamma_pdf, 0.008862, inf;
corr e_a, e_m, normal_pdf, 0, 0.007;
corr e_a, e_f, normal_pdf, 0, 0.007;
stderr gp_obs, inv_gamma_pdf, 0.008862, inf;
corr gp_obs, gy_obs, normal_pdf, 0, 0.007;
end;
varobs gp_obs gy_obs;
estimated_params_init(use_calibration);
//stderr e_a, 0.014000;
//stderr e_m, 0.005000;
//stderr gp_obs, 0.1;
corr e_m, e_a, 0.1;
//corr gy_obs, gp_obs, 0.000000;
//alp, 0.330000;
//bet, 0.990000;
gam, 0.003000;
mst, 1.011000;
//rho, 0.900000;
psi, 0.787000;
del, 0.020000;
end;
options_.plot_priors=0;
estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=2000, mh_nblocks=1, mh_jscale=0.8,prior_trunc=0);

View File

@ -0,0 +1,73 @@
% computes the steady state of fs2000 analyticaly
% largely inspired by the program of F. Schorfheide
% Copyright (C) 2004-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/>.
function [ys,check] = fs2000_steadystate(ys,exe)
global M_
alp = M_.params(1);
bet = M_.params(2);
gam = M_.params(3);
mst = M_.params(4);
rho = M_.params(5);
psi = M_.params(6);
del = M_.params(7);
check = 0;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
ist = y-c;
q = 1 - d;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
ys =[
m
P
c
e
W
R
k
d
n
l
gy_obs
gp_obs
y
dA ];

View File

@ -104,7 +104,7 @@ corr gp_obs, gy_obs,0;
end;
options_.TeX=1;
estimation(mode_compute=9,order=1,datafile=fsdat_simul,mode_check,smoother,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20) m P c e W R k d y gy_obs;
estimation(order=1,datafile=fsdat_simul,mode_check,smoother,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20) m P c e W R k d y gy_obs;