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

Merge branch 'master' into johannes_perfect_foresight_test

time-shift
Stéphane Adjemian(Charybdis) 2016-03-24 23:01:11 +01:00
parent d84bf3e107 e714dc1a9f
commit b688cbe5b2
94 changed files with 3953 additions and 470 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
doc/dynare.texi
View File

@ -3977,6 +3977,12 @@ Default value: @code{1e-6}.
Saves the contemporaneous correlation between the endogenous variables in @code{oo_.contemporaneous_correlation}. Saves the contemporaneous correlation between the endogenous variables in @code{oo_.contemporaneous_correlation}.
Requires the @code{nocorr}-option not to be set. Requires the @code{nocorr}-option not to be set.
@item spectral_density
@anchor{spectral_density}
Triggers the computation and display of the theoretical spectral density of the (filtered) model variables.
Results are stored in @code{oo_.SpectralDensity}, defined below.
Default: do not request spectral density estimates
@end table @end table
@outputhead @outputhead
@ -4112,6 +4118,13 @@ After a run of @code{stoch_simul} with the
and simulated contemporaneous correlations otherwise. The variables are arranged in declaration order. and simulated contemporaneous correlations otherwise. The variables are arranged in declaration order.
@end defvr @end defvr
@anchor{oo_.SpectralDensity}
@defvr {MATLAB/Octave variable} oo_.SpectralDensity
After a run of @code{stoch_simul} with option @code{spectral_density} contains the spectral density
of the model variables. There will be a @code{nvars} by @code{nfrequencies} subfield
@code{freqs} storing the respective frequency grid points ranging from 0 to 2*pi and a
same sized subfield @code{density} storing the corresponding density.
@end defvr
@defvr {MATLAB/Octave variable} oo_.irfs @defvr {MATLAB/Octave variable} oo_.irfs
After a run of @code{stoch_simul} with option @code{irf} different After a run of @code{stoch_simul} with option @code{irf} different
@ -4857,7 +4870,7 @@ first observation of the rolling window.
@item prefilter = @var{INTEGER} @item prefilter = @var{INTEGER}
@anchor{prefilter} @anchor{prefilter}
A value of @code{1} means that the estimation procedure will demean A value of @code{1} means that the estimation procedure will demean
each data series by its empirical mean. Default: @code{0}, @i{i.e.} no prefiltering each data series by its empirical mean. If the (@ref{loglinear}) option without the (@ref{logdata}) option is requested, the data will first be logged and then demeaned. Default: @code{0}, @i{i.e.} no prefiltering
@item presample = @var{INTEGER} @item presample = @var{INTEGER}
@anchor{presample} @anchor{presample}
@ -4918,12 +4931,18 @@ variables
@item 2 @item 2
For nonstationary models: a wide prior is used with an initial matrix For nonstationary models: a wide prior is used with an initial matrix
of variance of the error of forecast diagonal with 10 on the diagonal of variance of the error of forecast diagonal with 10 on the diagonal
(follows the suggestion of @cite{Harvey and Phillips(1979)})
@item 3 @item 3
For nonstationary models: use a diffuse filter (use rather the @code{diffuse_filter} option) For nonstationary models: use a diffuse filter (use rather the @code{diffuse_filter} option)
@item 4 @item 4
The filter is initialized with the fixed point of the Riccati equation The filter is initialized with the fixed point of the Riccati equation
@item 5
Use i) option 2 for the non-stationary elements by setting their initial variance in the
forecast error matrix to 10 on the diagonal and all covariances to 0 and ii) option 1 for the stationary elements.
@end table @end table
@noindent @noindent
@ -4986,8 +5005,12 @@ Metropolis-Hastings chain. Default: 2*@code{mh_scale}
@item mh_recover @item mh_recover
@anchor{mh_recover} Attempts to recover a Metropolis-Hastings @anchor{mh_recover} Attempts to recover a Metropolis-Hastings
simulation that crashed prematurely. Shouldn't be used together with simulation that crashed prematurely, starting with the last available saved
@code{load_mh_file} @code{mh}-file. Shouldn't be used together with
@code{load_mh_file} or a different @code{mh_replic} than in the crashed run.
To assure a neat continuation of the chain with the same proposal density, you should
provide the @code{mode_file} used in the previous
run or the same user-defined @code{mcmc_jumping_covariance} when using this option.
@item mh_mode = @var{INTEGER} @item mh_mode = @var{INTEGER}
@dots{} @dots{}
@ -5624,8 +5647,7 @@ using a standard Kalman filter.
@xref{irf}. Only used if @ref{bayesian_irf} is passed. @xref{irf}. Only used if @ref{bayesian_irf} is passed.
@item irf_shocks = ( @var{VARIABLE_NAME} [[,] @var{VARIABLE_NAME} @dots{}] ) @item irf_shocks = ( @var{VARIABLE_NAME} [[,] @var{VARIABLE_NAME} @dots{}] )
@xref{irf_shocks}. Only used if @ref{bayesian_irf} is passed. Cannot be used @xref{irf_shocks}. Only used if @ref{bayesian_irf} is passed.
with @ref{dsge_var}.
@item irf_plot_threshold = @var{DOUBLE} @item irf_plot_threshold = @var{DOUBLE}
@xref{irf_plot_threshold}. Only used if @ref{bayesian_irf} is passed. @xref{irf_plot_threshold}. Only used if @ref{bayesian_irf} is passed.
@ -13130,6 +13152,11 @@ Hansen, Nikolaus and Stefan Kern (2004): ``Evaluating the CMA Evolution Strategy
on Multimodal Test Functions''. In: @i{Eighth International Conference on Parallel on Multimodal Test Functions''. In: @i{Eighth International Conference on Parallel
Problem Solving from Nature PPSN VIII, Proceedings}, Berlin: Springer, 282--291 Problem Solving from Nature PPSN VIII, Proceedings}, Berlin: Springer, 282--291
@item
Harvey, Andrew C. and Garry D.A. Phillips (1979): ``Maximum likelihood estimation of
regression models with autoregressive-moving average disturbances,''
@i{Biometrika}, 66(1), 49--58
@item @item
Herbst, Edward (2015): Herbst, Edward (2015):
``Using the ``Chandrasekhar Recursions'' for Likelihood Evaluation of DSGE ``Using the ``Chandrasekhar Recursions'' for Likelihood Evaluation of DSGE

18
matlab/DsgeSmoother.m
View File

@ -1,4 +1,4 @@
function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value) function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,decomp,trend_addition] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value)
% Estimation of the smoothed variables and innovations. % Estimation of the smoothed variables and innovations.
% %
% INPUTS % INPUTS
@ -9,10 +9,10 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% o missing_value 1 if missing values, 0 otherwise % o missing_value 1 if missing values, 0 otherwise
% %
% OUTPUTS % OUTPUTS
% o alphahat [double] (m*T) matrix, smoothed endogenous variables. % o alphahat [double] (m*T) matrix, smoothed endogenous variables (a_{t|T})
% o etahat [double] (r*T) matrix, smoothed structural shocks (r>n is the umber of shocks). % o etahat [double] (r*T) matrix, smoothed structural shocks (r>n is the umber of shocks).
% o epsilonhat [double] (n*T) matrix, smoothed measurement errors. % o epsilonhat [double] (n*T) matrix, smoothed measurement errors.
% o ahat [double] (m*T) matrix, one step ahead filtered (endogenous) variables. % o ahat [double] (m*T) matrix, updated (endogenous) variables (a_{t|t})
% o SteadyState [double] (m*1) vector specifying the steady state level of each endogenous variable. % o SteadyState [double] (m*1) vector specifying the steady state level of each endogenous variable.
% o trend_coeff [double] (n*1) vector, parameters specifying the slope of the trend associated to each observed variable. % o trend_coeff [double] (n*1) vector, parameters specifying the slope of the trend associated to each observed variable.
% o aK [double] (K,n,T+K) array, k (k=1,...,K) steps ahead filtered (endogenous) variables. % o aK [double] (K,n,T+K) array, k (k=1,...,K) steps ahead filtered (endogenous) variables.
@ -23,6 +23,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% matrices (meaningless for periods 1:d) % matrices (meaningless for periods 1:d)
% o decomp 4D array of shock decomposition of k-step ahead % o decomp 4D array of shock decomposition of k-step ahead
% filtered variables % filtered variables
% o trend_addition [double] (n_observed_series*T) pure trend component; stored in options_.varobs order
% %
% ALGORITHM % ALGORITHM
% Diffuse Kalman filter (Durbin and Koopman) % Diffuse Kalman filter (Durbin and Koopman)
@ -93,15 +94,10 @@ else
end end
trend_coeff = zeros(vobs,1); trend_coeff = zeros(vobs,1);
if bayestopt_.with_trend == 1 if bayestopt_.with_trend == 1
trend_coeff = zeros(vobs,1); [trend_addition, trend_coeff] =compute_trend_coefficients(M_,options_,vobs,gend);
t = options_.trend_coeffs; trend = constant*ones(1,gend)+trend_addition;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = evalin('base',t{i});
end
end
trend = constant*ones(1,gend)+trend_coeff*(1:gend);
else else
trend_addition=zeros(size(constant,1),gend);
trend = constant*ones(1,gend); trend = constant*ones(1,gend);
end end
start = options_.presample+1; start = options_.presample+1;

3
matlab/TaRB_metropolis_hastings_core.m
View File

@ -236,7 +236,8 @@ for curr_chain = fblck:nblck,
dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_chain) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/blocked_draws_counter_this_chain)]); dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_chain) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/blocked_draws_counter_this_chain)]);
end end
if (draw_index_current_file == InitSizeArray(curr_chain)) || (draw_iter == nruns(curr_chain)) % Now I save the simulations, either because the current file is full or the chain is done if (draw_index_current_file == InitSizeArray(curr_chain)) || (draw_iter == nruns(curr_chain)) % Now I save the simulations, either because the current file is full or the chain is done
save([BaseName '_mh' int2str(NewFile(curr_chain)) '_blck' int2str(curr_chain) '.mat'],'x2','logpo2'); [LastSeeds.(['file' int2str(NewFile(curr_chain))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_chain))]).Normal] = get_dynare_random_generator_state();
save([BaseName '_mh' int2str(NewFile(curr_chain)) '_blck' int2str(curr_chain) '.mat'],'x2','logpo2','LastSeeds');
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,['\n']); fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_chain)) 'Blck' int2str(curr_chain) ' (' datestr(now,0) ')\n']); fprintf(fidlog,['%% Mh' int2str(NewFile(curr_chain)) 'Blck' int2str(curr_chain) ' (' datestr(now,0) ')\n']);

49
matlab/compute_trend_coefficients.m Normal file
View File

@ -0,0 +1,49 @@
function [trend_addition, trend_coeff]=compute_trend_coefficients(M_,DynareOptions,nvarobs,ntobs)
% [trend_addition, trend_coeff]=compute_trend_coefficients(DynareOptions,nvarobs,ntobs)
% Computes the trend coefficiencts and the trend, accounting for
% prefiltering
%
% INPUTS
% M_ [structure] describing the model; called in the eval
% statement
% DynareOptions [structure] describing the options
% nvarobs [scalar] number of observed variables
% ntobs [scalar] length of data sample for estimation
%
% OUTPUTS
% trend_addition [nvarobs by ntobs double] matrix storing deterministic
% trend component
% trend_coeff [nvarobs by 1] vector storing trend slope
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2014 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/>.
trend_coeff = zeros(nvarobs,1);
t = DynareOptions.trend_coeffs;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = eval(t{i});
end
end
trend_addition=trend_coeff*[DynareOptions.first_obs:DynareOptions.first_obs+ntobs-1];
if DynareOptions.prefilter
trend_addition = bsxfun(@minus,trend_addition,mean(trend_addition,2));
end

47
matlab/convert_dyn_45_to_44.m
View File

@ -52,7 +52,7 @@ if isfield(oo_,'PointForecast')
end end
%% change HPD-fields back to row vectors %% change HPD-fields back to row vectors
if isfield(oo_.PointForecast,'HPDinf') if isfield(oo_,'PointForecast') && isfield(oo_.PointForecast,'HPDinf')
names=fieldnames(oo_.PointForecast.HPDinf); names=fieldnames(oo_.PointForecast.HPDinf);
for ii=1:length(names) for ii=1:length(names)
oo_.PointForecast.HPDinf.(names{ii})=oo_.PointForecast.HPDinf.(names{ii})'; oo_.PointForecast.HPDinf.(names{ii})=oo_.PointForecast.HPDinf.(names{ii})';
@ -60,7 +60,7 @@ if isfield(oo_.PointForecast,'HPDinf')
end end
end end
if isfield(oo_.MeanForecast,'HPDinf') if isfield(oo_,'MeanForecast') && isfield(oo_.MeanForecast,'HPDinf')
names=fieldnames(oo_.MeanForecast.HPDinf); names=fieldnames(oo_.MeanForecast.HPDinf);
for ii=1:length(names) for ii=1:length(names)
oo_.MeanForecast.HPDinf.(names{ii})=oo_.MeanForecast.HPDinf.(names{ii})'; oo_.MeanForecast.HPDinf.(names{ii})=oo_.MeanForecast.HPDinf.(names{ii})';
@ -68,7 +68,7 @@ if isfield(oo_.MeanForecast,'HPDinf')
end end
end end
if isfield(oo_.UpdatedVariables,'HPDinf') if isfield(oo_,'UpdatedVariables') && isfield(oo_.UpdatedVariables,'HPDinf')
names=fieldnames(oo_.UpdatedVariables.HPDinf); names=fieldnames(oo_.UpdatedVariables.HPDinf);
for ii=1:length(names) for ii=1:length(names)
oo_.UpdatedVariables.HPDinf.(names{ii})=oo_.UpdatedVariables.HPDinf.(names{ii})'; oo_.UpdatedVariables.HPDinf.(names{ii})=oo_.UpdatedVariables.HPDinf.(names{ii})';
@ -76,7 +76,7 @@ if isfield(oo_.UpdatedVariables,'HPDinf')
end end
end end
if isfield(oo_.SmoothedVariables,'HPDinf') if isfield(oo_,'SmoothedVariables') && isfield(oo_.SmoothedVariables,'HPDinf')
names=fieldnames(oo_.SmoothedVariables.HPDinf); names=fieldnames(oo_.SmoothedVariables.HPDinf);
for ii=1:length(names) for ii=1:length(names)
oo_.SmoothedVariables.HPDinf.(names{ii})=oo_.SmoothedVariables.HPDinf.(names{ii})'; oo_.SmoothedVariables.HPDinf.(names{ii})=oo_.SmoothedVariables.HPDinf.(names{ii})';
@ -84,7 +84,7 @@ if isfield(oo_.SmoothedVariables,'HPDinf')
end end
end end
if isfield(oo_.FilteredVariables,'HPDinf') if isfield(oo_,'FilteredVariables') && isfield(oo_.FilteredVariables,'HPDinf')
names=fieldnames(oo_.FilteredVariables.HPDinf); names=fieldnames(oo_.FilteredVariables.HPDinf);
for ii=1:length(names) for ii=1:length(names)
oo_.FilteredVariables.HPDinf.(names{ii})=oo_.FilteredVariables.HPDinf.(names{ii})'; oo_.FilteredVariables.HPDinf.(names{ii})=oo_.FilteredVariables.HPDinf.(names{ii})';
@ -92,7 +92,7 @@ if isfield(oo_.FilteredVariables,'HPDinf')
end end
end end
if isfield(oo_.SmoothedShocks,'HPDinf') if isfield(oo_,'SmoothedShocks') && isfield(oo_.SmoothedShocks,'HPDinf')
names=fieldnames(oo_.SmoothedShocks.HPDinf); names=fieldnames(oo_.SmoothedShocks.HPDinf);
for ii=1:length(names) for ii=1:length(names)
oo_.SmoothedShocks.HPDinf.(names{ii})=oo_.SmoothedShocks.HPDinf.(names{ii})'; oo_.SmoothedShocks.HPDinf.(names{ii})=oo_.SmoothedShocks.HPDinf.(names{ii})';
@ -100,13 +100,44 @@ if isfield(oo_.SmoothedShocks,'HPDinf')
end end
end end
%% padd classical filtered variables with redundant zeros %% subtract mean from classical Updated variables
if isfield(oo_,'UpdatedVariables')
names=fieldnames(oo_.UpdatedVariables);
for ii=1:length(names)
%make sure Bayesian fields are not affected
if ~strcmp(names{ii},'Mean') && ~strcmp(names{ii},'Median') && ~strcmp(names{ii},'deciles') ...
&& ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup')
current_var_index=find(strmatch(names{ii},deblank(M_.endo_names),'exact'));
if options_.loglinear == 1 %logged steady state must be used
constant_current_variable=log(oo_.dr.ys(current_var_index));
elseif options_.loglinear == 0 %unlogged steady state must be used
constant_current_variable=oo_.dr.ys(current_var_index);
end
oo_.UpdatedVariables.(names{ii})=oo_.UpdatedVariables.(names{ii})-constant_current_variable;
if isfield(oo_.Smoother,'Trend') && isfield(oo_.Smoother.Trend,names{ii})
oo_.UpdatedVariables.(names{ii})=oo_.UpdatedVariables.(names{ii})-oo_.Smoother.Trend.(names{ii});
end
end
end
end
%% padd classical filtered variables with redundant zeros and subtract mean
if isfield(oo_,'FilteredVariables') if isfield(oo_,'FilteredVariables')
names=fieldnames(oo_.FilteredVariables); names=fieldnames(oo_.FilteredVariables);
for ii=1:length(names) for ii=1:length(names)
%make sure Bayesian fields are not affect %make sure Bayesian fields are not affected
if ~strcmp(names{ii},'Mean') && ~strcmp(names{ii},'Median') && ~strcmp(names{ii},'deciles') ... if ~strcmp(names{ii},'Mean') && ~strcmp(names{ii},'Median') && ~strcmp(names{ii},'deciles') ...
&& ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup') && ~strcmp(names{ii},'Var') && ~strcmp(names{ii},'HPDinf') && ~strcmp(names{ii},'HPDsup')
current_var_index=find(strmatch(names{ii},deblank(M_.endo_names),'exact'));
if options_.loglinear == 1 %logged steady state must be used
constant_current_variable=log(oo_.dr.ys(current_var_index));
elseif options_.loglinear == 0 %unlogged steady state must be used
constant_current_variable=oo_.dr.ys(current_var_index);
end
oo_.FilteredVariables.(names{ii})=oo_.FilteredVariables.(names{ii})-constant_current_variable;
if isfield(oo_.Smoother,'Trend') && isfield(oo_.Smoother.Trend,names{ii})
oo_.FilteredVariables.(names{ii})=oo_.FilteredVariables.(names{ii})-oo_.Smoother.Trend.(names{ii});
end
oo_.FilteredVariables.(names{ii})=[0; oo_.FilteredVariables.(names{ii}); zeros(options_.nk-1,1)]; oo_.FilteredVariables.(names{ii})=[0; oo_.FilteredVariables.(names{ii}); zeros(options_.nk-1,1)];
end end
end end

14
matlab/dsge_likelihood.m
View File

@ -307,14 +307,8 @@ end
% Define the deterministic linear trend of the measurement equation. % Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend if BayesInfo.with_trend
trend_coeff = zeros(DynareDataset.vobs,1); [trend_addition, trend_coeff]=compute_trend_coefficients(Model,DynareOptions,DynareDataset.vobs,DynareDataset.nobs);
t = DynareOptions.trend_coeffs; trend = repmat(constant,1,DynareDataset.nobs)+trend_addition;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = evalin('base',t{i});
end
end
trend = repmat(constant,1,DynareDataset.nobs)+trend_coeff*[1:DynareDataset.nobs];
else else
trend_coeff = zeros(DynareDataset.vobs,1); trend_coeff = zeros(DynareDataset.vobs,1);
trend = repmat(constant,1,DynareDataset.nobs); trend = repmat(constant,1,DynareDataset.nobs);
@ -496,7 +490,7 @@ switch DynareOptions.lik_init
indx_unstable = find(sum(abs(V),2)>1e-5); indx_unstable = find(sum(abs(V),2)>1e-5);
stable = find(sum(abs(V),2)<1e-5); stable = find(sum(abs(V),2)<1e-5);
nunit = length(eigenv) - nstable; nunit = length(eigenv) - nstable;
Pstar = options_.Harvey_scale_factor*eye(np); Pstar = DynareOptions.Harvey_scale_factor*eye(nunit);
if kalman_algo ~= 2 if kalman_algo ~= 2
kalman_algo = 1; kalman_algo = 1;
end end
@ -513,6 +507,8 @@ switch DynareOptions.lik_init
end end
Pstar(stable, stable) = Pstar_tmp; Pstar(stable, stable) = Pstar_tmp;
Pinf = []; Pinf = [];
a = zeros(mm,1);
Zflag = 0;
otherwise otherwise
error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!') error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!')
end end

57
matlab/dyn_forecast.m
View File

@ -1,12 +1,17 @@
function [forecast,info] = dyn_forecast(var_list,M,options,oo,task) function [forecast,info] = dyn_forecast(var_list,M,options,oo,task,dataset_info)
% function dyn_forecast(var_list,task) % function dyn_forecast(var_list,M,options,oo,task,dataset_info)
% computes mean forecast for a given value of the parameters % computes mean forecast for a given value of the parameters
% computes also confidence band for the forecast % compues also confidence band for the forecast
% %
% INPUTS % INPUTS
% var_list: list of variables (character matrix) % var_list: list of variables (character matrix)
% M: Dynare model structure
% options: Dynare options structure
% oo: Dynare results structure
% task: indicates how to initialize the forecast % task: indicates how to initialize the forecast
% either 'simul' or 'smoother' % either 'simul' or 'smoother'
% dataset_info: Various informations about the dataset (descriptive statistics and missing observations).
% OUTPUTS % OUTPUTS
% nothing is returned but the procedure saves output % nothing is returned but the procedure saves output
% in oo_.forecast.Mean % in oo_.forecast.Mean
@ -16,7 +21,7 @@ function [forecast,info] = dyn_forecast(var_list,M,options,oo,task)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2003-2015 Dynare Team % Copyright (C) 2003-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -33,9 +38,15 @@ function [forecast,info] = dyn_forecast(var_list,M,options,oo,task)
% 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/>.
if nargin<6 && options.prefilter
error('The prefiltering option is not allowed without providing a dataset')
elseif nargin==6
mean_varobs=dataset_info.descriptive.mean';
end
info = 0; info = 0;
make_ex_; oo=make_ex_(M,options,oo);
maximum_lag = M.maximum_lag; maximum_lag = M.maximum_lag;
@ -72,7 +83,20 @@ switch task
y0 = zeros(M.endo_nbr,maximum_lag); y0 = zeros(M.endo_nbr,maximum_lag);
for i = 1:M.endo_nbr for i = 1:M.endo_nbr
v_name = deblank(M.endo_names(i,:)); v_name = deblank(M.endo_names(i,:));
y0(i,:) = y_smoothed.(v_name)(end-maximum_lag+1:end)+oo.dr.ys(i); y0(i,:) = y_smoothed.(v_name)(end-maximum_lag+1:end); %includes steady state or mean, but simult_ will subtract only steady state
% 2. Subtract mean/steady state and add steady state; takes care of prefiltering
if isfield(oo.Smoother,'Constant') && isfield(oo.Smoother.Constant,v_name)
y0(i,:)=y0(i,:)-oo.Smoother.Constant.(v_name)(end-maximum_lag+1:end); %subtract mean or steady state
if options.loglinear
y0(i,:)=y0(i,:)+log(oo.dr.ys(strmatch(v_name,deblank(M.endo_names),'exact')));
else
y0(i,:)=y0(i,:)+oo.dr.ys(strmatch(v_name,deblank(M.endo_names),'exact'));
end
end
% 2. Subtract trend
if isfield(oo.Smoother,'Trend') && isfield(oo.Smoother.Trend,v_name)
y0(i,:)=y0(i,:)-oo.Smoother.Trend.(v_name)(end-maximum_lag+1:end); %subtract trend, which is not subtracted by simult_
end
end end
gend = options.nobs; gend = options.nobs;
if isfield(oo.Smoother,'TrendCoeffs') if isfield(oo.Smoother,'TrendCoeffs')
@ -89,12 +113,13 @@ switch task
end end
end end
if ~isempty(trend_coeffs) if ~isempty(trend_coeffs)
trend = trend_coeffs*(gend+(1-M.maximum_lag:horizon)); trend = trend_coeffs*(options.first_obs+gend-1+(1-M.maximum_lag:horizon));
if options.prefilter
trend = trend - repmat(mean(trend_coeffs*[options.first_obs:options.first_obs+gend-1],2),1,horizon+1); %subtract mean trend
end end
end end
global bayestopt_ else
if isfield(bayestopt_,'mean_varobs') trend_coeffs=zeros(length(options.varobs),1);
trend = trend + repmat(bayestopt_.mean_varobs,1,horizon+M.maximum_lag);
end end
otherwise otherwise
error('Wrong flag value') error('Wrong flag value')
@ -117,10 +142,20 @@ else
options.order,var_list,M,oo,options); options.order,var_list,M,oo,options);
end end
if ~isscalar(trend) if ~isscalar(trend) %add trend back to forecast
yf(i_var_obs,:) = yf(i_var_obs,:) + trend; yf(i_var_obs,:) = yf(i_var_obs,:) + trend;
end end
if options.loglinear == 1
if options.prefilter == 1 %subtract steady state and add mean for observables
yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(log(oo.dr.ys(i_var_obs)),1,horizon+M.maximum_lag)+ repmat(mean_varobs,1,horizon+M.maximum_lag);
end
else
if options.prefilter == 1 %subtract steady state and add mean for observables
yf(i_var_obs,:)=yf(i_var_obs,:)-repmat(oo.dr.ys(i_var_obs),1,horizon+M.maximum_lag)+ repmat(mean_varobs,1,horizon+M.maximum_lag);
end
end
for i=1:n_var for i=1:n_var
vname = deblank(var_list(i,:)); vname = deblank(var_list(i,:));
forecast.Mean.(vname) = yf(i,maximum_lag+(1:horizon))'; forecast.Mean.(vname) = yf(i,maximum_lag+(1:horizon))';

1
matlab/dynare_estimation.m
View File

@ -73,6 +73,7 @@ end
if options_.logged_steady_state if options_.logged_steady_state
oo_.dr.ys=exp(oo_.dr.ys); oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state); oo_.steady_state=exp(oo_.steady_state);
options_.logged_steady_state=0;
end end

129
matlab/dynare_estimation_1.m
View File

@ -170,38 +170,8 @@ end
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0 if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
if options_.smoother == 1 if options_.smoother == 1
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value); [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value);
oo_.Smoother.SteadyState = ys; [oo_]=write_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
oo_.Smoother.TrendCoeffs = trend_coeff;
if options_.filter_covariance
oo_.Smoother.Variance = P;
end
i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0
oo_.FilteredVariablesKStepAhead = ...
aK(options_.filter_step_ahead,i_endo,:);
if ~isempty(PK)
oo_.FilteredVariablesKStepAheadVariances = ...
PK(options_.filter_step_ahead,i_endo,i_endo,:);
end
if ~isempty(decomp)
oo_.FilteredVariablesShockDecomposition = ...
decomp(options_.filter_step_ahead,i_endo,:,:);
end
end
for i=bayestopt_.smoother_saved_var_list'
i1 = dr.order_var(bayestopt_.smoother_var_list(i));
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ...
' = atT(i,:)'';']);
if options_.nk > 0
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ...
' = squeeze(aK(1,i,2:end-(options_.nk-1)));']);
end
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ' = updated_variables(i,:)'';']);
end
for i=1:M_.exo_nbr
eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
end
end end
return return
end end
@ -513,39 +483,42 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
> 0) && options_.load_mh_file)) ... > 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && options_.partial_information == 0 % to be fixed || ~options_.smoother ) && options_.partial_information == 0 % to be fixed
%% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable %% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
oo_.Smoother.SteadyState = ys; [oo_,yf]=write_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
oo_.Smoother.TrendCoeffs = trend_coeff;
oo_.Smoother.Variance = P; % oo_.Smoother.SteadyState = ys;
i_endo = bayestopt_.smoother_saved_var_list; % oo_.Smoother.TrendCoeffs = trend_coeff;
if ~isempty(options_.nk) && options_.nk ~= 0 && (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))) % oo_.Smoother.Trend = Trend;
oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ... % oo_.Smoother.Variance = P;
i_endo,:); % i_endo = bayestopt_.smoother_saved_var_list;
if isfield(options_,'kalman_algo') % if ~isempty(options_.nk) && options_.nk ~= 0 && (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)))
if ~isempty(PK) % oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
oo_.FilteredVariablesKStepAheadVariances = ... % i_endo,:);
PK(options_.filter_step_ahead,i_endo,i_endo,:); % if isfield(options_,'kalman_algo')
end % if ~isempty(PK)
if ~isempty(decomp) % oo_.FilteredVariablesKStepAheadVariances = ...
oo_.FilteredVariablesShockDecomposition = ... % PK(options_.filter_step_ahead,i_endo,i_endo,:);
decomp(options_.filter_step_ahead,i_endo,:,:); % end
end % if ~isempty(decomp)
end % oo_.FilteredVariablesShockDecomposition = ...
end % decomp(options_.filter_step_ahead,i_endo,:,:);
for i=bayestopt_.smoother_saved_var_list' % end
i1 = dr.order_var(bayestopt_.smoother_var_list(i)); % end
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ' = ' ... % end
'atT(i,:)'';']); % for i=bayestopt_.smoother_saved_var_list'
if ~isempty(options_.nk) && options_.nk > 0 && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)) % i1 = dr.order_var(bayestopt_.smoother_var_list(i));
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ... % eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ' = ' ...
' = squeeze(aK(1,i,2:end-(options_.nk-1)));']); % 'atT(i,:)'';']);
end % if ~isempty(options_.nk) && options_.nk > 0 && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ... % eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ...
' = updated_variables(i,:)'';']); % ' = squeeze(aK(1,i,2:end-(options_.nk-1)));']);
end % end
for i=1:M_.exo_nbr % eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ...
eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']); % ' = updated_variables(i,:)'';']);
end % end
% for i=1:M_.exo_nbr
% eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
% end
if ~options_.nograph, if ~options_.nograph,
[nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr); [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
@ -616,18 +589,16 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fclose(fidTeX); fclose(fidTeX);
end end
end end
%% % %%
%% Smooth observational errors... % %% Smooth observational errors...
%% % %%
if options_.prefilter == 1 % if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
yf = atT(bayestopt_.mf,:)+repmat(dataset_info.descriptive.mean',1,gend); % yf = atT(bayestopt_.mf,:)+repmat(dataset_info.descriptive.mean',1,gend)+Trend;
elseif options_.loglinear == 1 % elseif options_.loglinear == 1
yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+... % yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+Trend;
trend_coeff*[1:gend]; % else
else % yf = atT(bayestopt_.mf,:)+repmat(ys(bayestopt_.mfys),1,gend)+Trend;
yf = atT(bayestopt_.mf,:)+repmat(ys(bayestopt_.mfys),1,gend)+... % end
trend_coeff*[1:gend];
end
if nvn if nvn
number_of_plots_to_draw = 0; number_of_plots_to_draw = 0;
index = []; index = [];
@ -636,7 +607,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
number_of_plots_to_draw = number_of_plots_to_draw + 1; number_of_plots_to_draw = number_of_plots_to_draw + 1;
index = cat(1,index,i); index = cat(1,index,i);
end end
eval(['oo_.SmoothedMeasurementErrors.' options_.varobs{i} ' = measurement_error(i,:)'';']); % eval(['oo_.SmoothedMeasurementErrors.' options_.varobs{i} ' = measurement_error(i,:)'';']);
end end
if ~options_.nograph if ~options_.nograph
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw); [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
@ -787,7 +758,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
end end
if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file
oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother'); oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
end end
if np > 0 if np > 0

14
matlab/dynare_estimation_init.m
View File

@ -375,13 +375,11 @@ if ~isfield(options_,'trend_coeffs') % No!
else% Yes! else% Yes!
bayestopt_.with_trend = 1; bayestopt_.with_trend = 1;
bayestopt_.trend_coeff = {}; bayestopt_.trend_coeff = {};
trend_coeffs = options_.trend_coeffs;
nt = length(trend_coeffs);
for i=1:options_.number_of_observed_variables for i=1:options_.number_of_observed_variables
if i > length(trend_coeffs) if i > length(options_.trend_coeffs)
bayestopt_.trend_coeff{i} = '0'; bayestopt_.trend_coeff{i} = '0';
else else
bayestopt_.trend_coeff{i} = trend_coeffs{i}; bayestopt_.trend_coeff{i} = options_.trend_coeffs{i};
end end
end end
end end
@ -530,9 +528,11 @@ end
[dataset_, dataset_info, newdatainterfaceflag] = makedataset(options_, options_.dsge_var*options_.dsge_varlag, gsa_flag); [dataset_, dataset_info, newdatainterfaceflag] = makedataset(options_, options_.dsge_var*options_.dsge_varlag, gsa_flag);
% Set options_.nobs %set options for old interface from the ones for new interface
if ~isempty(dataset_)
options_.nobs = dataset_.nobs; options_.nobs = dataset_.nobs;
options_.first_obs=double(dataset_.init);
end
% setting steadystate_check_flag option % setting steadystate_check_flag option
if options_.diffuse_filter if options_.diffuse_filter
steadystate_check_flag = 0; steadystate_check_flag = 0;
@ -556,7 +556,7 @@ if info(1)
print_info(info, 0, options_); print_info(info, 0, options_);
end end
if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9) if (~options_.loglinear && all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)) || (options_.loglinear && all(abs(log(oo_.steady_state(bayestopt_.mfys)))<1e-9))
options_.noconstant = 1; options_.noconstant = 1;
else else
options_.noconstant = 0; options_.noconstant = 0;

4
matlab/ep/extended_path.m
View File

@ -42,6 +42,10 @@ options_.simul.maxit = ep.maxit;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver % Prepare a structure needed by the matlab implementation of the perfect foresight model solver
pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_); pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_);
if M_.exo_det_nbr~=0
error('ep: Extended path does not support varexo_det.')
end
endo_nbr = M_.endo_nbr; endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr; exo_nbr = M_.exo_nbr;
maximum_lag = M_.maximum_lag; maximum_lag = M_.maximum_lag;

33
matlab/evaluate_smoother.m
View File

@ -87,35 +87,6 @@ if ischar(parameters)
end end
end end
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = ... [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend] = ...
DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state); DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
[oo_]=write_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;
if options_.filter_covariance
oo_.Smoother.Variance = P;
end
i_endo = bayestopt_.smoother_saved_var_list;
if options_.nk ~= 0
oo_.FilteredVariablesKStepAhead = ...
aK(options_.filter_step_ahead,i_endo,:);
if ~isempty(PK)
oo_.FilteredVariablesKStepAheadVariances = ...
PK(options_.filter_step_ahead,i_endo,i_endo,:);
end
if ~isempty(decomp)
oo_.FilteredVariablesShockDecomposition = ...
decomp(options_.filter_step_ahead,i_endo,:,:);
end
end
for i=bayestopt_.smoother_saved_var_list'
i1 = oo_.dr.order_var(bayestopt_.smoother_var_list(i));
eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ' = atT(i,:)'';']);
if options_.nk>0
eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ' = squeeze(aK(1,i,:));']);
end
eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ' = updated_variables(i,:)'';']);
end
for i=1:M_.exo_nbr
eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
end

2
matlab/evaluate_steady_state.m
View File

@ -60,7 +60,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
options); options);
%test whether it solves model conditional on the instruments %test whether it solves model conditional on the instruments
resids = evaluate_static_model(ys,exo_ss,params,M,options); resids = evaluate_static_model(ys,exo_ss,params,M,options);
n_multipliers=M.endo_nbr-M.orig_endo_nbr; n_multipliers=M.ramsey_eq_nbr;
nan_indices=find(isnan(resids(n_multipliers+1:end))); nan_indices=find(isnan(resids(n_multipliers+1:end)));
if ~isempty(nan_indices) if ~isempty(nan_indices)

4
matlab/forcst.m
View File

@ -17,7 +17,7 @@ function [yf,int_width]=forcst(dr,y0,horizon,var_list)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 2003-2012 Dynare Team % Copyright (C) 2003-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -36,7 +36,7 @@ function [yf,int_width]=forcst(dr,y0,horizon,var_list)
global M_ oo_ options_ global M_ oo_ options_
make_ex_; oo_=make_ex_(M_,options_,oo_);
yf = simult_(y0,dr,zeros(horizon,M_.exo_nbr),1); yf = simult_(y0,dr,zeros(horizon,M_.exo_nbr),1);
nstatic = M_.nstatic; nstatic = M_.nstatic;
nspred = M_.nspred; nspred = M_.nspred;

17
matlab/forcst2.m
View File

@ -1,5 +1,16 @@
function yf=forcst2(y0,horizon,dr,n) function yf=forcst2(y0,horizon,dr,n)
% function yf=forcst2(y0,horizon,dr,n)
%
% computes forecasts based on first order model solution, given shocks drawn from the shock distribution
% Inputs:
% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
% - dr [structure] structure with Dynare decision rules
% - e [horizon by exo_nbr] matrix with shocks
% - n [scalar] number of repetitions
%
% Outputs:
% - yf [horizon+ykmin_ by endo_nbr by n] array of forecasts
%
% Copyright (C) 2008-2010 Dynare Team % Copyright (C) 2008-2010 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
@ -24,14 +35,10 @@ endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr; exo_nbr = M_.exo_nbr;
ykmin_ = M_.maximum_endo_lag; ykmin_ = M_.maximum_endo_lag;
order = options_.order;
k1 = [ykmin_:-1:1]; k1 = [ykmin_:-1:1];
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]); k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr; k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
it_ = ykmin_ + 1 ;
% eliminate shocks with 0 variance % eliminate shocks with 0 variance
i_exo_var = setdiff([1:exo_nbr],find(diag(Sigma_e_) == 0)); i_exo_var = setdiff([1:exo_nbr],find(diag(Sigma_e_) == 0));
nxs = length(i_exo_var); nxs = length(i_exo_var);

16
matlab/forcst2a.m
View File

@ -1,6 +1,15 @@
function yf=forcst2a(y0,dr,e) function yf=forcst2a(y0,dr,e)
% function yf=forcst2a(y0,dr,e)
% Copyright (C) 2008-2010 Dynare Team % computes forecasts based on first order model solution, assuming the absence of shocks
% Inputs:
% - y0 [endo_nbr by maximum_endo_lag] matrix of starting values
% - dr [structure] structure with Dynare decision rules
% - e [horizon by exo_nbr] matrix with shocks
%
% Outputs:
% - yf [horizon+maximum_endo_lag,endo_nbr] matrix of forecasts
%
% Copyright (C) 2008-2015 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -19,13 +28,10 @@ function yf=forcst2a(y0,dr,e)
global M_ options_ global M_ options_
Sigma_e_ = M_.Sigma_e;
endo_nbr = M_.endo_nbr; endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
ykmin_ = M_.maximum_endo_lag; ykmin_ = M_.maximum_endo_lag;
horizon = size(e,1); horizon = size(e,1);
order = options_.order;
k1 = [ykmin_:-1:1]; k1 = [ykmin_:-1:1];
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]); k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);

5
matlab/gsa/read_data.m
View File

@ -24,7 +24,7 @@ function [gend, data] = read_data()
% 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_ bayestopt_ global options_
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
@ -36,8 +36,7 @@ if options_.loglinear == 1 & ~options_.logdata
rawdata = log(rawdata); rawdata = log(rawdata);
end end
if options_.prefilter == 1 if options_.prefilter == 1
bayestopt_.mean_varobs = mean(rawdata,1); data = transpose(rawdata-ones(gend,1)* mean(rawdata,1));
data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
else else
data = transpose(rawdata); data = transpose(rawdata);
end end

4
matlab/initial_estimation_checks.m
View File

@ -161,11 +161,13 @@ if info(1) > 0
print_info(info, DynareOptions.noprint, DynareOptions) print_info(info, DynareOptions.noprint, DynareOptions)
end end
if any(abs(DynareResults.steady_state(BayesInfo.mfys))>1e-9) && (DynareOptions.prefilter==1) if DynareOptions.prefilter==1
if (~DynareOptions.loglinear && any(abs(DynareResults.steady_state(BayesInfo.mfys))>1e-9)) || (DynareOptions.loglinear && any(abs(log(DynareResults.steady_state(BayesInfo.mfys)))>1e-9))
disp(['You are trying to estimate a model with a non zero steady state for the observed endogenous']) disp(['You are trying to estimate a model with a non zero steady state for the observed endogenous'])
disp(['variables using demeaned data!']) disp(['variables using demeaned data!'])
error('You should change something in your mod file...') error('You should change something in your mod file...')
end end
end
if ~isequal(DynareOptions.mode_compute,11) if ~isequal(DynareOptions.mode_compute,11)
disp(['Initial value of the log posterior (or likelihood): ' num2str(-fval)]); disp(['Initial value of the log posterior (or likelihood): ' num2str(-fval)]);

6
matlab/load_last_mh_history_file.m
View File

@ -10,7 +10,7 @@ function info = load_last_mh_history_file(MetropolisFolder, ModelName)
% Notes: The record structure is written to the caller workspace via an % Notes: The record structure is written to the caller workspace via an
% assignin statement. % assignin statement.
% Copyright (C) 2013 Dynare Team % Copyright (C) 2013-16 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -53,6 +53,7 @@ if format_mh_history_file %needed to preserve backward compatibility
record.AcceptanceRatio = record.AcceptationRates; record.AcceptanceRatio = record.AcceptationRates;
record.InitialSeeds = NaN; % This information is forever lost... record.InitialSeeds = NaN; % This information is forever lost...
record.MCMCConcludedSuccessfully = NaN; % This information is forever lost... record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...
record.MAX_nruns=NaN(size(record.MhDraws,1),1); % This information is forever lost...
record = rmfield(record,'LastLogLiK'); record = rmfield(record,'LastLogLiK');
record = rmfield(record,'InitialLogLiK'); record = rmfield(record,'InitialLogLiK');
record = rmfield(record,'Seeds'); record = rmfield(record,'Seeds');
@ -64,6 +65,9 @@ else
if ~isfield(record,'MCMCConcludedSuccessfully') if ~isfield(record,'MCMCConcludedSuccessfully')
record.MCMCConcludedSuccessfully = NaN; % This information is forever lost... record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...
end end
if ~isfield(record,'MAX_nruns')
record.MAX_nruns=NaN(size(record.MhDraws,1),1); % This information is forever lost...
end
end end
if isequal(nargout,0) if isequal(nargout,0)

225
matlab/metropolis_hastings_initialization.m
View File

@ -1,6 +1,6 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ...
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
%function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ... % function [ ix2, ilogpo2, ModelName, MetropolisFolder, FirstBlock, FirstLine, npar, NumberOfBlocks, nruns, NewFile, MAX_nruns, d ] = ...
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) % metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
% Metropolis-Hastings initialization. % Metropolis-Hastings initialization.
% %
@ -19,16 +19,16 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck,
% o oo_ outputs structure % o oo_ outputs structure
% %
% OUTPUTS % OUTPUTS
% o ix2 [double] (nblck*npar) vector of starting points for different chains % o ix2 [double] (NumberOfBlocks*npar) vector of starting points for different chains
% o ilogpo2 [double] (nblck*1) vector of initial posterior values for different chains % o ilogpo2 [double] (NumberOfBlocks*1) vector of initial posterior values for different chains
% o ModelName [string] name of the mod-file % o ModelName [string] name of the mod-file
% o MetropolisFolder [string] path to the Metropolis subfolder % o MetropolisFolder [string] path to the Metropolis subfolder
% o fblck [scalar] number of the first MH chain to be run (not equal to 1 in case of recovery) % o FirstBlock [scalar] number of the first MH chain to be run (not equal to 1 in case of recovery)
% o fline [double] (nblck*1) vector of first draw in each chain (not equal to 1 in case of recovery) % o FirstLine [double] (NumberOfBlocks*1) vector of first draw in each chain (not equal to 1 in case of recovery)
% o npar [scalar] number of parameters estimated % o npar [scalar] number of parameters estimated
% o nblck [scalar] Number of MCM chains requested % o NumberOfBlocks [scalar] Number of MCM chains requested
% o nruns [double] (nblck*1) number of draws in each chain % o nruns [double] (NumberOfBlocks*1) number of draws in each chain
% o NewFile [scalar] (nblck*1) vector storing the number % o NewFile [scalar] (NumberOfBlocks*1) vector storing the number
% of the first MH-file to created for each chain when saving additional % of the first MH-file to created for each chain when saving additional
% draws % draws
% o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk % o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk
@ -59,10 +59,10 @@ ix2 = [];
ilogpo2 = []; ilogpo2 = [];
ModelName = []; ModelName = [];
MetropolisFolder = []; MetropolisFolder = [];
fblck = []; FirstBlock = [];
fline = []; FirstLine = [];
npar = []; npar = [];
nblck = []; NumberOfBlocks = [];
nruns = []; nruns = [];
NewFile = []; NewFile = [];
MAX_nruns = []; MAX_nruns = [];
@ -76,15 +76,15 @@ end
MetropolisFolder = CheckPath('metropolis',M_.dname); MetropolisFolder = CheckPath('metropolis',M_.dname);
BaseName = [MetropolisFolder filesep ModelName]; BaseName = [MetropolisFolder filesep ModelName];
nblck = options_.mh_nblck; NumberOfBlocks = options_.mh_nblck;
nruns = ones(nblck,1)*options_.mh_replic; nruns = ones(NumberOfBlocks,1)*options_.mh_replic;
npar = length(xparam1); npar = length(xparam1);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
d = chol(vv); d = chol(vv);
if ~options_.load_mh_file && ~options_.mh_recover if ~options_.load_mh_file && ~options_.mh_recover
% Here we start a new Metropolis-Hastings, previous draws are discarded. % Here we start a new Metropolis-Hastings, previous draws are discarded.
if nblck > 1 if NumberOfBlocks > 1
disp('Estimation::mcmc: Multiple chains mode.') disp('Estimation::mcmc: Multiple chains mode.')
else else
disp('Estimation::mcmc: One Chain mode.') disp('Estimation::mcmc: One Chain mode.')
@ -108,16 +108,17 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n\n'); fprintf(fidlog,' \n\n');
fprintf(fidlog,'%% Session 1.\n'); fprintf(fidlog,'%% Session 1.\n');
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
fprintf(fidlog,[' Number of blocks...............: ' int2str(nblck) '\n']); fprintf(fidlog,[' Number of blocks...............: ' int2str(NumberOfBlocks) '\n']);
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']); fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
% Find initial values for the nblck chains... % Find initial values for the NumberOfBlocks chains...
if nblck > 1% Case 1: multiple chains if NumberOfBlocks > 1% Case 1: multiple chains
set_dynare_seed('default');
fprintf(fidlog,[' Initial values of the parameters:\n']); fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::mcmc: Searching for initial values...') disp('Estimation::mcmc: Searching for initial values...')
ix2 = zeros(nblck,npar); ix2 = zeros(NumberOfBlocks,npar);
ilogpo2 = zeros(nblck,1); ilogpo2 = zeros(NumberOfBlocks,1);
for j=1:nblck for j=1:NumberOfBlocks
validate = 0; validate = 0;
init_iter = 0; init_iter = 0;
trial = 1; trial = 1;
@ -183,22 +184,23 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
end end
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
fblck = 1; FirstBlock = 1;
fline = ones(nblck,1); FirstLine = ones(NumberOfBlocks,1);
NewFile = ones(nblck,1); NewFile = ones(NumberOfBlocks,1);
% Delete the mh-history files % Delete the mh-history files
delete_mh_history_files(MetropolisFolder, ModelName); delete_mh_history_files(MetropolisFolder, ModelName);
% Create a new record structure % Create a new record structure
fprintf(['Estimation::mcmc: Write details about the MCMC... ']); fprintf(['Estimation::mcmc: Write details about the MCMC... ']);
AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns); AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns; AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
record.Nblck = nblck; record.Nblck = NumberOfBlocks;
record.MhDraws = zeros(1,3); record.MhDraws = zeros(1,3);
record.MhDraws(1,1) = nruns(1); record.MhDraws(1,1) = nruns(1);
record.MhDraws(1,2) = AnticipatedNumberOfFiles; record.MhDraws(1,2) = AnticipatedNumberOfFiles;
record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile; record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
record.AcceptanceRatio = zeros(1,nblck); record.MAX_nruns=MAX_nruns;
for j=1:nblck record.AcceptanceRatio = zeros(1,NumberOfBlocks);
for j=1:NumberOfBlocks
% we set a different seed for the random generator for each block then we record the corresponding random generator state (vector) % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
set_dynare_seed(options_.DynareRandomStreams.seed+j); set_dynare_seed(options_.DynareRandomStreams.seed+j);
% record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name % record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
@ -206,8 +208,8 @@ if ~options_.load_mh_file && ~options_.mh_recover
end end
record.InitialParameters = ix2; record.InitialParameters = ix2;
record.InitialLogPost = ilogpo2; record.InitialLogPost = ilogpo2;
record.LastParameters = zeros(nblck,npar); record.LastParameters = zeros(NumberOfBlocks,npar);
record.LastLogPost = zeros(nblck,1); record.LastLogPost = zeros(NumberOfBlocks,1);
record.LastFileNumber = AnticipatedNumberOfFiles ; record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile; record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
record.MCMCConcludedSuccessfully = 0; record.MCMCConcludedSuccessfully = 0;
@ -219,7 +221,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']); fprintf(fidlog,[' Expected number of files per block.......: ' int2str(AnticipatedNumberOfFiles) '.\n']);
fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']); fprintf(fidlog,[' Expected number of lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\n']);
fprintf(fidlog,['\n']); fprintf(fidlog,['\n']);
for j = 1:nblck, for j = 1:NumberOfBlocks,
fprintf(fidlog,[' Initial state of the Gaussian random number generator for chain number ',int2str(j),':\n']); fprintf(fidlog,[' Initial state of the Gaussian random number generator for chain number ',int2str(j),':\n']);
for i=1:length(record.InitialSeeds(j).Normal) for i=1:length(record.InitialSeeds(j).Normal)
fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Normal(i)') '\n']); fprintf(fidlog,[' ' num2str(record.InitialSeeds(j).Normal(i)') '\n']);
@ -247,30 +249,30 @@ elseif options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,'\n'); fprintf(fidlog,'\n');
fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']); fprintf(fidlog,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\n']);
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
fprintf(fidlog,[' Number of blocks...............: ' int2str(nblck) '\n']); fprintf(fidlog,[' Number of blocks...............: ' int2str(NumberOfBlocks) '\n']);
fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']); fprintf(fidlog,[' Number of simulations per block: ' int2str(nruns(1)) '\n']);
fprintf(fidlog,' \n'); fprintf(fidlog,' \n');
past_number_of_blocks = record.Nblck; past_number_of_blocks = record.Nblck;
if past_number_of_blocks ~= nblck if past_number_of_blocks ~= NumberOfBlocks
disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!') disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!')
disp(['Estimation::mcmc: You declared ' int2str(nblck) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.']) disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ]) disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
nblck = past_number_of_blocks; NumberOfBlocks = past_number_of_blocks;
options_.mh_nblck = nblck; options_.mh_nblck = NumberOfBlocks;
end end
% I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations: % I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations:
LastFileNumber = record.LastFileNumber; LastFileNumber = record.LastFileNumber;
LastLineNumber = record.LastLineNumber; LastLineNumber = record.LastLineNumber;
if LastLineNumber < MAX_nruns if LastLineNumber < MAX_nruns
NewFile = ones(nblck,1)*LastFileNumber; NewFile = ones(NumberOfBlocks,1)*LastFileNumber;
fline = ones(nblck,1)*(LastLineNumber+1); FirstLine = ones(NumberOfBlocks,1)*(LastLineNumber+1);
else else
NewFile = ones(nblck,1)*(LastFileNumber+1); NewFile = ones(NumberOfBlocks,1)*(LastFileNumber+1);
fline = ones(nblck,1); FirstLine = ones(NumberOfBlocks,1);
end end
ilogpo2 = record.LastLogPost; ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters; ix2 = record.LastParameters;
fblck = 1; FirstBlock = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1); NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
fprintf('Estimation::mcmc: I am writing a new mh-history file... '); fprintf('Estimation::mcmc: I am writing a new mh-history file... ');
record.MhDraws = [record.MhDraws;zeros(1,3)]; record.MhDraws = [record.MhDraws;zeros(1,3)];
@ -283,6 +285,7 @@ elseif options_.load_mh_file && ~options_.mh_recover
record.MhDraws(end,1) = nruns(1); record.MhDraws(end,1) = nruns(1);
record.MhDraws(end,2) = AnticipatedNumberOfFiles; record.MhDraws(end,2) = AnticipatedNumberOfFiles;
record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile; record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
record.MAX_nruns = [record.MAX_nruns;MAX_nruns];
record.InitialSeeds = record.LastSeeds; record.InitialSeeds = record.LastSeeds;
write_mh_history_file(MetropolisFolder, ModelName, record); write_mh_history_file(MetropolisFolder, ModelName, record);
fprintf('Done.\n') fprintf('Done.\n')
@ -293,83 +296,143 @@ elseif options_.mh_recover
% The previous metropolis-hastings crashed before the end! I try to recover the saved draws... % The previous metropolis-hastings crashed before the end! I try to recover the saved draws...
disp('Estimation::mcmc: Recover mode!') disp('Estimation::mcmc: Recover mode!')
load_last_mh_history_file(MetropolisFolder, ModelName); load_last_mh_history_file(MetropolisFolder, ModelName);
nblck = record.Nblck;% Number of "parallel" mcmc chains. NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains.
options_.mh_nblck = nblck; options_.mh_nblck = NumberOfBlocks;
%% check consistency of options
if record.MhDraws(end,1)~=options_.mh_replic
fprintf('\nEstimation::mcmc: You cannot specify a different mh_replic than in the chain you are trying to recover\n')
fprintf('Estimation::mcmc: I am resetting mh_replic to %u\n',record.MhDraws(end,1))
options_.mh_replic=record.MhDraws(end,1);
nruns = ones(NumberOfBlocks,1)*options_.mh_replic;
end
if ~isnan(record.MAX_nruns(end,1)) %field exists
if record.MAX_nruns(end,1)~=MAX_nruns
fprintf('\nEstimation::mcmc: You cannot specify a different MaxNumberOfBytes than in the chain you are trying to recover\n')
fprintf('Estimation::mcmc: I am resetting MAX_nruns to %u\n',record.MAX_nruns(end,1))
MAX_nruns=record.MAX_nruns(end,1);
end
end
%% do tentative initialization as if full last run of MCMC were to be redone
if size(record.MhDraws,1) == 1 if size(record.MhDraws,1) == 1
OldMh = 0;% The crashed metropolis was the first session. OldMhExists = 0;% The crashed metropolis was the first session.
else else
OldMh = 1;% The crashed metropolis wasn't the first session. OldMhExists = 1;% The crashed metropolis wasn't the first session.
end end
% Default initialization: % Default initialization:
if OldMh if OldMhExists
ilogpo2 = record.LastLogPost; ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters; ix2 = record.LastParameters;
else else
ilogpo2 = record.InitialLogPost; ilogpo2 = record.InitialLogPost;
ix2 = record.InitialParameters; ix2 = record.InitialParameters;
end end
% Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers. % Set NewFile, a NumberOfBlocks*1 vector of integers, and FirstLine (first line), a NumberOfBlocks*1 vector of integers.
if OldMh % Relevant for starting at the correct file and potentially filling a file from the previous run that was not yet full
if OldMhExists
LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session. LastLineNumberInThePreviousMh = record.MhDraws(end-1,3);% Number of lines in the last mh files of the previous session.
LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions. LastFileNumberInThePreviousMh = sum(record.MhDraws(1:end-1,2),1);% Number of mh files in the the previous sessions.
if LastLineNumberInThePreviousMh < MAX_nruns% Test if the last mh files of the previous session are not complete (yes) %Test if the last mh files of the previous session were not full yet
NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh; if LastLineNumberInThePreviousMh < MAX_nruns%not full
fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1); %store starting point if whole chain needs to be redone
else% The last mh files of the previous session are complete. NewFile = ones(NumberOfBlocks,1)*LastFileNumberInThePreviousMh;
NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1); FirstLine = ones(NumberOfBlocks,1)*(LastLineNumberInThePreviousMh+1);
fline = ones(nblck,1); LastFileFullIndicator=0;
else% The last mh files of the previous session were full
%store starting point if whole chain needs to be redone
NewFile = ones(NumberOfBlocks,1)*(LastFileNumberInThePreviousMh+1);
FirstLine = ones(NumberOfBlocks,1);
LastFileFullIndicator=1;
end end
else else
LastLineNumberInThePreviousMh = 0; LastLineNumberInThePreviousMh = 0;
LastFileNumberInThePreviousMh = 0; LastFileNumberInThePreviousMh = 0;
NewFile = ones(nblck,1); NewFile = ones(NumberOfBlocks,1);
fline = ones(nblck,1); FirstLine = ones(NumberOfBlocks,1);
LastFileFullIndicator=1;
end end
% Set fblck (First block), an integer targeting the crashed mcmc chain.
fblck = 1; %% Now find out what exactly needs to be redone
% 1. Check if really something needs to be done
% How many mh files should we have ? % How many mh files should we have ?
ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1); ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck; ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*NumberOfBlocks;
% I count the total number of saved mh files... % How many mh files do we actually have ?
AllMhFiles = dir([BaseName '_mh*_blck*.mat']); AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
TotalNumberOfMhFiles = length(AllMhFiles); TotalNumberOfMhFiles = length(AllMhFiles);
% And I quit if I can't find a crashed mcmc chain. % Quit if no crashed mcmc chain can be found as there are as many files as expected
if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles) if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!') disp('Estimation::mcmc: It appears that you don''t need to use the mh_recover option!')
disp(' You have to edit the mod file and remove the mh_recover option') disp(' You have to edit the mod file and remove the mh_recover option')
disp(' in the estimation command') disp(' in the estimation command')
error('Estimation::mcmc: mh_recover option not required!') error('Estimation::mcmc: mh_recover option not required!')
end end
% I count the number of saved mh files per block. % 2. Something needs to be done; find out what
NumberOfMhFilesPerBlock = zeros(nblck,1); % Count the number of saved mh files per block.
for b = 1:nblck NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
for b = 1:NumberOfBlocks
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']); BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
NumberOfMhFilesPerBlock(b) = length(BlckMhFiles); NumberOfMhFilesPerBlock(b) = length(BlckMhFiles);
end end
% Is there a chain with less mh files than expected ? % Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
while fblck <= nblck FirstBlock = 1; %initialize
if NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock while FirstBlock <= NumberOfBlocks
disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is not complete!']) if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is not complete!'])
break break
% The mh_recover session will start from chain fblck. % The mh_recover session will start from chain FirstBlock.
else else
disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is complete!']) disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is complete!'])
end end
fblck = fblck+1; FirstBlock = FirstBlock+1;
end end
%% 3. Overwrite default settings for
% How many mh-files are saved in this block? % How many mh-files are saved in this block?
NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck); NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(FirstBlock);
% Correct the number of saved mh files if the crashed Metropolis was not the first session (so ExistingDrawsInLastMCFile=0; %initialize: no MCMC draws of current MCMC are in file from last run
% that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain % Check whether last present file is a file included in the last MCMC run
% of the current session). if ~LastFileFullIndicator
if OldMh if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only that last file exists, but no files from current MCMC
NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh; loaded_results=load([BaseName '_mh' int2str(NewFile(FirstBlock)) '_blck' int2str(FirstBlock) '.mat']);
%check whether that last file was filled
if size(loaded_results.x2,1)==MAX_nruns %file is full
NewFile(FirstBlock)=NewFile(FirstBlock)+1; %set first file to be created to next one
FirstLine(FirstBlock) = 1; %use first line of next file
ExistingDrawsInLastMCFile=MAX_nruns-record.MhDraws(end-1,3);
else
ExistingDrawsInLastMCFile=0;
end end
NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh; end
elseif LastFileFullIndicator
ExistingDrawsInLastMCFile=0;
if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only the last file exists, but no files from current MCMC
NewFile(FirstBlock)=NewFile(FirstBlock)+1; %set first file to be created to next one
end
end
% % Correct the number of saved mh files if the crashed Metropolis was not the first session (so
% % that NumberOfSavedMhFilesInTheCrashedBlck is the number of saved mh files in the crashed chain
% % of the current session).
% if OldMhExists
% NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
% end
% NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
% Correct initial conditions. % Correct initial conditions.
if NumberOfSavedMhFiles if NumberOfSavedMhFilesInTheCrashedBlck<ExpectedNumberOfMhFilesPerBlock
load([BaseName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']); loaded_results=load([BaseName '_mh' int2str(NumberOfSavedMhFilesInTheCrashedBlck) '_blck' int2str(FirstBlock) '.mat']);
ilogpo2(fblck) = logpo2(end); ilogpo2(FirstBlock) = loaded_results.logpo2(end);
ix2(fblck,:) = x2(end,:); ix2(FirstBlock,:) = loaded_results.x2(end,:);
nruns(FirstBlock)=nruns(FirstBlock)-ExistingDrawsInLastMCFile-(NumberOfSavedMhFilesInTheCrashedBlck-LastFileNumberInThePreviousMh)*MAX_nruns;
%reset seed if possible
if isfield(loaded_results,'LastSeeds')
record.InitialSeeds(FirstBlock).Unifor=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Unifor;
record.InitialSeeds(FirstBlock).Normal=loaded_results.LastSeeds.(['file' int2str(NumberOfSavedMhFilesInTheCrashedBlck)]).Normal;
else
fprintf('Estimation::mcmc: You are trying to recover a chain generated with an older Dynare version.\n')
fprintf('Estimation::mcmc: I am using the default seeds to continue the chain.\n')
end
write_mh_history_file(MetropolisFolder, ModelName, record);
end end
end end

15
matlab/non_linear_dsge_likelihood.m
View File

@ -227,6 +227,21 @@ end
% Define a vector of indices for the observed variables. Is this really usefull?... % Define a vector of indices for the observed variables. Is this really usefull?...
BayesInfo.mf = BayesInfo.mf1; BayesInfo.mf = BayesInfo.mf1;
% Define the deterministic linear trend of the measurement equation.
if DynareOptions.noconstant
constant = zeros(DynareDataset.vobs,1);
else
constant = SteadyState(BayesInfo.mfys);
end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
[trend_addition, trend_coeff]=compute_trend_coefficients(Model,DynareOptions,DynareDataset.vobs,DynareDataset.nobs);
trend = repmat(constant,1,DynareDataset.info.ntobs)+trend_addition;
else
trend = repmat(constant,1,DynareDataset.nobs);
end
% Get needed informations for kalman filter routines. % Get needed informations for kalman filter routines.
start = DynareOptions.presample+1; start = DynareOptions.presample+1;
np = size(T,1); np = size(T,1);

4
matlab/osr.m
View File

@ -21,7 +21,7 @@ function osr_res = osr(var_list,params,i_var,W)
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none. % none.
% %
% Copyright (C) 2001-2012 Dynare Team % Copyright (C) 2001-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -46,7 +46,7 @@ if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6; options_.qz_criterium = 1+1e-6;
end end
make_ex_; oo_=make_ex_(M_,options_,oo_);
np = size(params,1); np = size(params,1);
i_params = zeros(np,1); i_params = zeros(np,1);

10
matlab/perfect-foresight-models/det_cond_forecast.m
View File

@ -12,7 +12,7 @@ function data_set = det_cond_forecast(varargin)
% dataset [dseries] Returns a dseries containing the forecasted endgenous variables and shocks % dataset [dseries] Returns a dseries containing the forecasted endgenous variables and shocks
% %
% %
% Copyright (C) 2013-2014 Dynare Team % Copyright (C) 2013-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -412,8 +412,8 @@ if pf && ~surprise
indx_endo(col_count : col_count + constrained_periods - 1) = constrained_vars(j) + (time_index_constraint - 1) * ny; indx_endo(col_count : col_count + constrained_periods - 1) = constrained_vars(j) + (time_index_constraint - 1) * ny;
col_count = col_count + constrained_periods; col_count = col_count + constrained_periods;
end; end;
make_ex_; oo_=make_ex_(M_,options_,oo_);
make_y_; oo_=make_y_(M_,options_,oo_);
it = 1; it = 1;
convg = 0; convg = 0;
normra = 1e+50; normra = 1e+50;
@ -613,13 +613,13 @@ else
disp(['t=' int2str(t) ' conditional (surprise=' int2str(surprise) ' perfect foresight=' int2str(pf) ') unconditional (surprise=' int2str(b_surprise) ' perfect foresight=' int2str(b_pf) ')']); disp(['t=' int2str(t) ' conditional (surprise=' int2str(surprise) ' perfect foresight=' int2str(pf) ') unconditional (surprise=' int2str(b_surprise) ' perfect foresight=' int2str(b_pf) ')']);
disp('==============================================================================================='); disp('===============================================================================================');
if t == 1 if t == 1
make_ex_; oo_=make_ex_(M_,options_,oo_);
if maximum_lag > 0 if maximum_lag > 0
exo_init = oo_.exo_simul; exo_init = oo_.exo_simul;
else else
exo_init = zeros(size(oo_.exo_simul)); exo_init = zeros(size(oo_.exo_simul));
end end
make_y_; oo_=make_y_(M_,options_,oo_);
end; end;
%exo_init %exo_init
oo_.exo_simul = exo_init; oo_.exo_simul = exo_init;

12
matlab/perfect-foresight-models/make_ex_.m
View File

@ -1,18 +1,20 @@
function make_ex_() function oo_=make_ex_(M_,options_,oo_)
% forms oo_.exo_simul and oo_.exo_det_simul % forms oo_.exo_simul and oo_.exo_det_simul
% %
% INPUTS % INPUTS
% none % M_: Dynare model structure
% options_: Dynare options structure
% oo_: Dynare results structure
% %
% OUTPUTS % OUTPUTS
% none % oo_: Dynare results structure
% %
% ALGORITHM % ALGORITHM
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% %
% Copyright (C) 1996-2014 Dynare Team % Copyright (C) 1996-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -29,7 +31,7 @@ function make_ex_()
% 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_ ex0_ global ex0_
if isempty(oo_.exo_steady_state) if isempty(oo_.exo_steady_state)
oo_.exo_steady_state = zeros(M_.exo_nbr,1); oo_.exo_steady_state = zeros(M_.exo_nbr,1);

16
matlab/perfect-foresight-models/make_y_.m
View File

@ -1,18 +1,22 @@
function make_y_() function oo_=make_y_(M_,options_,oo_)
% function make_y_ % function oo_=make_y_(M_,options_,oo_)
% forms oo_.endo_simul as guess values for deterministic simulations % forms oo_.endo_simul as guess values for deterministic simulations
% %
% INPUTS % INPUTS
% ... % M_: Dynare model structure
% options_: Dynare options structure
% oo_: Dynare results structure
%
% OUTPUTS % OUTPUTS
% ... % oo_: Dynare results structure
%
% ALGORITHM % ALGORITHM
% ... % ...
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% %
% Copyright (C) 1996-2014 Dynare Team % Copyright (C) 1996-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -29,7 +33,7 @@ function make_y_()
% 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_ ys0_ global ys0_
if options_.steadystate_flag if options_.steadystate_flag
[oo_.steady_state,M_.params,check] = ... [oo_.steady_state,M_.params,check] = ...

6
matlab/perfect-foresight-models/perfect_foresight_setup.m
View File

@ -12,7 +12,7 @@ function perfect_foresight_setup()
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 1996-2014 Dynare Team % Copyright (C) 1996-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -53,8 +53,8 @@ end
if ~options_.initval_file if ~options_.initval_file
if isempty(options_.datafile) if isempty(options_.datafile)
make_ex_; oo_=make_ex_(M_,options_,oo_);
make_y_; oo_=make_y_(M_,options_,oo_);
else else
read_data_; read_data_;
end end

32
matlab/perfect-foresight-models/perfect_foresight_solver.m
View File

@ -12,7 +12,7 @@ function perfect_foresight_solver()
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% none % none
% Copyright (C) 1996-2015 Dynare Team % Copyright (C) 1996-2016 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -62,15 +62,17 @@ oo_ = perfect_foresight_solver_core(M_,options_,oo_);
% If simulation failed try homotopy. % If simulation failed try homotopy.
if ~oo_.deterministic_simulation.status && ~options_.no_homotopy if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
skipline() skipline()
disp('Simulation of the perfect foresight model failed!') disp('Simulation of the perfect foresight model failed!')
disp('Switching to a homotopy method...') disp('Switching to a homotopy method...')
skipline()
if ~M_.maximum_lag if ~M_.maximum_lag
disp('Homotopy not implemented for purely forward models!') disp('Homotopy not implemented for purely forward models!')
disp('Failed to solve the model!') disp('Failed to solve the model!')
disp('Return with empty oo_.endo_simul.') disp('Return with empty oo_.endo_simul.')
oo_.endo_simul = []; oo_.endo_simul = [];
skipline()
return return
end end
if ~M_.maximum_lead if ~M_.maximum_lead
@ -78,10 +80,9 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
disp('Failed to solve the model!') disp('Failed to solve the model!')
disp('Return with empty oo_.endo_simul.') disp('Return with empty oo_.endo_simul.')
oo_.endo_simul = []; oo_.endo_simul = [];
skipline()
return return
end end
skipline()
% Disable warnings if homotopy % Disable warnings if homotopy
warning_old_state = warning; warning_old_state = warning;
warning off all warning off all
@ -89,14 +90,16 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
oldverbositylevel = options_.verbosity; oldverbositylevel = options_.verbosity;
options_.verbosity = 0; options_.verbosity = 0;
exosim = oo_.exo_simul; % Set initial paths for the endogenous and exogenous variables.
endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead);
exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1); exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
% Copy the current paths for the exogenous and endogenous variables.
exosim = oo_.exo_simul;
endosim = oo_.endo_simul; endosim = oo_.endo_simul;
endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead);
current_weight = 0; % Current weight of target point in convex combination current_weight = 0; % Current weight of target point in convex combination.
step = .5; step = .5; % Set default step size.
success_counter = 0; success_counter = 0;
iteration = 0; iteration = 0;
@ -120,21 +123,24 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
% Compute convex combination for exo path and initial/terminal endo conditions % Compute convex combination for exo path and initial/terminal endo conditions
% But take care of not overwriting the computed part of oo_.endo_simul % But take care of not overwriting the computed part of oo_.endo_simul
oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight); oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight);
oo_.endo_simul(:,[initperiods, lastperiods]) = new_weight*endosim(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
oo_.endo_simul(:,[initperiods, lastperiods]) = ... % Detect Nans or complex numbers in the solution.
new_weight*oo_.endo_simul(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
path_with_nans = any(any(isnan(oo_.endo_simul))); path_with_nans = any(any(isnan(oo_.endo_simul)));
path_with_cplx = any(any(~isreal(oo_.endo_simul))); path_with_cplx = any(any(~isreal(oo_.endo_simul)));
if isequal(iteration, 1) if isequal(iteration, 1)
% First iteration, same initial guess as in the first call to perfect_foresight_solver_core routine.
oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = endoinit(:,1:options_.periods); oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = endoinit(:,1:options_.periods);
elseif path_with_nans || path_with_cplx elseif path_with_nans || path_with_cplx
% If solver failed with NaNs or complex number, use previous solution as an initial guess.
oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = saved_endo_simul(:,1+M_.maximum_lag:end-M_.maximum_lead); oo_.endo_simul(:,M_.maximum_lag+1:end-M_.maximum_lead) = saved_endo_simul(:,1+M_.maximum_lag:end-M_.maximum_lead);
end end
% Make a copy of the paths.
saved_endo_simul = oo_.endo_simul; saved_endo_simul = oo_.endo_simul;
% Solve for the paths of the endogenous variables.
[oo_,me] = perfect_foresight_solver_core(M_,options_,oo_); [oo_,me] = perfect_foresight_solver_core(M_,options_,oo_);
if oo_.deterministic_simulation.status == 1 if oo_.deterministic_simulation.status == 1
@ -150,6 +156,7 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
end end
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me) fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me)
else else
% If solver failed, then go back.
oo_.endo_simul = saved_endo_simul; oo_.endo_simul = saved_endo_simul;
success_counter = 0; success_counter = 0;
step = step / 2; step = step / 2;
@ -168,11 +175,12 @@ end
if oo_.deterministic_simulation.status == 1 if oo_.deterministic_simulation.status == 1
disp('Perfect foresight solution found.') disp('Perfect foresight solution found.')
skipline()
else else
warning('Failed to solve perfect foresight model') warning('Failed to solve perfect foresight model')
end end
skipline()
dyn2vec; dyn2vec;
if ~isdates(options_.initial_period) && isnan(options_.initial_period) if ~isdates(options_.initial_period) && isnan(options_.initial_period)

9
matlab/prior_posterior_statistics.m
View File

@ -44,7 +44,8 @@ Y = transpose(dataset.data);
gend = dataset.nobs; gend = dataset.nobs;
data_index = dataset_info.missing.aindex; data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state; missing_value = dataset_info.missing.state;
bayestopt_.mean_varobs = dataset_info.descriptive.mean'; mean_varobs = dataset_info.descriptive.mean;
nvx = estim_params_.nvx; nvx = estim_params_.nvx;
nvn = estim_params_.nvn; nvn = estim_params_.nvn;
@ -63,6 +64,7 @@ nvobs = length(options_.varobs);
iendo = 1:endo_nbr; iendo = 1:endo_nbr;
horizon = options_.forecast; horizon = options_.forecast;
filtered_vars = options_.filtered_vars; filtered_vars = options_.filtered_vars;
IdObs = bayestopt_.mfys;
if horizon if horizon
i_last_obs = gend+(1-M_.maximum_endo_lag:0); i_last_obs = gend+(1-M_.maximum_endo_lag:0);
end end
@ -106,8 +108,6 @@ if horizon
MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8)); MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ... MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
8)); 8));
IdObs = bayestopt_.mfys;
end end
MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8))); MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
@ -161,15 +161,16 @@ localVars.Y=Y;
localVars.data_index=data_index; localVars.data_index=data_index;
localVars.missing_value=missing_value; localVars.missing_value=missing_value;
localVars.varobs=options_.varobs; localVars.varobs=options_.varobs;
localVars.mean_varobs=mean_varobs;
localVars.irun=irun; localVars.irun=irun;
localVars.endo_nbr=endo_nbr; localVars.endo_nbr=endo_nbr;
localVars.nvn=nvn; localVars.nvn=nvn;
localVars.naK=naK; localVars.naK=naK;
localVars.horizon=horizon; localVars.horizon=horizon;
localVars.iendo=iendo; localVars.iendo=iendo;
localVars.IdObs=IdObs;
if horizon if horizon
localVars.i_last_obs=i_last_obs; localVars.i_last_obs=i_last_obs;
localVars.IdObs=IdObs;
localVars.MAX_nforc1=MAX_nforc1; localVars.MAX_nforc1=MAX_nforc1;
localVars.MAX_nforc2=MAX_nforc2; localVars.MAX_nforc2=MAX_nforc2;
end end

62
matlab/prior_posterior_statistics_core.m
View File

@ -57,15 +57,16 @@ Y=myinputs.Y;
data_index=myinputs.data_index; data_index=myinputs.data_index;
missing_value=myinputs.missing_value; missing_value=myinputs.missing_value;
varobs=myinputs.varobs; varobs=myinputs.varobs;
mean_varobs=myinputs.mean_varobs;
irun=myinputs.irun; irun=myinputs.irun;
endo_nbr=myinputs.endo_nbr; endo_nbr=myinputs.endo_nbr;
nvn=myinputs.nvn; nvn=myinputs.nvn;
naK=myinputs.naK; naK=myinputs.naK;
horizon=myinputs.horizon; horizon=myinputs.horizon;
iendo=myinputs.iendo; iendo=myinputs.iendo;
IdObs=myinputs.IdObs; %index of observables
if horizon if horizon
i_last_obs=myinputs.i_last_obs; i_last_obs=myinputs.i_last_obs;
IdObs=myinputs.IdObs;
MAX_nforc1=myinputs.MAX_nforc1; MAX_nforc1=myinputs.MAX_nforc1;
MAX_nforc2=myinputs.MAX_nforc2; MAX_nforc2=myinputs.MAX_nforc2;
end end
@ -166,10 +167,10 @@ for b=fpar:B
if run_smoother if run_smoother
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK] = ... [alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,junk3,junk4,junk5,trend_addition] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value); DsgeSmoother(deep,gend,Y,data_index,missing_value);
if options_.loglinear if options_.loglinear %reads values from smoother results, which are in dr-order and put them into declaration order
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ... stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(log(SteadyState(dr.order_var)),1,gend); repmat(log(SteadyState(dr.order_var)),1,gend);
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ... stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
@ -180,22 +181,61 @@ for b=fpar:B
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ... stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
repmat(SteadyState(dr.order_var),1,gend); repmat(SteadyState(dr.order_var),1,gend);
end end
%% Compute constant for observables
if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
constant_part=repmat(mean_varobs',1,gend);
elseif options_.prefilter == 0 && options_.loglinear == 1 %logged steady state must be used
constant_part=repmat(log(SteadyState(IdObs)),1,gend);
elseif options_.prefilter == 0 && options_.loglinear == 0 %unlogged steady state must be used
constant_part=repmat(SteadyState(IdObs),1,gend);
end
%add trend to observables
if options_.prefilter
%do correction for prefiltering for observed variables
if options_.loglinear
mean_correction=-repmat(log(SteadyState(IdObs)),1,gend)+constant_part;
else
mean_correction=-repmat(SteadyState(IdObs),1,gend)+constant_part;
end
%smoothed variables are E_T(y_t) so no trend shift is required
stock_smooth(IdObs,:,irun(1))=stock_smooth(IdObs,:,irun(1))+trend_addition+mean_correction;
%updated variables are E_t(y_t) so no trend shift is required
stock_update(IdObs,:,irun(1))=stock_update(IdObs,:,irun(1))+trend_addition+mean_correction;
else
stock_smooth(IdObs,:,irun(1))=stock_smooth(IdObs,:,irun(1))+trend_addition;
stock_update(IdObs,:,irun(1))=stock_update(IdObs,:,irun(1))+trend_addition;
end
stock_innov(:,:,irun(2)) = etahat; stock_innov(:,:,irun(2)) = etahat;
if nvn if nvn
stock_error(:,:,irun(3)) = epsilonhat; stock_error(:,:,irun(3)) = epsilonhat;
end end
if naK if naK
%filtered variable E_t(y_t+k) requires to shift trend by k periods
%write percentage deviation of variables into declaration order
stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:); stock_filter_step_ahead(:,dr.order_var,:,irun(4)) = aK(options_.filter_step_ahead,1:endo_nbr,:);
%now add trend and constant to filtered variables
for ii=1:length(options_.filter_step_ahead)
stock_filter_step_ahead(ii,IdObs,:,irun(4)) = squeeze(stock_filter_step_ahead(ii,IdObs,:,irun(4)))...
+repmat(constant_part(:,1),1,gend+max(options_.filter_step_ahead))... %constant
+[trend_addition repmat(trend_addition(:,end),1,max(options_.filter_step_ahead))+trend_coeff*[1:max(options_.filter_step_ahead)]]; %trend
end
end end
if horizon if horizon
yyyy = alphahat(iendo,i_last_obs); yyyy = alphahat(iendo,i_last_obs);
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr)); yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
if options_.prefilter if options_.prefilter
yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ... % add mean
yf(:,IdObs) = yf(:,IdObs)+repmat(mean_varobs, ...
horizon+maxlag,1); horizon+maxlag,1);
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf(:,IdObs) = yf(:,IdObs)+((options_.first_obs-1)+gend+[1-maxlag:horizon]')*trend_coeff'-...
repmat(mean(trend_coeff*[options_.first_obs:options_.first_obs+gend-1],2)',length(1-maxlag:horizon),1); %center trend
else
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf(:,IdObs) = yf(:,IdObs)+((options_.first_obs-1)+gend+[1-maxlag:horizon]')*trend_coeff';
end end
yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
if options_.loglinear if options_.loglinear
yf = yf+repmat(log(SteadyState'),horizon+maxlag,1); yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
else else
@ -203,11 +243,17 @@ for b=fpar:B
end end
yf1 = forcst2(yyyy,horizon,dr,1); yf1 = forcst2(yyyy,horizon,dr,1);
if options_.prefilter == 1 if options_.prefilter == 1
% add mean
yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ... yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]); repmat(mean_varobs,[horizon+maxlag,1,1]);
end % add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ... yf1(:,IdObs) = yf1(:,IdObs)+((options_.first_obs-1)+gend+[1-maxlag:horizon]')*trend_coeff'-...
repmat(mean(trend_coeff*[options_.first_obs:options_.first_obs+gend-1],2)',length(1-maxlag:horizon),1); %center trend
else
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat(((options_.first_obs-1)+gend+[1-maxlag:horizon]')* ...
trend_coeff',[1,1,1]); trend_coeff',[1,1,1]);
end
if options_.loglinear if options_.loglinear
yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]); yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
else else

4
matlab/random_walk_metropolis_hastings_core.m
View File

@ -109,6 +109,7 @@ proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
block_iter=0; block_iter=0;
for curr_block = fblck:nblck, for curr_block = fblck:nblck,
LastSeeds=[];
block_iter=block_iter+1; block_iter=block_iter+1;
try try
% This will not work if the master uses a random number generator not % This will not work if the master uses a random number generator not
@ -171,7 +172,8 @@ for curr_block = fblck:nblck,
dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]); dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_block) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/draw_iter)]);
end end
if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done if (draw_index_current_file == InitSizeArray(curr_block)) || (draw_iter == nruns(curr_block)) % Now I save the simulations, either because the current file is full or the chain is done
save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2'); [LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal] = get_dynare_random_generator_state();
save([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'],'x2','logpo2','LastSeeds');
fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); fidlog = fopen([MetropolisFolder '/metropolis.log'],'a');
fprintf(fidlog,['\n']); fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']); fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']);

6
matlab/reversed_extended_path.m
View File

@ -14,7 +14,7 @@ function innovation_paths = reversed_extended_path(controlled_variable_names, co
% %
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% Copyright (C) 2010-2011 Dynare Team. % Copyright (C) 2010-2016 Dynare Team.
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
@ -56,10 +56,10 @@ options_.order = old_options_order;
options_.periods = 100; options_.periods = 100;
% Set-up oo_.exo_simul. % Set-up oo_.exo_simul.
make_ex_; oo_=make_ex_(M_,options_,oo_);
% Set-up oo_.endo_simul. % Set-up oo_.endo_simul.
make_y_; oo_=make_y_(M_,options_,oo_);
% Get indices of the controlled endogenous variables in endo_simul. % Get indices of the controlled endogenous variables in endo_simul.
n = length(controlled_variable_names); n = length(controlled_variable_names);

6
matlab/simult_.m
View File

@ -4,7 +4,7 @@ function y_=simult_(y0,dr,ex_,iorder)
% %
% INPUTS % INPUTS
% y0 [double] n*1 vector, initial value (n is the number of declared endogenous variables plus the number % y0 [double] n*1 vector, initial value (n is the number of declared endogenous variables plus the number
% of auxilliary variables for lags and leads) % of auxilliary variables for lags and leads); must be in declaration order, i.e. as in M_.endo_names
% dr [struct] matlab's structure where the reduced form solution of the model is stored. % dr [struct] matlab's structure where the reduced form solution of the model is stored.
% ex_ [double] T*q matrix of innovations. % ex_ [double] T*q matrix of innovations.
% iorder [integer] order of the taylor approximation. % iorder [integer] order of the taylor approximation.
@ -41,6 +41,10 @@ exo_nbr = M_.exo_nbr;
y_ = zeros(size(y0,1),iter+M_.maximum_lag); y_ = zeros(size(y0,1),iter+M_.maximum_lag);
y_(:,1) = y0; y_(:,1) = y0;
if options_.loglinear && ~options_.logged_steady_state
dr.ys=log(dr.ys);
end
if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if k_order_pert is not used or if we do not use Dynare++ with k_order_pert if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if k_order_pert is not used or if we do not use Dynare++ with k_order_pert
if iorder==1 if iorder==1
y_(:,1) = y_(:,1)-dr.ys; y_(:,1) = y_(:,1)-dr.ys;

6
matlab/stoch_simul.m
View File

@ -78,14 +78,16 @@ else
if options_.logged_steady_state %if steady state was previously logged, undo this if options_.logged_steady_state %if steady state was previously logged, undo this
oo_.dr.ys=exp(oo_.dr.ys); oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state); oo_.steady_state=exp(oo_.steady_state);
options_.logged_steady_state=0;
end end
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
end end
if options_.loglinear && isfield(oo_.dr,'ys') %log steady state for correct display of decision rules and simulations if options_.loglinear && isfield(oo_.dr,'ys') && options_.logged_steady_state==0 %log steady state for correct display of decision rule
oo_.dr.ys=log(oo_.dr.ys); oo_.dr.ys=log(oo_.dr.ys);
oo_.steady_state=log(oo_.steady_state); oo_.steady_state=log(oo_.steady_state);
options_old.logged_steady_state = 1; options_old.logged_steady_state = 1; %make sure option is preserved outside of stoch_simul
options_.logged_steady_state=1; %set option for use in stoch_simul
end end
if info(1) if info(1)
options_ = options_old; options_ = options_old;

154
matlab/write_smoother_results.m Normal file
View File

@ -0,0 +1,154 @@
function [oo_, yf]=write_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend)
% oo_=write_smoother_results(M_,oo_,options_,bayestopt_,dataset_,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend)
% Writes the smoother results into respective fields in oo_
%
% Inputs:
% M_ [structure] storing the model information
% oo_ [structure] storing the results
% options_ [structure] storing the options
% bayestopt_ [structure] storing information about priors
% dataset_ [structure] storing the dataset
% atT [double] (m*T) matrix, smoothed endogenous variables (a_{t|T})
% innov [double] (r*T) matrix, smoothed structural shocks (r>n is the umber of shocks).
% measurement_error [double] (n*T) matrix, smoothed measurement errors.
% updated_variables [double] (m*T) matrix, updated (endogenous) variables (a_{t|T})
% ys [double] (m*1) vector specifying the steady state level of each endogenous variable.
% trend_coeff [double] (n*1) vector, parameters specifying the slope of the trend associated to each observed variable.
% aK [double] (K,n,T+K) array, k (k=1,...,K) steps ahead filtered (endogenous) variables.
% P [3D array] of one-step ahead forecast error variance
% matrices
% PK [4D array] of k-step ahead forecast error variance
% matrices (meaningless for periods 1:d)
% decomp
% Trend [double] [nvarobs*T] matrix of trends in observables
%
% Outputs:
% oo_ [structure] storing the results
% yf [double] (nvarobs*T) matrix storing the smoothed observed variables
%
% Notes: first all smoothed variables are saved without trend and constant.
% Then trend and constant are added for the observed variables.
%
% Copyright (C) 2014 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/>.
gend=dataset_.nobs;
if nargin<16
Trend=zeros(options_.number_of_observed_variables,gend);
end
%% write variable steady state
oo_.Smoother.SteadyState = ys;
%% write trend coefficients and trend
oo_.Smoother.TrendCoeffs = trend_coeff; %are in order of options_.varobs
if ~isempty(Trend)
for var_iter=1:size(options_.varobs,2)
oo_.Smoother.Trend.(deblank(options_.varobs{1,var_iter})) = Trend(var_iter,:)';
end
end
%% Compute constant for observables
if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
constant_part=repmat(dataset_info.descriptive.mean',1,gend);
elseif options_.prefilter == 0 && options_.loglinear == 1 %logged steady state must be used
constant_part=repmat(log(ys(bayestopt_.mfys)),1,gend);
elseif options_.prefilter == 0 && options_.loglinear == 0 %unlogged steady state must be used
constant_part=repmat(ys(bayestopt_.mfys),1,gend);
end
%% get observed variables including trend and constant
trend_constant_observables=constant_part+Trend;
yf = atT(bayestopt_.mf,:)+trend_constant_observables;
if options_.nk > 0
%filtered variable E_t(y_t+k) requires to shift trend by k periods
filter_steps_required=union(1,options_.filter_step_ahead); % 1 is required for standard filtered variables
for filter_iter=1:length(filter_steps_required)
filter_step=filter_steps_required(filter_iter);
trend_constant_observables_filtered.(['filter_ahead_' num2str(filter_step)])=constant_part+[Trend+repmat(filter_step*trend_coeff,1,gend)];
end
end
%% write smoother variance
if options_.filter_covariance
oo_.Smoother.Variance = P;
end
%get indicees of smoothed variables
i_endo = bayestopt_.smoother_saved_var_list;
if ~isempty(options_.nk) && options_.nk ~= 0 && (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)))
%write deviations from steady state, add constant for observables later
oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead,i_endo,:);
if ~isempty(PK) %get K-step ahead variances
oo_.FilteredVariablesKStepAheadVariances = ...
PK(options_.filter_step_ahead,i_endo,i_endo,:);
end
if ~isempty(decomp) %get decomposition
oo_.FilteredVariablesShockDecomposition = ...
decomp(options_.filter_step_ahead,i_endo,:,:);
end
end
for i=bayestopt_.smoother_saved_var_list'
i1 = oo_.dr.order_var(bayestopt_.smoother_var_list(i)); %get indices of smoothed variables in name vector
%% Compute constant
if options_.loglinear == 1 %logged steady state must be used
constant_current_variable=repmat(log(ys(i1)),gend,1);
elseif options_.loglinear == 0 %unlogged steady state must be used
constant_current_variable=repmat((ys(i1)),gend,1);
end
oo_.SmoothedVariables.(deblank(M_.endo_names(i1,:)))=atT(i,:)'+constant_current_variable;
if ~isempty(options_.nk) && options_.nk > 0 && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))
oo_.FilteredVariables.(deblank(M_.endo_names(i1,:)))=squeeze(aK(1,i,2:end-(options_.nk-1)));
end
oo_.UpdatedVariables.(deblank(M_.endo_names(i1,:)))=updated_variables(i,:)'+constant_current_variable;
end
%% Add trend and constant for observed variables
for pos_iter=1:length(bayestopt_.mf)
oo_.Smoother.Constant.(deblank(M_.endo_names(bayestopt_.mfys(pos_iter),:)))=constant_part(pos_iter,:);
oo_.SmoothedVariables.(deblank(M_.endo_names(bayestopt_.mfys(pos_iter),:)))=yf(pos_iter,:)';
if ~isempty(options_.nk) && options_.nk > 0 && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))
%filtered variable E_t(y_t+1) requires to shift trend by 1 period
oo_.FilteredVariables.(deblank(M_.endo_names(bayestopt_.mfys(pos_iter),:)))=...
squeeze(aK(1,bayestopt_.mf(pos_iter),2:end-(options_.nk-1)))...
+trend_constant_observables_filtered.filter_ahead_1(pos_iter,:)';
for filter_iter=1:length(options_.filter_step_ahead)
filter_step=options_.filter_step_ahead(filter_iter);
oo_.FilteredVariablesKStepAhead(filter_iter,bayestopt_.mf(pos_iter),1+filter_step:end-(max(options_.filter_step_ahead)-filter_step)) = ...
squeeze(aK(filter_step,bayestopt_.mf(pos_iter),1+filter_step:end-(max(options_.filter_step_ahead)-filter_step)))...
+trend_constant_observables_filtered.(['filter_ahead_' num2str(filter_step)])(pos_iter,:)';
end
end
%updated variables are E_t(y_t) so no trend shift is required
oo_.UpdatedVariables.(deblank(M_.endo_names(bayestopt_.mfys(pos_iter),:)))=...
updated_variables(bayestopt_.mf(pos_iter),:)'+trend_constant_observables(pos_iter,:)';
end
%% get smoothed shocks
for i=1:M_.exo_nbr
oo_.SmoothedShocks.(deblank(M_.exo_names(i,:)))=innov(i,:)';
end
%% Smoothed measurement errors
if ~isequal(M_.H,0)
% measurement_error_indices=find(diag(M_.H)~=0);
for meas_error_iter=1:length(options_.varobs)
oo_.SmoothedMeasurementErrors.(options_.varobs{meas_error_iter})= measurement_error(meas_error_iter,:)';
end
end

9
preprocessor/ComputingTasks.cc
View File

@ -559,13 +559,6 @@ EstimationStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
it = options_list.num_options.find("dsge_var"); it = options_list.num_options.find("dsge_var");
if (it != options_list.num_options.end()) if (it != options_list.num_options.end())
// Ensure that irf_shocks & dsge_var have not both been passed
if (options_list.symbol_list_options.find("irf_shocks") != options_list.symbol_list_options.end())
{
cerr << "The irf_shocks and dsge_var options may not both be passed to estimation." << endl;
exit(EXIT_FAILURE);
}
else
// Fill in mod_file_struct.dsge_var_calibrated // Fill in mod_file_struct.dsge_var_calibrated
mod_file_struct.dsge_var_calibrated = it->second; mod_file_struct.dsge_var_calibrated = it->second;
@ -1021,7 +1014,7 @@ ObservationTrendsStatement::ObservationTrendsStatement(const trend_elements_t &t
void void
ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
{ {
output << "options_.trend_coeff_ = {};" << endl; output << "options_.trend_coeff = {};" << endl;
trend_elements_t::const_iterator it; trend_elements_t::const_iterator it;

1
preprocessor/DynareBison.yy
View File

@ -2481,6 +2481,7 @@ calib_smoother_option : o_filtered_vars
| o_prefilter | o_prefilter
| o_loglinear | o_loglinear
| o_first_obs | o_first_obs
| o_filter_decomposition
; ;
extended_path : EXTENDED_PATH ';' extended_path : EXTENDED_PATH ';'

2
preprocessor/ModFile.cc
View File

@ -798,6 +798,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
<< " save('" << basename << "_results.mat', 'dataset_', '-append');" << endl << "end" << endl << " save('" << basename << "_results.mat', 'dataset_', '-append');" << endl << "end" << endl
<< "if exist('estimation_info', 'var') == 1" << endl << "if exist('estimation_info', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'estimation_info', '-append');" << endl << "end" << endl << " save('" << basename << "_results.mat', 'estimation_info', '-append');" << endl << "end" << endl
<< "if exist('dataset_info', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'dataset_info', '-append');" << endl << "end" << endl
<< "if exist('oo_recursive_', 'var') == 1" << endl << "if exist('oo_recursive_', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'oo_recursive_', '-append');" << endl << "end" << endl; << " save('" << basename << "_results.mat', 'oo_recursive_', '-append');" << endl << "end" << endl;

58
tests/Makefile.am
View File

@ -7,6 +7,9 @@ MODFILES = \
estimation/fs2000_model_comparison.mod \ estimation/fs2000_model_comparison.mod \
estimation/fs2000_fast.mod \ estimation/fs2000_fast.mod \
estimation/MH_recover/fs2000_recover.mod \ estimation/MH_recover/fs2000_recover.mod \
estimation/MH_recover/fs2000_recover_2.mod \
estimation/MH_recover/fs2000_recover_3.mod \
estimation/MH_recover/fs2000_recover_tarb.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_var_decomp.mod \
@ -166,6 +169,14 @@ MODFILES = \
ms-sbvar/test_ms_variances_repeated_runs.mod \ ms-sbvar/test_ms_variances_repeated_runs.mod \
ms-dsge/test_ms_dsge.mod \ ms-dsge/test_ms_dsge.mod \
kalman/lyapunov/fs2000_lyap.mod \ kalman/lyapunov/fs2000_lyap.mod \
kalman/lik_init/fs2000_ns_lik_init_2.mod \
kalman/lik_init/fs2000_ns_lik_init_3.mod \
kalman/lik_init/fs2000_ns_lik_init_5.mod \
kalman/lik_init/fs2000_lik_init_1.mod \
kalman/lik_init/fs2000_lik_init_2.mod \
kalman/lik_init/fs2000_lik_init_3.mod \
kalman/lik_init/fs2000_lik_init_4.mod \
kalman/lik_init/fs2000_lik_init_5.mod \
kalman_filter_smoother/gen_data.mod \ kalman_filter_smoother/gen_data.mod \
kalman_filter_smoother/algo1.mod \ kalman_filter_smoother/algo1.mod \
kalman_filter_smoother/algo2.mod \ kalman_filter_smoother/algo2.mod \
@ -188,6 +199,11 @@ MODFILES = \
kalman/likelihood_from_dynare/fs2000_uncorr_ME.mod \ kalman/likelihood_from_dynare/fs2000_uncorr_ME.mod \
kalman/likelihood_from_dynare/fs2000_uncorr_ME_missing.mod \ kalman/likelihood_from_dynare/fs2000_uncorr_ME_missing.mod \
second_order/burnside_1.mod \ second_order/burnside_1.mod \
kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod \
kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod \
kalman_filter_smoother/compare_results_simulation/fs2000.mod \
kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod \
second_order/burnside_1.mod \
second_order/ds1.mod \ second_order/ds1.mod \
second_order/ds2.mod \ second_order/ds2.mod \
ep/rbc.mod \ ep/rbc.mod \
@ -260,6 +276,31 @@ MODFILES = \
differentiate_forward_vars/RBC_differentiate_forward.mod \ differentiate_forward_vars/RBC_differentiate_forward.mod \
TeX/fs2000_corr_ME.mod \ TeX/fs2000_corr_ME.mod \
prior_posterior_function/fs2000_prior_posterior_function.mod prior_posterior_function/fs2000_prior_posterior_function.mod
observation_trends_and_prefiltering/MCMC/Trend_loglin_no_prefilt_first_obs_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_loglin_prefilt_first_obs_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_loglinear_no_prefilter_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_loglinear_prefilter_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_first_obs_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_prefilter_first_obs_MC.mod \
observation_trends_and_prefiltering/MCMC/Trend_prefilter_MC.mod \
observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter.mod \
observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter_first_obs.mod \
observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter.mod \
observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter_first_obs.mod \
observation_trends_and_prefiltering/ML/Trend_no_prefilter.mod \
observation_trends_and_prefiltering/ML/Trend_no_prefilter_first_obs.mod \
observation_trends_and_prefiltering/ML/Trend_prefilter.mod \
observation_trends_and_prefiltering/ML/Trend_prefilter_first_obs.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_calib_smoother.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilt_first_obs_cal_smooth.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_calib_smoother.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_prefilt_first_obs_cal_smooth.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_loglin_calib_smoother.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_no_prefil_f_obs_loglin_cal_smoother.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod \
observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod \
reporting/example1.mod
XFAIL_MODFILES = ramst_xfail.mod \ XFAIL_MODFILES = ramst_xfail.mod \
estim_param_in_shock_value_xfail.mod \ estim_param_in_shock_value_xfail.mod \
@ -458,6 +499,9 @@ EXTRA_DIST = \
kalman/likelihood_from_dynare/fsdat_simul_corr_ME_missing.m \ kalman/likelihood_from_dynare/fsdat_simul_corr_ME_missing.m \
kalman/likelihood_from_dynare/fsdat_simul_uncorr_ME.m \ kalman/likelihood_from_dynare/fsdat_simul_uncorr_ME.m \
kalman/likelihood_from_dynare/fsdat_simul_uncorr_ME_missing.m \ kalman/likelihood_from_dynare/fsdat_simul_uncorr_ME_missing.m \
kalman/lik_init/fs2000_common.inc \
kalman/lik_init/fs2000_ns_common.inc \
kalman_filter_smoother/compare_results_simulation/fsdat_simul_logged.m \
identification/kim/kim2_steadystate.m \ identification/kim/kim2_steadystate.m \
identification/as2007/as2007_steadystate.m \ identification/as2007/as2007_steadystate.m \
estimation/fsdat_simul.m \ estimation/fsdat_simul.m \
@ -469,8 +513,22 @@ EXTRA_DIST = \
smoother2histval/fsdat_simul.m \ smoother2histval/fsdat_simul.m \
optimal_policy/Ramsey/find_c.m \ optimal_policy/Ramsey/find_c.m \
optimal_policy/Ramsey/oo_ramsey_policy_initval.mat \ optimal_policy/Ramsey/oo_ramsey_policy_initval.mat \
observation_trends_and_prefiltering/Trend_diagnostics_MCMC_common.inc \
observation_trends_and_prefiltering/Trend_diagnostics_calib_common.inc \
observation_trends_and_prefiltering/Trend_diagnostics_ML_common.inc \
observation_trends_and_prefiltering/Trend_exp_model_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_exp_model_no_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_model_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_model_no_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_exp_model_calib_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_exp_model_calib_no_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_model_calib_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_model_calib_no_prefilter_common.inc \
observation_trends_and_prefiltering/Trend_load_data_common.inc
optimal_policy/Ramsey/oo_ramsey_policy_initval.mat \
optimizers/optimizer_function_wrapper.m \ optimizers/optimizer_function_wrapper.m \
optimizers/fs2000.common.inc \ optimizers/fs2000.common.inc \
estimation/MH_recover/fs2000.common.inc \
prior_posterior_function/posterior_function_demo.m prior_posterior_function/posterior_function_demo.m

116
tests/estimation/MH_recover/fs2000.common.inc Normal file
View File

@ -0,0 +1,116 @@
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
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;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady_state_model;
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;
end;
steady;
check;
estimated_params;
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;
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;
end;
varobs gp_obs gy_obs;
options_.plot_priors=0;

146
tests/estimation/MH_recover/fs2000_recover.mod
View File

@ -1,139 +1,35 @@
/* //Test mh_recover function for RW-MH
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/* @#include "fs2000.common.inc"
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA; options_.MaxNumberOfBytes=1000*11*8/2;
varexo e_a e_m; estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=1000, mh_nblocks=2, mh_jscale=0.8);
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
parameters alp bet gam mst rho psi del; copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'],[M_.dname '_mh2_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'])
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;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady_state_model;
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;
end;
steady;
check;
estimated_params;
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;
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;
end;
varobs gp_obs gy_obs;
options_.MaxNumberOfBytes=2000*11*8/4;
estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8);
copyfile([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh1_blck1.mat'],'fs2000_mh1_blck1.mat')
copyfile([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh3_blck2.mat'],'fs2000_mh3_blck2.mat')
delete([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh4_blck2.mat'])
estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover); estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
%check first unaffected chain %check first unaffected chain
temp1=load('fs2000_mh1_blck1.mat'); temp1=load([M_.dname '_mh1_blck1.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh1_blck1.mat']); temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
if max(max(temp1.x2-temp2.x2))>1e-10 if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of unaffected chain are not the same') error('Draws of unaffected chain are not the same')
end end
%check second, affected chain with last unaffected file %check second, affected chain with last unaffected file
temp1=load('fs2000_mh3_blck2.mat'); temp1=load([M_.dname '_mh1_blck1.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh3_blck2.mat']); temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
if max(max(temp1.x2-temp2.x2))>1e-10 if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s unaffected files are not the same') error('Draws of affected chain''s unaffected files are not the same')
end end
%check second, affected chain with affected file
temp1=load([M_.dname '_mh2_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s affected files are not the same')
end

49
tests/estimation/MH_recover/fs2000_recover_2.mod Normal file
View File

@ -0,0 +1,49 @@
//Test mh_recover function for RW-MH when load_mh_file was also used
//Previous chain did not fill last mh_file so mh_recover needs to fill that file
@#include "fs2000.common.inc"
options_.MaxNumberOfBytes=2000*11*8/4;
estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=999, mh_nblocks=2, mh_jscale=0.8);
estimation(order=1,mode_compute=0,mode_file=fs2000_recover_2_mode, datafile='../fsdat_simul',nobs=192, loglinear, load_mh_file,mh_replic=1002, mh_nblocks=2, mh_jscale=0.8);
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'],[M_.dname '_mh3_blck2.mat'])
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'],[M_.dname '_mh4_blck2.mat'])
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh5_blck2.mat'],[M_.dname '_mh5_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh5_blck2.mat'])
estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_2_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
%check first unaffected chain
temp1=load([M_.dname '_mh1_blck1.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of unaffected chain are not the same')
end
%check second, affected chain with last unaffected file
temp1=load([M_.dname '_mh3_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s unaffected files are not the same')
end
%check second, affected chain with affected file
temp1=load([M_.dname '_mh4_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s affected files are not the same')
end
%check second, affected chain with affected file
temp1=load([M_.dname '_mh5_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh5_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s affected files are not the same')
end

39
tests/estimation/MH_recover/fs2000_recover_3.mod Normal file
View File

@ -0,0 +1,39 @@
//Test mh_recover function for RW-MH when load_mh_file was also used
//Previous chain did fill last mh_file so mh_recover does not need to fill that file
@#include "fs2000.common.inc"
options_.MaxNumberOfBytes=2000*11*8/4;
estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=1000, mh_nblocks=2, mh_jscale=0.8);
estimation(order=1,mode_compute=0,mode_file=fs2000_recover_3_mode, datafile='../fsdat_simul',nobs=192, loglinear, load_mh_file,mh_replic=1000, mh_nblocks=2, mh_jscale=0.8);
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'],[M_.dname '_mh3_blck2.mat'])
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'],[M_.dname '_mh4_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat'])
estimation(order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_3_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
%check first unaffected chain
temp1=load([M_.dname '_mh1_blck1.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of unaffected chain are not the same')
end
%check second, affected chain with last unaffected file
temp1=load([M_.dname '_mh3_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh3_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s unaffected files are not the same')
end
%check second, affected chain with affected file
temp1=load([M_.dname '_mh4_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh4_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s affected files are not the same')
end

35
tests/estimation/MH_recover/fs2000_recover_tarb.mod Normal file
View File

@ -0,0 +1,35 @@
//Test mh_recover function for RW-MH when use_tarb was also used
@#include "fs2000.common.inc"
options_.MaxNumberOfBytes=10*11*8/2;
estimation(use_tarb,order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=10, mh_nblocks=2, mh_jscale=0.8);
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat'],[M_.dname '_mh1_blck1.mat'])
copyfile([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'],[M_.dname '_mh2_blck2.mat'])
delete([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat'])
estimation(use_tarb,order=1, datafile='../fsdat_simul',mode_compute=0,mode_file=fs2000_recover_mode, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8,mh_recover);
%check first unaffected chain
temp1=load([M_.dname '_mh1_blck1.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of unaffected chain are not the same')
end
%check second, affected chain with last unaffected file
temp1=load([M_.dname '_mh1_blck1.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh1_blck1.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s unaffected files are not the same')
end
%check second, affected chain with affected file
temp1=load([M_.dname '_mh2_blck2.mat']);
temp2=load([M_.dname filesep 'metropolis' filesep M_.dname '_mh2_blck2.mat']);
if max(max(abs(temp1.x2-temp2.x2)))>1e-10
error('Draws of affected chain''s affected files are not the same')
end

117
tests/kalman/lik_init/fs2000_common.inc Normal file
View File

@ -0,0 +1,117 @@
/*
* This file is based on the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* Copyright (C) 2004-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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del theta;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
theta=0;
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;
steady_state_model;
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;
end;
varobs gp_obs gy_obs;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
corr gy_obs,gp_obs = 0.5;
end;
steady;
estimated_params;
alp, 0.356;
gam, 0.0085;
del, 0.01;
stderr e_a, 0.035449;
stderr e_m, 0.008862;
corr e_m, e_a, 0;
stderr gp_obs, 1;
stderr gy_obs, 1;
corr gp_obs, gy_obs,0;
end;
data(file='../../fs2000/fsdat_simul.m');
[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = ...
dynare_estimation_init(char(), M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_);

24
tests/kalman/lik_init/fs2000_lik_init_1.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_common.inc"
options_.lik_init = 1;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

24
tests/kalman/lik_init/fs2000_lik_init_2.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_common.inc"
options_.lik_init = 2;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

24
tests/kalman/lik_init/fs2000_lik_init_3.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_common.inc"
options_.lik_init = 3;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

24
tests/kalman/lik_init/fs2000_lik_init_4.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_common.inc"
options_.lik_init = 4;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

24
tests/kalman/lik_init/fs2000_lik_init_5.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_common.inc"
options_.lik_init = 5;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

124
tests/kalman/lik_init/fs2000_ns_common.inc Normal file
View File

@ -0,0 +1,124 @@
/*
* This file is based on the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model. This version estimates the model in level rather than in growth rates.
*/
/*
* Copyright (C) 2004-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/>.
*/
var m P c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA;
varexo e_a e_m;
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;
Y_obs/Y_obs(-1) = gy_obs;
P_obs/P_obs(-1) = gp_obs;
end;
steady_state_model;
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;
gy_obs = exp(gam);
gp_obs = exp(-gam);
Y_obs=gy_obs;
P_obs=gp_obs;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
check(nocheck);
estimated_params;
alp, 0.356;
gam, 0.0085;
del, 0.01;
stderr e_a, 0.035449;
stderr e_m, 0.008862;
corr e_m, e_a, 0;
stderr P_obs, 1;
stderr Y_obs, 1;
corr P_obs, Y_obs,0;
end;
varobs P_obs Y_obs;
observation_trends;
P_obs (log(mst)-gam);
Y_obs (gam);
end;
data(file='../../fs2000/fsdat_simul.m');
[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = ...
dynare_estimation_init(char(), M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_);

24
tests/kalman/lik_init/fs2000_ns_lik_init_2.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_ns_common.inc"
options_.lik_init = 2;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

24
tests/kalman/lik_init/fs2000_ns_lik_init_3.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_ns_common.inc"
options_.lik_init = 3;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

24
tests/kalman/lik_init/fs2000_ns_lik_init_5.mod Normal file
View File

@ -0,0 +1,24 @@
/*
* 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/>.
*/
@#include "fs2000_ns_common.inc"
options_.lik_init = 5;
oo_ = initial_estimation_checks(str2func('dsge_likelihood'),xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);

158
tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod Normal file
View File

@ -0,0 +1,158 @@
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
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(+1)*P(+1)*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));
exp(gy_obs) = dA*y/y(-1);
exp(gp_obs) = (P/P(-1))*m(-1)/dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady_state_model;
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 = log(m/dA);
gy_obs = log(dA);
end;
steady;
check;
estimated_params;
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;
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;
end;
varobs gp_obs gy_obs;
estimation(order=1,datafile=fsdat_simul_logged,consider_all_endogenous,nobs=192,mh_replic=2000, mh_nblocks=1,smoother, mh_jscale=0.8);
ex_=[];
for shock_iter=1:M_.exo_nbr
ex_=[ex_ oo_.SmoothedShocks.Mean.(deblank(M_.exo_names(shock_iter,:)))];
end
ex_ = ex_(2:end,:);
% ex_ = zeros(size(ex_));
y0=[];
for endo_iter=1:M_.endo_nbr
y0 = [y0;
oo_.SmoothedVariables.Mean.(deblank(M_.endo_names(endo_iter,:)))(1)];
end;
%make sure decision rules were updated
[oo_.dr,info,M_,options_] = resol(0,M_,options_,oo_);
dr = oo_.dr;
iorder=1;
y_=simult_(y0,dr,ex_,iorder);
fsdat_simul_logged;
%Needs bigger tolerance than ML, because transformation from parameters to steady states is not linear and steady state at mean parameters is not mean of steady states
if mean(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-(gy_obs(1:options_.nobs))))>1e-3 ||...
mean(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.Mean.gy_obs))>1e-3 ||...
mean(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-(gp_obs(1:options_.nobs))))>1e-1 ||...
mean(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.Mean.gp_obs))>1e-2
error('Smoother is wrong')
end
% figure
% plot((gy_obs))
% hold on
% plot(y_(strmatch('gy_obs',M_.endo_names,'exact'),:),'r--')
%
% figure
% plot((gp_obs))
% hold on
% plot(y_(strmatch('gp_obs',M_.endo_names,'exact'),:),'r--')

163
tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod Normal file
View File

@ -0,0 +1,163 @@
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del theta;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
theta=0;
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));
exp(gy_obs) = dA*y/y(-1);
exp(gp_obs) = (P/P(-1))*m(-1)/dA;
end;
steady_state_model;
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 = log(m/dA);
gy_obs = log(dA);
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
varobs gp_obs gy_obs;
steady;
check;
estimated_params;
alp, 0.356;
gam, 0.0085;
mst, 1.0002;
rho, 0.129;
psi, 0.65;
del, 0.02;
stderr e_a, 0.035449;
stderr e_m, 0.008862;
end;
estimation(order=1,datafile='fsdat_simul_logged', nobs=192, forecast=8,smoother,filtered_vars,filter_step_ahead=[1,2,4],filter_decomposition,selected_variables_only) m P c e W R k d y gy_obs;
% write shock matrix
ex_=[];
for shock_iter=1:M_.exo_nbr
ex_=[ex_ oo_.SmoothedShocks.(deblank(M_.exo_names(shock_iter,:)))];
end
%select shocks happening after initial period
ex_ = ex_(2:end,:);
%get state variables at t=0
y0=[];
for endo_iter=1:M_.endo_nbr
y0 = [y0;
oo_.SmoothedVariables.(deblank(M_.endo_names(endo_iter,:)))(1)];
end;
%make sure decision rules were updated
[oo_.dr,info,M_,options_] = resol(0,M_,options_,oo_);
dr = oo_.dr;
iorder=1;
%run simulation
y_=simult_(y0,dr,ex_,iorder);
fsdat_simul_logged;
if max(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-gy_obs(1:options_.nobs)))>1e-10 ||...
max(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.gy_obs))>1e-10 ||...
max(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-gp_obs(1:options_.nobs)))>1e-10 ||...
max(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.gp_obs))>1e-10
error('Smoother is wrong')
end
% figure
% subplot(2,1,1)
% plot(log(gy_obs))
% hold on
% plot(y_(strmatch('gy_obs',M_.endo_names,'exact'),:),'r--')
%
% figure
% subplot(2,1,2)
% plot(log(gp_obs))
% hold on
% plot(y_(strmatch('gp_obs',M_.endo_names,'exact'),:),'r--')

150
tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod Normal file
View File

@ -0,0 +1,150 @@
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
parameters alp bet gam mst rho psi del theta;
alp = 0.33;
bet = 0.99;
gam = 0.003;
mst = 1.011;
rho = 0.7;
psi = 0.787;
del = 0.02;
theta=0;
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;
steady_state_model;
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;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
varobs gp_obs gy_obs;
steady;
check;
estimated_params;
alp, 0.356;
gam, 0.0085;
mst, 1.0002;
rho, 0.129;
psi, 0.65;
del, 0.02;
stderr e_a, 0.035449;
stderr e_m, 0.008862;
end;
estimation(order=1,datafile='../fsdat_simul',loglinear, nobs=192, forecast=8,smoother,filtered_vars,filter_step_ahead=[1,2,4],filter_decomposition,selected_variables_only) m P c e W R k d y gy_obs;
% write shock matrix
ex_=[];
for shock_iter=1:M_.exo_nbr
ex_=[ex_ oo_.SmoothedShocks.(deblank(M_.exo_names(shock_iter,:)))];
end
%select shocks happening after initial period
ex_ = ex_(2:end,:);
%get state variables at t=0
y0=[];
for endo_iter=1:M_.endo_nbr
y0 = [y0;
oo_.SmoothedVariables.(deblank(M_.endo_names(endo_iter,:)))(1)];
end;
%make sure decision rules were updated
[oo_.dr,info,M_,options_] = resol(0,M_,options_,oo_);
dr = oo_.dr;
iorder=1;
%run simulation
y_=simult_(y0,dr,ex_,iorder);
fsdat_simul_logged;
if max(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-gy_obs(1:options_.nobs)))>1e-10 ||...
max(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.gy_obs))>1e-10 ||...
max(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-gp_obs(1:options_.nobs)))>1e-10 ||...
max(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.gp_obs))>1e-10
error('Smoother is wrong')
end

176
tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod Normal file
View File

@ -0,0 +1,176 @@
/*
* This file replicates the estimation of the cash in advance model described
* Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
* Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
*
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
* implications of long-run neutrality for monetary business cycle models",
* Journal of Applied Econometrics, 9, S37-S70.
* Note that there is an initial minus sign missing in equation (A1), p. S63.
*
* This implementation was written by Michel Juillard. Please note that the
* following copyright notice only applies to this Dynare implementation of the
* model.
*/
/*
* 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/>.
*/
var m P c e W R k d n l gy_obs gp_obs y dA;
varexo e_a e_m;
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(+1)*P(+1)*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;
end;
steady_state_model;
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;
end;
steady;
check;
estimated_params;
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;
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;
end;
varobs gp_obs gy_obs;
estimation(order=1, datafile='../fsdat_simul', nobs=192, loglinear, mh_replic=2000, mh_nblocks=1,smoother, mh_jscale=0.8,consider_all_endogenous);
ex_=[];
for shock_iter=1:M_.exo_nbr
ex_=[ex_ oo_.SmoothedShocks.Mean.(deblank(M_.exo_names(shock_iter,:)))];
end
ex_ = ex_(2:end,:);
% ex_ = zeros(size(ex_));
y0=[];
for endo_iter=1:M_.endo_nbr
y0 = [y0;
oo_.SmoothedVariables.Mean.(deblank(M_.endo_names(endo_iter,:)))(1)];
end;
%make sure decision rules were updated
[oo_.dr,info,M_,options_] = resol(0,M_,options_,oo_);
dr = oo_.dr;
iorder=1;
% if options_.loglinear
% y0=exp(y0);
% end
y_=simult_(y0,dr,ex_,iorder);
fsdat_simul_logged;
if mean(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-gy_obs(1:options_.nobs)))>1e-3 ||...
mean(abs(y_(strmatch('gy_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.Mean.gy_obs))>1e-3 ||...
mean(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-gp_obs(1:options_.nobs)))>1e-3 ||...
mean(abs(y_(strmatch('gp_obs',M_.endo_names,'exact'),:)'-oo_.SmoothedVariables.Mean.gp_obs))>1e-3
error('Smoother is wrong')
end
% figure
% plot(log(gy_obs))
% hold on
% plot(y_(strmatch('gy_obs',M_.endo_names,'exact'),:),'r--')
%
% figure
% plot(log(gp_obs))
% hold on
% plot(y_(strmatch('gp_obs',M_.endo_names,'exact'),:),'r--')

832
tests/kalman_filter_smoother/compare_results_simulation/fsdat_simul_logged.m Normal file
View File

@ -0,0 +1,832 @@
gy_obs =[
1.0030045
0.99990934
1.0172778
0.99464043
1.0253423
1.0150215
0.97772557
0.97832186
1.0159561
1.0085937
1.0102649
1.0007604
1.0112596
1.0163279
1.0173204
1.0103896
1.0006493
0.99447124
1.0196405
1.0089304
0.99650737
1.0139707
0.97865842
1.0192225
0.99139628
1.0141362
1.0196612
0.97483476
0.99686151
0.99594464
1.0000642
1.0172243
1.0025773
0.97199728
1.0217815
1.0219949
0.99490252
1.0190728
1.0111337
1.0003792
0.98969164
1.010438
1.0216309
1.0016671
1.0357588
0.98803787
1.0093457
1.0177035
0.98548204
1.0274294
1.0141377
1.0091174
0.96427632
1.0083272
1.0007882
0.99038262
1.0031336
0.99500213
0.98203716
0.9889452
1.011632
0.99451949
0.97291047
0.98750871
0.99992418
0.97657318
0.99930448
1.0008515
1.0044064
0.98133792
1.0091702
1.0087023
1.0119876
1.0143019
1.0311061
0.99340471
1.0057428
0.99197259
1.0071019
0.99448853
1.0061819
1.0070088
0.9950913
1.0302318
0.9817693
1.0072885
0.97355282
0.98782586
1.0136674
0.99863956
1.0205668
0.99611384
1.0073805
0.99691529
1.0089194
1.0030467
1.0112006
1.0260523
0.97803331
0.99423374
1.0043727
1.0140173
1.0111473
0.99524348
0.99775943
0.9958619
0.9982344
1.0210212
1.0022288
1.0014801
1.011456
1.0124871
0.99843599
0.99324886
0.99912838
1.003327
1.0072071
1.0115223
1.009266
1.0070554
1.0129916
1.0053413
1.0051638
0.99212952
1.0214422
0.98716707
0.99905788
0.98877357
0.98568476
0.99767393
1.0061791
0.98423439
0.99492949
0.98786999
0.99754239
1.0168619
0.99472384
1.0041658
0.98123181
1.0112882
0.99245422
1.0010255
1.0017799
1.0089968
1.0072824
0.99768475
1.0044726
1.0118678
1.0056385
1.0276965
1.0025122
1.0065161
1.0234338
0.99760167
0.98922272
1.0101918
1.011615
1.0085286
1.0074455
0.98866757
0.99959012
1.0129881
0.99127881
0.97971901
1.0185314
1.020054
1.0132605
0.98063643
0.99490253
1.0101531
1.0004526
1.0059109
0.98974491
1.0062391
1.0216488
0.99398446
0.97786609
1.0019274
0.99587153
1.0095881
1.0111887
0.99457649
0.97896734
1.000172
1.0142951
1.0034224
1.0037242
1.0016059
1.016556
0.99687023
1.0117844
1.0059212
0.98083159
0.98638851
1.0128713
1.0096232
1.0115891
1.0011213
1.0147105
1.0066344
1.0164429
0.99825038
0.99403411
];
gp_obs =[
1.0079715
1.0074573
1.0153107
1.0152677
1.0011653
0.99950061
1.0328311
1.0192317
1.009827
0.99588916
1.007474
1.0113061
0.98696624
0.99978663
0.98240542
0.98861723
0.99008763
1.0185076
1.0052452
0.99447194
1.0092685
1.01208
1.0105237
0.98513875
1.0165628
0.99485934
1.0050255
1.0140756
1.0093128
1.0155868
1.0107023
0.99212762
1.0095465
1.0028435
1.0069437
1.0070473
1.0145902
1.0186922
1.0059917
1.0113072
1.0107386
0.99769196
0.99793444
1.0050791
0.98307821
1.0107594
0.99689982
0.98667064
0.9991662
0.98274722
0.98422032
0.99393016
1.0118567
0.99912781
1.0023744
1.0086662
1.0164773
1.0169327
1.0372478
1.0314242
1.0004256
1.0110541
1.0076575
1.0119851
1.0055188
1.0213959
1.0234416
1.0264917
1.0292725
1.0385184
1.0200999
1.0107697
1.008583
1.0200332
1.0030413
1.0108659
1.0185145
1.0168619
1.0180462
1.0239657
1.0205509
1.0189973
1.0246446
1.0135089
1.0352973
1.0099289
1.0266474
1.0279829
1.0101653
1.041216
1.0103861
1.0114727
1.0054605
1.0190722
1.0114837
1.0179213
1.006082
1.0049696
1.0143629
0.9971036
1.0005602
1.0078403
1.0240222
1.0195063
1.0355136
1.0218743
1.0171331
1.0049817
1.0140974
1.0168431
1.0049966
1.0045568
1.0156414
1.0273055
1.0197653
1.0030624
1.0154993
0.99782084
0.99711648
1.014408
1.0057417
0.99936837
1.0096934
1.0095138
1.0057734
1.0114497
1.0059784
1.0328889
1.0098032
1.0041114
1.0101247
1.0181588
1.0115712
1.0227509
1.0065104
1.0110902
1.0298169
1.0089532
1.0368733
1.0123033
1.0060763
1.0150937
1.0239325
0.99555536
0.99861271
1.0076201
0.99941535
1.0119522
1.0129183
0.99288924
1.0260784
1.0144982
1.0121985
1.0234916
1.02215
1.0190118
1.0172679
1.0118398
1.0002123
1.0092124
1.0071943
0.99508468
1.0019303
1.0030733
0.9964198
1.0027298
0.99797614
1.006942
0.99793928
1.0083214
1.0283732
1.0111102
1.016936
1.0229061
0.98846454
1.0015387
1.0201769
1.0079822
1.0064007
1.0095543
1.0092207
1.0135485
1.0198974
1.0140252
1.0128686
1.0092903
1.0141974
1.0023492
0.99731455
1.0026598
0.99303643
1.0036469
1.0160975
1.0368378
1.0139625
1.01493
1.0113531
1.0114548
0.99833441
0.99648401
0.97645361
1.0154053
1.01703
];
Y_obs =[
1
0.99690484
1.0111781
1.0028141
1.0251518
1.0371688
1.0118899
0.98720726
1.0001589
1.0057481
1.0130085
1.0107643
1.0190194
1.0323428
1.0466587
1.0540438
1.0516886
1.0431553
1.0597913
1.0657172
1.0592201
1.0701863
1.0458402
1.0620582
1.0504499
1.0615817
1.0782384
1.0500687
1.0439257
1.0368658
1.0339255
1.0481453
1.0477181
1.0167109
1.0354878
1.0544782
1.0463762
1.0624445
1.0705737
1.0679484
1.0546356
1.0620691
1.0806955
1.0793581
1.1121124
1.0971458
1.1034869
1.1181859
1.1006634
1.1250883
1.1362214
1.1423343
1.1036061
1.1089288
1.1067125
1.0940906
1.0942197
1.0862174
1.06525
1.0511907
1.0598182
1.0513331
1.0212391
1.0057433
1.002663
0.97623167
0.97253165
0.97037865
0.97178055
0.95011397
0.95627969
0.96197747
0.97096053
0.98225794
1.0103595
1.0007597
1.003498
0.99246608
0.99656347
0.98804749
0.99122491
0.99522926
0.98731605
1.0145434
0.99330816
0.99759216
0.96814048
0.95296183
0.96362471
0.95925977
0.97682205
0.96993138
0.9743074
0.96821818
0.97413308
0.9741753
0.98237142
1.0054193
0.98044807
0.9716773
0.9730455
0.98405828
0.99220103
0.98444001
0.97919493
0.97205233
0.96728223
0.98529893
0.98452324
0.98299888
0.99145042
1.000933
0.99636447
0.98660883
0.98273271
0.98305518
0.98725774
0.99577549
1.002037
1.0060879
1.016075
1.0184118
1.0205711
1.0096961
1.0281337
1.0122963
1.0083497
0.99411874
0.976799
0.97146842
0.97464304
0.95587292
0.94779791
0.93266339
0.92720128
0.94105864
0.93277798
0.93393927
0.91216657
0.92045028
0.9099
0.90792098
0.90669634
0.91268867
0.91696661
0.91164685
0.91311495
0.92197825
0.92461222
0.94930422
0.9488119
0.95232353
0.97275278
0.96734995
0.95356817
0.96075548
0.96936594
0.97489002
0.97933106
0.96499412
0.96157973
0.97156334
0.95983765
0.93655215
0.95207909
0.96912862
0.97938462
0.95701655
0.94891457
0.95606317
0.95351125
0.95641767
0.94315807
0.94639265
0.96503697
0.95601693
0.93087851
0.92980141
0.92266844
0.92925206
0.93743628
0.92900826
0.9049711
0.90213859
0.91342916
0.91384707
0.91456681
0.91316822
0.92671976
0.92058549
0.92936541
0.93228212
0.91010921
0.89349322
0.90336005
0.90997873
0.91856328
0.91668007
0.92838606
0.932016
0.94545438
0.94070026
0.93172987
];
P_obs =[
1
0.99948573
1.0068249
1.0141211
1.0073149
0.99884398
1.0237035
1.0349636
1.036819
1.0247366
1.0242391
1.0275737
1.0065684
0.99838346
0.97281734
0.95346302
0.9355791
0.9461152
0.94338882
0.92988921
0.9311862
0.93529467
0.93784681
0.91501401
0.92360522
0.91049302
0.90754698
0.91365103
0.91499228
0.92260749
0.92533824
0.90949431
0.91106924
0.90594116
0.90491334
0.9039891
0.91060772
0.92132842
0.91934854
0.92268418
0.92545127
0.91517169
0.90513459
0.90224212
0.87734878
0.88013667
0.86906494
0.84776403
0.83895869
0.81373437
0.78998314
0.77594176
0.77982695
0.77098321
0.76538611
0.76608075
0.77458654
0.78354767
0.81282389
0.83627649
0.82873051
0.83181309
0.83149903
0.83551261
0.83305985
0.84648418
0.86195421
0.88047436
0.90177533
0.93232215
0.94445051
0.9472487
0.94786015
0.95992178
0.95499149
0.95788581
0.9684288
0.97731917
0.98739379
1.0033879
1.0159673
1.0269931
1.0436661
1.0492034
1.0765292
1.0784865
1.0971624
1.1171737
1.1193675
1.1526119
1.1550265
1.1585277
1.1560166
1.1671172
1.1706294
1.1805791
1.1786896
1.1756876
1.1820789
1.171211
1.1637997
1.1636684
1.179719
1.1912538
1.2187959
1.2326986
1.2418602
1.2388704
1.2449963
1.2538678
1.2508929
1.2474781
1.255148
1.274482
1.2862757
1.2813665
1.2888943
1.2787436
1.2678886
1.274325
1.2720952
1.263492
1.2652139
1.2667561
1.264558
1.2680362
1.2660431
1.2909605
1.2927921
1.288932
1.2910852
1.3012725
1.3048721
1.3196515
1.3181903
1.321309
1.3431543
1.344136
1.3730377
1.3773695
1.3754742
1.3825964
1.3985574
1.3861412
1.3767823
1.3764309
1.3678747
1.3718554
1.3768022
1.3617199
1.3798267
1.3863533
1.3905803
1.4061004
1.4202788
1.4313191
1.4406155
1.4444837
1.4367244
1.4379653
1.4371881
1.4243012
1.41826
1.4133617
1.40181
1.3965683
1.3865729
1.3855433
1.3755111
1.3758609
1.3962625
1.3994012
1.4083656
1.4233002
1.4037932
1.3973604
1.4095657
1.4095764
1.4080055
1.4095882
1.4108374
1.4164143
1.4283402
1.4343939
1.4392909
1.4406097
1.4468355
1.4412132
1.4305562
1.4252445
1.4103094
1.4059847
1.4141106
1.4429769
1.4489679
1.4559263
1.4593079
1.4627911
1.453154
1.4416665
1.4101485
1.4175823
1.4266407
];
gp_obs=log(gp_obs);
gy_obs=log(gy_obs);
P_obs=log(P_obs);
Y_obs=log(Y_obs);

19
tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_no_prefilt_first_obs_MC.mod Normal file
View File

@ -0,0 +1,19 @@
@#include "../Trend_exp_model_no_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=2000,
mode_compute=4,first_obs=1000,loglinear,smoother,forecast=100,prefilter=0,
mcmc_jumping_covariance='MCMC_jump_covar',
filtered_vars, filter_step_ahead = [1,2,4],
mh_nblocks=1,mh_jscale=0.3) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

21
tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_prefilt_first_obs_MC.mod Normal file
View File

@ -0,0 +1,21 @@
@#include "../Trend_exp_model_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=2000,
mode_compute=4,first_obs=1000,loglinear,smoother,forecast=100,prefilter=1,
mcmc_jumping_covariance='MCMC_jump_covar_prefilter',
filtered_vars, filter_step_ahead = [1,2,4],
mh_nblocks=1,mh_jscale=1e-4) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params_prefilter');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
loaded_par_full=load('../orig_params');
y_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_y',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_y',loaded_par_full.param_names,'exact'));
p_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_p',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_p',loaded_par_full.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

18
tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_no_prefilter_MC.mod Normal file
View File

@ -0,0 +1,18 @@
@#include "../Trend_exp_model_no_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=2000,
mode_compute=4,first_obs=1,loglinear,diffuse_filter,smoother,forecast=100,prefilter=0,
mcmc_jumping_covariance='MCMC_jump_covar',
filtered_vars, filter_step_ahead = [1,2,4],
mh_nblocks=1,mh_jscale=0.3) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

21
tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_prefilter_MC.mod Normal file
View File

@ -0,0 +1,21 @@
@#include "../Trend_exp_model_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=2000,
mode_compute=4,first_obs=1,loglinear,smoother,forecast=100,prefilter=1,
mcmc_jumping_covariance='MCMC_jump_covar_prefilter',
filtered_vars, filter_step_ahead = [1,2,4],
mh_nblocks=1,mh_jscale=1e-4) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params_prefilter');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
loaded_par_full=load('../orig_params');
y_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_y',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_y',loaded_par_full.param_names,'exact'));
p_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_p',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_p',loaded_par_full.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

18
tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_MC.mod Normal file
View File

@ -0,0 +1,18 @@
@#include "../Trend_model_no_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',mh_replic=2000,
mode_compute=4,first_obs=1,smoother,mh_nblocks=1,mh_jscale=0.3,
filtered_vars, filter_step_ahead = [1,2,4],
mcmc_jumping_covariance='MCMC_jump_covar',forecast=100,prefilter=0) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

19
tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_first_obs_MC.mod Normal file
View File

@ -0,0 +1,19 @@
@#include "../Trend_model_no_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',
mh_replic=2000,mode_compute=4,first_obs=1000,smoother,forecast=100,prefilter=0,
mcmc_jumping_covariance='MCMC_jump_covar',
filtered_vars, filter_step_ahead = [1,2,4],
mh_nblocks=1,mh_jscale=0.3) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

21
tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_MC.mod Normal file
View File

@ -0,0 +1,21 @@
@#include "../Trend_model_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',mh_replic=2000,mode_compute=4,
first_obs=1,smoother,prefilter=1,
mh_nblocks=1,mh_jscale=1e-4,
filtered_vars, filter_step_ahead = [1,2,4],
mcmc_jumping_covariance='MCMC_jump_covar_prefilter',forecast=100) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params_prefilter');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
loaded_par_full=load('../orig_params');
y_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_y',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_y',loaded_par_full.param_names,'exact'));
p_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_p',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_p',loaded_par_full.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

21
tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_first_obs_MC.mod Normal file
View File

@ -0,0 +1,21 @@
@#include "../Trend_model_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',mh_replic=2000,mode_compute=4,
first_obs=1000,smoother,prefilter=1,
mh_nblocks=1,mh_jscale=1e-4,
filtered_vars, filter_step_ahead = [1,2,4],
mcmc_jumping_covariance='MCMC_jump_covar_prefilter',forecast=100) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params_prefilter');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
loaded_par_full=load('../orig_params');
y_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_y',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_y',loaded_par_full.param_names,'exact'));
p_forecast_100_periods=loaded_par_full.orig_params(strmatch('const_p',loaded_par_full.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par_full.orig_params(strmatch('g_p',loaded_par_full.param_names,'exact'));
@#include "../Trend_diagnostics_MCMC_common.inc"

19
tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter.mod Normal file
View File

@ -0,0 +1,19 @@
@#include "../Trend_exp_model_no_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=0,
mode_compute=4,first_obs=1,
filtered_vars, filter_step_ahead = [1,2,4],
loglinear,smoother,forecast=100,prefilter=0) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',M_.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',M_.param_names,'exact'));
@#include "../Trend_diagnostics_ML_common.inc"

19
tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter_first_obs.mod Normal file
View File

@ -0,0 +1,19 @@
@#include "../Trend_exp_model_no_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=0,
mode_compute=4,first_obs=1000,
filtered_vars, filter_step_ahead = [1,2,4],
loglinear,smoother,forecast=100,prefilter=0) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',M_.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',M_.param_names,'exact'));
@#include "../Trend_diagnostics_ML_common.inc"

20
tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter.mod Normal file
View File

@ -0,0 +1,20 @@
@#include "../Trend_exp_model_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=0,mode_compute=4,
first_obs=1,smoother,loglinear,
filtered_vars, filter_step_ahead = [1,2,4],
forecast=100,prefilter=1) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'));
@#include "../Trend_diagnostics_ML_common.inc"

20
tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter_first_obs.mod Normal file
View File

@ -0,0 +1,20 @@
@#include "../Trend_exp_model_prefilter_common.inc"
estimation(order=1,datafile='../Exp_AR1_trend_data_with_constant',mh_replic=0,mode_compute=4,
first_obs=1000,smoother,loglinear,
filtered_vars, filter_step_ahead = [1,2,4],
forecast=100,prefilter=1) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'));
@#include "../Trend_diagnostics_ML_common.inc"

18
tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter.mod Normal file
View File

@ -0,0 +1,18 @@
@#include "../Trend_model_no_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',mh_replic=0,
mode_compute=4,first_obs=1,
filtered_vars, filter_step_ahead = [1,2,4],
diffuse_filter,smoother,forecast=100,prefilter=0) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',M_.param_names,'exact'))
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',M_.param_names,'exact'))
@#include "../Trend_diagnostics_ML_common.inc"

19
tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter_first_obs.mod Normal file
View File

@ -0,0 +1,19 @@
@#include "../Trend_model_no_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',mh_replic=0,
mode_compute=4,first_obs=1000,
filtered_vars, filter_step_ahead = [1,2,4],
smoother,forecast=100,prefilter=0) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs(M_.params-loaded_par.orig_params)./loaded_par.orig_params)>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',M_.param_names,'exact'));
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',M_.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',M_.param_names,'exact'));
@#include "../Trend_diagnostics_ML_common.inc"

20
tests/observation_trends_and_prefiltering/ML/Trend_prefilter.mod Normal file
View File

@ -0,0 +1,20 @@
@#include "../Trend_model_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',mh_replic=0,mode_compute=4,
first_obs=1,
filtered_vars, filter_step_ahead = [1,2,4],
smoother,forecast=100,prefilter=1) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'))
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'))
@#include "../Trend_diagnostics_ML_common.inc"

21
tests/observation_trends_and_prefiltering/ML/Trend_prefilter_first_obs.mod Normal file
View File

@ -0,0 +1,21 @@
@#include "../Trend_model_prefilter_common.inc"
estimation(order=1,datafile='../AR1_trend_data_with_constant',
mh_replic=0,mode_compute=4,
filtered_vars, filter_step_ahead = [1,2,4],
first_obs=1000,diffuse_filter,smoother,forecast=100,prefilter=1) P_obs Y_obs junk2;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameter estimates do not match')
end
y_forecast_100_periods=loaded_par.orig_params(strmatch('const_y',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_y',loaded_par.param_names,'exact'))
p_forecast_100_periods=loaded_par.orig_params(strmatch('const_p',loaded_par.param_names,'exact'))+(options_.first_obs+options_.nobs-1+options_.forecast)*loaded_par.orig_params(strmatch('g_p',loaded_par.param_names,'exact'))
@#include "../Trend_diagnostics_ML_common.inc"

50
tests/observation_trends_and_prefiltering/Trend_diagnostics_MCMC_common.inc Normal file
View File

@ -0,0 +1,50 @@
if max(abs(oo_.SmoothedVariables.Mean.Y_obs-Y_obs'))>1e-5 ||...
max(abs(oo_.SmoothedVariables.Mean.P_obs-P_obs'))>1e-5 || ...
max(abs(oo_.SmoothedVariables.Mean.junk2-junk2'))>1e-5
error('Smoothed Variables are wrong')
end
if max(abs(oo_.UpdatedVariables.Mean.Y_obs-Y_obs'))>1e-5 ||...
max(abs(oo_.UpdatedVariables.Mean.P_obs-P_obs'))>1e-5 || ...
max(abs(oo_.UpdatedVariables.Mean.junk2-junk2'))>1e-5
error('Updated Variables are wrong')
end
if mean(abs(oo_.FilteredVariables.Mean.Y_obs(1:end-1)-Y_obs(2:end)'))>1e-3 ||...
mean(abs(oo_.FilteredVariables.Mean.P_obs(1:end-1)-P_obs(2:end)'))>1e-3
error('Filtered Variables are wrong')
end
if abs(corr(oo_.FilteredVariables.Mean.Y_obs(2:end-1)-Y_obs(3:end)',oo_.FilteredVariables.Mean.Y_obs(1:end-2)-Y_obs(2:end-1)'))>2e-2 ||...
abs(corr(oo_.FilteredVariables.Mean.P_obs(2:end-1)-P_obs(3:end)',oo_.FilteredVariables.Mean.P_obs(1:end-2)-P_obs(2:end-1)'))>2e-2 ||...
abs(corr(oo_.FilteredVariables.Mean.junk2(2:end-1)-junk2(3:end)',oo_.FilteredVariables.Mean.junk2(1:end-2)-junk2(2:end-1)'))>2e-2
error('Filtered Variables are wrong')
end
if max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,1,2:end-(options_.nk-1)))-oo_.FilteredVariables.Mean.P_obs))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,2,2:end-(options_.nk-1)))-oo_.FilteredVariables.Mean.Y_obs))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,3,2:end-(options_.nk-1)))-oo_.FilteredVariables.Mean.junk2))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(2,1,3:end-options_.nk))-oo_.FilteredVariables.Mean.P_obs(3:end)))>1e-2 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(2,2,3:end-options_.nk))-oo_.FilteredVariables.Mean.Y_obs(3:end)))>1e-2 ||...
mean(squeeze(oo_.FilteredVariablesKStepAhead(2,3,3:end-options_.nk)))>1e-1
error('FilteredVariablesKStepAhead is wrong')
end
if abs(oo_.PointForecast.Mean.Y_obs(end)- y_forecast_100_periods)>2e-4 || abs(oo_.PointForecast.Mean.P_obs(end)- p_forecast_100_periods)>2e-4
error('Mean Point Forecasts do not match')
end
if abs(oo_.PointForecast.Median.Y_obs(end)- y_forecast_100_periods)>2e-4 || abs(oo_.PointForecast.Median.P_obs(end)- p_forecast_100_periods)>2e-4
error('Median Point Forecasts do not match')
end
if abs(oo_.MeanForecast.Mean.Y_obs(end)- y_forecast_100_periods)>2e-4 || abs(oo_.MeanForecast.Mean.P_obs(end)- p_forecast_100_periods)>2e-4
error('Mean Mean Forecasts do not match')
end
if abs(oo_.MeanForecast.Median.Y_obs(end)- y_forecast_100_periods)>2e-4 || abs(oo_.MeanForecast.Median.P_obs(end)- p_forecast_100_periods)>1e-3
error('Median Mean Forecasts do not match')
end
if abs(mean(oo_.SmoothedShocks.Mean.e_y))>1e-2 || abs(mean(oo_.SmoothedShocks.Mean.e_p))>1e-2 || abs(mean(oo_.SmoothedShocks.Median.e_y))>1e-2 || abs(mean(oo_.SmoothedShocks.Median.e_p))>1e-2
error('Residuals are not mean 0')
end

40
tests/observation_trends_and_prefiltering/Trend_diagnostics_ML_common.inc Normal file
View File

@ -0,0 +1,40 @@
if abs(oo_.forecast.Mean.Y_obs(end)- y_forecast_100_periods)>1e-4 || abs(oo_.forecast.Mean.P_obs(end)- p_forecast_100_periods)>1e-4
error('Forecasts do not match')
end
if max(abs(oo_.SmoothedVariables.Y_obs-Y_obs'))>1e-5 ||...
max(abs(oo_.SmoothedVariables.P_obs-P_obs'))>1e-5 || ...
max(abs(oo_.SmoothedVariables.junk2-junk2'))>1e-5
error('Smoothed Variables are wrong')
end
if max(abs(oo_.UpdatedVariables.Y_obs-Y_obs'))>1e-5 ||...
max(abs(oo_.UpdatedVariables.P_obs-P_obs'))>1e-5 || ...
max(abs(oo_.UpdatedVariables.junk2-junk2'))>1e-5
error('Updated Variables are wrong')
end
if mean(abs(oo_.FilteredVariables.Y_obs(1:end-1)-Y_obs(2:end)'))>1e-3 ||...
mean(abs(oo_.FilteredVariables.P_obs(1:end-1)-P_obs(2:end)'))>1e-3
error('Smoothed Variables are wrong')
end
if abs(corr(oo_.FilteredVariables.Y_obs(2:end-1)-Y_obs(3:end)',oo_.FilteredVariables.Y_obs(1:end-2)-Y_obs(2:end-1)'))>2e-2 ||...
abs(corr(oo_.FilteredVariables.P_obs(2:end-1)-P_obs(3:end)',oo_.FilteredVariables.P_obs(1:end-2)-P_obs(2:end-1)'))>2e-1 ||...
abs(corr(oo_.FilteredVariables.junk2(2:end-1)-junk2(3:end)',oo_.FilteredVariables.junk2(1:end-2)-junk2(2:end-1)'))>2e-2
error('Filtered Variables are wrong')
end
if max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,1,2:end-(options_.nk-1)))-oo_.FilteredVariables.Y_obs))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,2,2:end-(options_.nk-1)))-oo_.FilteredVariables.P_obs))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,3,2:end-(options_.nk-1)))-oo_.FilteredVariables.junk2))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(2,1,3:end-options_.nk))-oo_.FilteredVariables.Y_obs(3:end)))>1e-2 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(2,2,3:end-options_.nk))-oo_.FilteredVariables.P_obs(3:end)))>1e-2 ||...
mean(squeeze(oo_.FilteredVariablesKStepAhead(2,3,3:end-options_.nk)))>1e-1 ||...
mean(squeeze(oo_.FilteredVariablesKStepAhead(2,4,3:end-options_.nk)))>1e-1
error('FilteredVariablesKStepAhead is wrong')
end
if abs(mean(oo_.SmoothedShocks.e_y))>1e-1 || abs(mean(oo_.SmoothedShocks.e_p))>1e-1
error('Residuals are not mean 0')
end

36
tests/observation_trends_and_prefiltering/Trend_diagnostics_calib_common.inc Normal file
View File

@ -0,0 +1,36 @@
if max(abs(oo_.SmoothedVariables.Y_obs-Y_obs'))>1e-5 ||...
max(abs(oo_.SmoothedVariables.P_obs-P_obs'))>1e-5 || ...
max(abs(oo_.SmoothedVariables.junk2-junk2'))>1e-5
error('Smoothed Variables are wrong')
end
if max(abs(oo_.UpdatedVariables.Y_obs-Y_obs'))>1e-5 ||...
max(abs(oo_.UpdatedVariables.P_obs-P_obs'))>1e-5 || ...
max(abs(oo_.UpdatedVariables.junk2-junk2'))>1e-5
error('Updated Variables are wrong')
end
if mean(abs(oo_.FilteredVariables.Y_obs(1:end-1)-Y_obs(2:end)'))>1e-3 ||...
mean(abs(oo_.FilteredVariables.P_obs(1:end-1)-P_obs(2:end)'))>1e-3
error('Smoothed Variables are wrong')
end
if abs(corr(oo_.FilteredVariables.Y_obs(2:end-1)-Y_obs(3:end)',oo_.FilteredVariables.Y_obs(1:end-2)-Y_obs(2:end-1)'))>2e-2 ||...
abs(corr(oo_.FilteredVariables.P_obs(2:end-1)-P_obs(3:end)',oo_.FilteredVariables.P_obs(1:end-2)-P_obs(2:end-1)'))>2e-1 ||...
abs(corr(oo_.FilteredVariables.junk2(2:end-1)-junk2(3:end)',oo_.FilteredVariables.junk2(1:end-2)-junk2(2:end-1)'))>2e-2
error('Filtered Variables are wrong')
end
if max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,1,2:end-(options_.nk-1)))-oo_.FilteredVariables.Y_obs))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,2,2:end-(options_.nk-1)))-oo_.FilteredVariables.P_obs))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(1,3,2:end-(options_.nk-1)))-oo_.FilteredVariables.junk2))>1e-5 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(2,1,3:end-options_.nk))-oo_.FilteredVariables.Y_obs(3:end)))>1e-2 ||...
max(abs(squeeze(oo_.FilteredVariablesKStepAhead(2,2,3:end-options_.nk))-oo_.FilteredVariables.P_obs(3:end)))>1e-2 ||...
mean(squeeze(oo_.FilteredVariablesKStepAhead(2,3,3:end-options_.nk)))>1e-1 ||...
mean(squeeze(oo_.FilteredVariablesKStepAhead(2,4,3:end-options_.nk)))>1e-1
error('FilteredVariablesKStepAhead is wrong')
end
if abs(mean(oo_.SmoothedShocks.e_y))>0.05 || abs(mean(oo_.SmoothedShocks.e_p))>0.05
error('Residuals are not mean 0')
end

43
tests/observation_trends_and_prefiltering/Trend_exp_model_calib_no_prefilter_common.inc Normal file
View File

@ -0,0 +1,43 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p const_y const_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
const_y=2;
const_p=2;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = exp(const_y)^(1-rho_y)*Y_obs(-1)^rho_y*exp(sigma_y*e_y);
P_obs = exp(const_p)^(1-rho_p)*P_obs(-1)^rho_p*exp(sigma_p*e_p);
junk1 = (junk1(+1))^0.9;
junk2 = (junk2(-1))^0.9*exp(eps_junk);
end;
steady_state_model;
Y_obs = exp(const_y);
P_obs = exp(const_p);
junk1=1;
junk2=1;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady(nocheck);
check;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;

41
tests/observation_trends_and_prefiltering/Trend_exp_model_calib_prefilter_common.inc Normal file
View File

@ -0,0 +1,41 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = Y_obs(-1)^rho_y*exp(sigma_y*e_y);
P_obs = P_obs(-1)^rho_p*exp(sigma_p*e_p);
junk1 = (junk1(+1))^0.9;
junk2 = (junk2(-1))^0.9*exp(eps_junk);
end;
steady_state_model;
Y_obs = 1;
P_obs = 1;
junk1=1;
junk2=1;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady(nocheck);
check;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;

62
tests/observation_trends_and_prefiltering/Trend_exp_model_no_prefilter_common.inc Normal file
View File

@ -0,0 +1,62 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p const_y const_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
const_y=2;
const_p=2;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = exp(const_y)^(1-rho_y)*Y_obs(-1)^rho_y*exp(sigma_y*e_y);
P_obs = exp(const_p)^(1-rho_p)*P_obs(-1)^rho_p*exp(sigma_p*e_p);
junk1 = (junk1(+1))^0.9;
junk2 = (junk2(-1))^0.9*exp(eps_junk);
end;
steady_state_model;
Y_obs = exp(const_y);
P_obs = exp(const_p);
junk1=1;
junk2=1;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady(nocheck);
check;
estimated_params;
const_y, normal_pdf, 2, 1;
const_p, normal_pdf, 2, 1;
g_y, normal_pdf, 0.0001, 1;
g_p, normal_pdf, -0.0001, 1;
rho_y, normal_pdf, 0.5, 1;
rho_p, normal_pdf, 0.5, 1;
sigma_y, inv_gamma_pdf, 0.001, inf;
sigma_p, inv_gamma_pdf, 0.001, inf;
end;
estimated_params_init(use_calibration);
end;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;
// estimated_params_init(use_calibration);
// end;
options_.plot_priors=0;

59
tests/observation_trends_and_prefiltering/Trend_exp_model_prefilter_common.inc Normal file
View File

@ -0,0 +1,59 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = Y_obs(-1)^rho_y*exp(sigma_y*e_y);
P_obs = P_obs(-1)^rho_p*exp(sigma_p*e_p);
junk1 = (junk1(+1))^0.9;
junk2 = (junk2(-1))^0.9*exp(eps_junk);
end;
steady_state_model;
Y_obs = 1;
P_obs = 1;
junk1=1;
junk2=1;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady(nocheck);
check;
estimated_params;
g_y, normal_pdf, 0.0001, 1;
g_p, normal_pdf, -0.0001, 1;
rho_y, normal_pdf, 0.5, 1;
rho_p, normal_pdf, 0.5, 1;
sigma_y, inv_gamma_pdf, 0.001, inf;
sigma_p, inv_gamma_pdf, 0.001, inf;
end;
estimated_params_init(use_calibration);
end;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;
// estimated_params_init(use_calibration);
// end;
options_.plot_priors=0;

11
tests/observation_trends_and_prefiltering/Trend_load_data_common.inc Normal file
View File

@ -0,0 +1,11 @@
verbatim;
if options_.loglinear
Y_obs=log(Y_obs(options_.first_obs:end));
P_obs=log(P_obs(options_.first_obs:end));
junk2=log(junk2(options_.first_obs:end));
else
Y_obs=Y_obs(options_.first_obs:end);
P_obs=P_obs(options_.first_obs:end);
junk2=junk2(options_.first_obs:end);
end
end;

43
tests/observation_trends_and_prefiltering/Trend_model_calib_no_prefilter_common.inc Normal file
View File

@ -0,0 +1,43 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p const_y const_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
const_y=2;
const_p=2;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = (1-rho_y)*const_y + rho_y*Y_obs(-1)+sigma_y*e_y;
P_obs = (1-rho_p)*const_p + rho_p*P_obs(-1)+sigma_p*e_p;
junk1 = 0.9*junk1(+1);
junk2 = 0.9*junk2(-1)+eps_junk;
end;
steady_state_model;
Y_obs = const_y;
P_obs = const_p;
junk1=0;
junk2=0;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;

43
tests/observation_trends_and_prefiltering/Trend_model_calib_prefilter_common.inc Normal file
View File

@ -0,0 +1,43 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = rho_y*Y_obs(-1)+sigma_y*e_y;
P_obs = rho_p*P_obs(-1)+sigma_p*e_p;
junk1 = 0.9*junk1(+1);
junk2 = 0.9*junk2(-1)+eps_junk;
end;
steady_state_model;
Y_obs = 0;
P_obs = 0;
junk1=0;
junk2=0;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;

58
tests/observation_trends_and_prefiltering/Trend_model_no_prefilter_common.inc Normal file
View File

@ -0,0 +1,58 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p const_y const_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
const_y=2;
const_p=2;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = (1-rho_y)*const_y + rho_y*Y_obs(-1)+sigma_y*e_y;
P_obs = (1-rho_p)*const_p + rho_p*P_obs(-1)+sigma_p*e_p;
junk1 = 0.9*junk1(+1);
junk2 = 0.9*junk2(-1)+eps_junk;
end;
steady_state_model;
Y_obs = const_y;
P_obs = const_p;
junk1=0;
junk2=0;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady;
estimated_params;
const_y, normal_pdf, 2, 1;
const_p, normal_pdf, 2, 1;
g_y, normal_pdf, 0.0001, 1;
g_p, normal_pdf, -0.0001, 1;
rho_y, normal_pdf, 0.5, 1;
rho_p, normal_pdf, 0.5, 1;
sigma_y, inv_gamma_pdf, 0.001, inf;
sigma_p, inv_gamma_pdf, 0.001, inf;
end;
varobs P_obs Y_obs junk2;
observation_trends;
Y_obs (g_y);
P_obs (g_p);
end;
estimated_params_init(use_calibration);
end;
options_.plot_priors=0;

60
tests/observation_trends_and_prefiltering/Trend_model_prefilter_common.inc Normal file
View File

@ -0,0 +1,60 @@
var Y_obs P_obs junk1 junk2;
varexo e_y e_p eps_junk;
parameters rho_y rho_p g_y g_p sigma_y sigma_p;
rho_y=0.5;
rho_p=0.5;
g_y=0.0001;
g_p=-0.0001;
sigma_y=0.001;
sigma_p=0.001;
model;
Y_obs = rho_y*Y_obs(-1)+sigma_y*e_y;
P_obs = rho_p*P_obs(-1)+sigma_p*e_p;
junk1 = 0.9*junk1(+1);
junk2 = 0.9*junk2(-1)+eps_junk;
end;
steady_state_model;
Y_obs = 0;
P_obs = 0;
junk1=0;
junk2=0;
end;
shocks;
var e_p; stderr 1;
var e_y; stderr 1;
var eps_junk; stderr 1;
end;
steady(nocheck);
check;
estimated_params;
g_y, normal_pdf, 0.0001, 1;
g_p, normal_pdf, -0.0001, 1;
rho_y, normal_pdf, 0.5, 1;
rho_p, normal_pdf, 0.5, 1;
sigma_y, inv_gamma_pdf, 0.001, inf;
sigma_p, inv_gamma_pdf, 0.001, inf;
end;
estimated_params_init(use_calibration);
end;
varobs P_obs Y_obs junk2;
observation_trends;
P_obs (g_p);
Y_obs (g_y);
end;
// estimated_params_init(use_calibration);
// end;
options_.plot_priors=0;

15
tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefil_f_obs_loglin_cal_smoother.mod Normal file
View File

@ -0,0 +1,15 @@
@#include "../Trend_exp_model_calib_no_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../Exp_AR1_trend_data_with_constant',prefilter=0,loglinear,first_obs=1000,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

16
tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilt_first_obs_cal_smooth.mod Normal file
View File

@ -0,0 +1,16 @@
@#include "../Trend_model_calib_no_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../AR1_trend_data_with_constant',prefilter=0,first_obs=1000,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

14
tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_calib_smoother.mod Normal file
View File

@ -0,0 +1,14 @@
@#include "../Trend_model_calib_no_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../AR1_trend_data_with_constant',prefilter=0,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

16
tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_loglin_calib_smoother.mod Normal file
View File

@ -0,0 +1,16 @@
@#include "../Trend_exp_model_calib_no_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../Exp_AR1_trend_data_with_constant',prefilter=0,loglinear,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params)./loaded_par.orig_params))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

15
tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod Normal file
View File

@ -0,0 +1,15 @@
@#include "../Trend_exp_model_calib_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../Exp_AR1_trend_data_with_constant',prefilter=1,loglinear,first_obs=1000,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

15
tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilt_first_obs_cal_smooth.mod Normal file
View File

@ -0,0 +1,15 @@
@#include "../Trend_model_calib_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../AR1_trend_data_with_constant',prefilter=1,first_obs=1000,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

15
tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_calib_smoother.mod Normal file
View File

@ -0,0 +1,15 @@
@#include "../Trend_model_calib_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../AR1_trend_data_with_constant',prefilter=1,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs;
load('../AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"

15
tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod Normal file
View File

@ -0,0 +1,15 @@
@#include "../Trend_exp_model_calib_prefilter_common.inc"
options_.filter_decomposition=1;
calib_smoother(datafile='../Exp_AR1_trend_data_with_constant',prefilter=1,loglinear,
// filter_decomposition,
filtered_vars, filter_step_ahead = [1,2,4]) P_obs Y_obs junk2;
load('../Exp_AR1_trend_data_with_constant');
@#include "../Trend_load_data_common.inc"
loaded_par=load('../orig_params');
if max(abs((M_.params-loaded_par.orig_params([1:4,7:8]))./loaded_par.orig_params([1:4,7:8])))>0.03
error('Parameters do not match')
end
@#include "../Trend_diagnostics_calib_common.inc"