diff --git a/doc/dynare.texi b/doc/dynare.texi
index bfdeb6be6..c06b434c0 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3977,6 +3977,12 @@ Default value: @code{1e-6}.
Saves the contemporaneous correlation between the endogenous variables in @code{oo_.contemporaneous_correlation}.
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
@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.
@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
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}
@anchor{prefilter}
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}
@anchor{presample}
@@ -4917,13 +4930,19 @@ variables
@item 2
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
For nonstationary models: use a diffuse filter (use rather the @code{diffuse_filter} option)
@item 4
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
@noindent
@@ -4986,8 +5005,12 @@ Metropolis-Hastings chain. Default: 2*@code{mh_scale}
@item mh_recover
@anchor{mh_recover} Attempts to recover a Metropolis-Hastings
-simulation that crashed prematurely. Shouldn't be used together with
-@code{load_mh_file}
+simulation that crashed prematurely, starting with the last available saved
+@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}
@dots{}
@@ -5624,8 +5647,7 @@ using a standard Kalman filter.
@xref{irf}. Only used if @ref{bayesian_irf} is passed.
@item irf_shocks = ( @var{VARIABLE_NAME} [[,] @var{VARIABLE_NAME} @dots{}] )
-@xref{irf_shocks}. Only used if @ref{bayesian_irf} is passed. Cannot be used
-with @ref{dsge_var}.
+@xref{irf_shocks}. Only used if @ref{bayesian_irf} is passed.
@item irf_plot_threshold = @var{DOUBLE}
@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
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
Herbst, Edward (2015):
``Using the ``Chandrasekhar Recursions'' for Likelihood Evaluation of DSGE
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 0140ecf12..53cfc1098 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -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.
%
% 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
%
% 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 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 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.
@@ -23,7 +23,8 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,de
% matrices (meaningless for periods 1:d)
% o decomp 4D array of shock decomposition of k-step ahead
% filtered variables
-%
+% o trend_addition [double] (n_observed_series*T) pure trend component; stored in options_.varobs order
+%
% ALGORITHM
% Diffuse Kalman filter (Durbin and Koopman)
%
@@ -93,15 +94,10 @@ else
end
trend_coeff = zeros(vobs,1);
if bayestopt_.with_trend == 1
- trend_coeff = zeros(vobs,1);
- t = options_.trend_coeffs;
- 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);
+ [trend_addition, trend_coeff] =compute_trend_coefficients(M_,options_,vobs,gend);
+ trend = constant*ones(1,gend)+trend_addition;
else
+ trend_addition=zeros(size(constant,1),gend);
trend = constant*ones(1,gend);
end
start = options_.presample+1;
diff --git a/matlab/TaRB_metropolis_hastings_core.m b/matlab/TaRB_metropolis_hastings_core.m
index 7259c8ff8..f8e5d8915 100644
--- a/matlab/TaRB_metropolis_hastings_core.m
+++ b/matlab/TaRB_metropolis_hastings_core.m
@@ -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)]);
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
- 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');
fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_chain)) 'Blck' int2str(curr_chain) ' (' datestr(now,0) ')\n']);
diff --git a/matlab/compute_trend_coefficients.m b/matlab/compute_trend_coefficients.m
new file mode 100644
index 000000000..9686867e2
--- /dev/null
+++ b/matlab/compute_trend_coefficients.m
@@ -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 .
+
+
+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
diff --git a/matlab/convert_dyn_45_to_44.m b/matlab/convert_dyn_45_to_44.m
index 3de191c90..afd795a47 100644
--- a/matlab/convert_dyn_45_to_44.m
+++ b/matlab/convert_dyn_45_to_44.m
@@ -52,7 +52,7 @@ if isfield(oo_,'PointForecast')
end
%% 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);
for ii=1:length(names)
oo_.PointForecast.HPDinf.(names{ii})=oo_.PointForecast.HPDinf.(names{ii})';
@@ -60,7 +60,7 @@ if isfield(oo_.PointForecast,'HPDinf')
end
end
-if isfield(oo_.MeanForecast,'HPDinf')
+if isfield(oo_,'MeanForecast') && isfield(oo_.MeanForecast,'HPDinf')
names=fieldnames(oo_.MeanForecast.HPDinf);
for ii=1:length(names)
oo_.MeanForecast.HPDinf.(names{ii})=oo_.MeanForecast.HPDinf.(names{ii})';
@@ -68,7 +68,7 @@ if isfield(oo_.MeanForecast,'HPDinf')
end
end
-if isfield(oo_.UpdatedVariables,'HPDinf')
+if isfield(oo_,'UpdatedVariables') && isfield(oo_.UpdatedVariables,'HPDinf')
names=fieldnames(oo_.UpdatedVariables.HPDinf);
for ii=1:length(names)
oo_.UpdatedVariables.HPDinf.(names{ii})=oo_.UpdatedVariables.HPDinf.(names{ii})';
@@ -76,7 +76,7 @@ if isfield(oo_.UpdatedVariables,'HPDinf')
end
end
-if isfield(oo_.SmoothedVariables,'HPDinf')
+if isfield(oo_,'SmoothedVariables') && isfield(oo_.SmoothedVariables,'HPDinf')
names=fieldnames(oo_.SmoothedVariables.HPDinf);
for ii=1:length(names)
oo_.SmoothedVariables.HPDinf.(names{ii})=oo_.SmoothedVariables.HPDinf.(names{ii})';
@@ -84,7 +84,7 @@ if isfield(oo_.SmoothedVariables,'HPDinf')
end
end
-if isfield(oo_.FilteredVariables,'HPDinf')
+if isfield(oo_,'FilteredVariables') && isfield(oo_.FilteredVariables,'HPDinf')
names=fieldnames(oo_.FilteredVariables.HPDinf);
for ii=1:length(names)
oo_.FilteredVariables.HPDinf.(names{ii})=oo_.FilteredVariables.HPDinf.(names{ii})';
@@ -92,7 +92,7 @@ if isfield(oo_.FilteredVariables,'HPDinf')
end
end
-if isfield(oo_.SmoothedShocks,'HPDinf')
+if isfield(oo_,'SmoothedShocks') && isfield(oo_.SmoothedShocks,'HPDinf')
names=fieldnames(oo_.SmoothedShocks.HPDinf);
for ii=1:length(names)
oo_.SmoothedShocks.HPDinf.(names{ii})=oo_.SmoothedShocks.HPDinf.(names{ii})';
@@ -100,13 +100,44 @@ if isfield(oo_.SmoothedShocks,'HPDinf')
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')
names=fieldnames(oo_.FilteredVariables);
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') ...
&& ~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)];
end
end
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index ade4a0d0d..64f972160 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -307,14 +307,8 @@ end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
- trend_coeff = zeros(DynareDataset.vobs,1);
- t = DynareOptions.trend_coeffs;
- 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];
+ [trend_addition, trend_coeff]=compute_trend_coefficients(Model,DynareOptions,DynareDataset.vobs,DynareDataset.nobs);
+ trend = repmat(constant,1,DynareDataset.nobs)+trend_addition;
else
trend_coeff = zeros(DynareDataset.vobs,1);
trend = repmat(constant,1,DynareDataset.nobs);
@@ -496,7 +490,7 @@ switch DynareOptions.lik_init
indx_unstable = find(sum(abs(V),2)>1e-5);
stable = find(sum(abs(V),2)<1e-5);
nunit = length(eigenv) - nstable;
- Pstar = options_.Harvey_scale_factor*eye(np);
+ Pstar = DynareOptions.Harvey_scale_factor*eye(nunit);
if kalman_algo ~= 2
kalman_algo = 1;
end
@@ -513,6 +507,8 @@ switch DynareOptions.lik_init
end
Pstar(stable, stable) = Pstar_tmp;
Pinf = [];
+ a = zeros(mm,1);
+ Zflag = 0;
otherwise
error('dsge_likelihood:: Unknown initialization approach for the Kalman filter!')
end
diff --git a/matlab/dyn_forecast.m b/matlab/dyn_forecast.m
index 5d8098458..58c33c96d 100644
--- a/matlab/dyn_forecast.m
+++ b/matlab/dyn_forecast.m
@@ -1,12 +1,17 @@
-function [forecast,info] = dyn_forecast(var_list,M,options,oo,task)
-% function dyn_forecast(var_list,task)
+function [forecast,info] = dyn_forecast(var_list,M,options,oo,task,dataset_info)
+% function dyn_forecast(var_list,M,options,oo,task,dataset_info)
% 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
% 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
% either 'simul' or 'smoother'
+% dataset_info: Various informations about the dataset (descriptive statistics and missing observations).
+
% OUTPUTS
% nothing is returned but the procedure saves output
% in oo_.forecast.Mean
@@ -16,7 +21,7 @@ function [forecast,info] = dyn_forecast(var_list,M,options,oo,task)
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2003-2015 Dynare Team
+% Copyright (C) 2003-2016 Dynare Team
%
% 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
% along with Dynare. If not, see .
+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;
-make_ex_;
+oo=make_ex_(M,options,oo);
maximum_lag = M.maximum_lag;
@@ -72,7 +83,20 @@ switch task
y0 = zeros(M.endo_nbr,maximum_lag);
for i = 1:M.endo_nbr
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
gend = options.nobs;
if isfield(oo.Smoother,'TrendCoeffs')
@@ -89,12 +113,13 @@ switch task
end
end
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
- global bayestopt_
- if isfield(bayestopt_,'mean_varobs')
- trend = trend + repmat(bayestopt_.mean_varobs,1,horizon+M.maximum_lag);
+ else
+ trend_coeffs=zeros(length(options.varobs),1);
end
otherwise
error('Wrong flag value')
@@ -117,10 +142,20 @@ else
options.order,var_list,M,oo,options);
end
-if ~isscalar(trend)
+if ~isscalar(trend) %add trend back to forecast
yf(i_var_obs,:) = yf(i_var_obs,:) + trend;
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
vname = deblank(var_list(i,:));
forecast.Mean.(vname) = yf(i,maximum_lag+(1:horizon))';
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 4aaf0d18b..bb1b47c34 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -73,6 +73,7 @@ end
if options_.logged_steady_state
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
+ options_.logged_steady_state=0;
end
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 58ee4fef4..77a65da5f 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -170,38 +170,8 @@ end
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
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);
- 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 = 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
+ [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_]=write_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend);
end
return
end
@@ -513,39 +483,42 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
> 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && options_.partial_information == 0 % to be fixed
%% 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);
- oo_.Smoother.SteadyState = ys;
- oo_.Smoother.TrendCoeffs = trend_coeff;
- oo_.Smoother.Variance = P;
- 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)))
- oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
- i_endo,:);
- if isfield(options_,'kalman_algo')
- 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
- 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 ~isempty(options_.nk) && options_.nk > 0 && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))
- 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
+ [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_,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.SteadyState = ys;
+% oo_.Smoother.TrendCoeffs = trend_coeff;
+% oo_.Smoother.Trend = Trend;
+% oo_.Smoother.Variance = P;
+% 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)))
+% oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
+% i_endo,:);
+% if isfield(options_,'kalman_algo')
+% 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
+% 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 ~isempty(options_.nk) && options_.nk > 0 && ~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file))
+% 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
if ~options_.nograph,
[nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
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);
end
end
- %%
- %% Smooth observational errors...
- %%
- if options_.prefilter == 1
- yf = atT(bayestopt_.mf,:)+repmat(dataset_info.descriptive.mean',1,gend);
- elseif options_.loglinear == 1
- yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+...
- trend_coeff*[1:gend];
- else
- yf = atT(bayestopt_.mf,:)+repmat(ys(bayestopt_.mfys),1,gend)+...
- trend_coeff*[1:gend];
- end
+% %%
+% %% Smooth observational errors...
+% %%
+% 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)+Trend;
+% elseif options_.loglinear == 1
+% yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+Trend;
+% else
+% yf = atT(bayestopt_.mf,:)+repmat(ys(bayestopt_.mfys),1,gend)+Trend;
+% end
if nvn
number_of_plots_to_draw = 0;
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;
index = cat(1,index,i);
end
- eval(['oo_.SmoothedMeasurementErrors.' options_.varobs{i} ' = measurement_error(i,:)'';']);
+% eval(['oo_.SmoothedMeasurementErrors.' options_.varobs{i} ' = measurement_error(i,:)'';']);
end
if ~options_.nograph
[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
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
if np > 0
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index 02d83fb62..3166f1f3a 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -375,13 +375,11 @@ if ~isfield(options_,'trend_coeffs') % No!
else% Yes!
bayestopt_.with_trend = 1;
bayestopt_.trend_coeff = {};
- trend_coeffs = options_.trend_coeffs;
- nt = length(trend_coeffs);
for i=1:options_.number_of_observed_variables
- if i > length(trend_coeffs)
+ if i > length(options_.trend_coeffs)
bayestopt_.trend_coeff{i} = '0';
else
- bayestopt_.trend_coeff{i} = trend_coeffs{i};
+ bayestopt_.trend_coeff{i} = options_.trend_coeffs{i};
end
end
end
@@ -530,9 +528,11 @@ end
[dataset_, dataset_info, newdatainterfaceflag] = makedataset(options_, options_.dsge_var*options_.dsge_varlag, gsa_flag);
-% Set options_.nobs
-options_.nobs = dataset_.nobs;
-
+%set options for old interface from the ones for new interface
+if ~isempty(dataset_)
+ options_.nobs = dataset_.nobs;
+ options_.first_obs=double(dataset_.init);
+end
% setting steadystate_check_flag option
if options_.diffuse_filter
steadystate_check_flag = 0;
@@ -556,7 +556,7 @@ if info(1)
print_info(info, 0, options_);
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;
else
options_.noconstant = 0;
diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index f18bd4e61..f610a53e8 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -42,6 +42,10 @@ options_.simul.maxit = ep.maxit;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
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;
exo_nbr = M_.exo_nbr;
maximum_lag = M_.maximum_lag;
diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m
index 76b8444c0..4348a65f2 100644
--- a/matlab/evaluate_smoother.m
+++ b/matlab/evaluate_smoother.m
@@ -87,35 +87,6 @@ if ischar(parameters)
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);
-
-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
\ No newline at end of file
+[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);
\ No newline at end of file
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index 23b820601..aec3a1b3f 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -60,7 +60,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
options);
%test whether it solves model conditional on the instruments
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)));
if ~isempty(nan_indices)
diff --git a/matlab/forcst.m b/matlab/forcst.m
index dc269e511..2684df8ee 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -17,7 +17,7 @@ function [yf,int_width]=forcst(dr,y0,horizon,var_list)
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2003-2012 Dynare Team
+% Copyright (C) 2003-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -34,9 +34,9 @@ function [yf,int_width]=forcst(dr,y0,horizon,var_list)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
-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);
nstatic = M_.nstatic;
nspred = M_.nspred;
diff --git a/matlab/forcst2.m b/matlab/forcst2.m
index 41bf18f48..efa5d463b 100644
--- a/matlab/forcst2.m
+++ b/matlab/forcst2.m
@@ -1,5 +1,16 @@
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
%
% This file is part of Dynare.
@@ -24,14 +35,10 @@ endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
ykmin_ = M_.maximum_endo_lag;
-order = options_.order;
-
k1 = [ykmin_:-1:1];
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
k2 = k2(:,1)+(ykmin_+1-k2(:,2))*endo_nbr;
-it_ = ykmin_ + 1 ;
-
% eliminate shocks with 0 variance
i_exo_var = setdiff([1:exo_nbr],find(diag(Sigma_e_) == 0));
nxs = length(i_exo_var);
diff --git a/matlab/forcst2a.m b/matlab/forcst2a.m
index ecfb665af..16fcdc143 100644
--- a/matlab/forcst2a.m
+++ b/matlab/forcst2a.m
@@ -1,6 +1,15 @@
function yf=forcst2a(y0,dr,e)
-
-% Copyright (C) 2008-2010 Dynare Team
+% function yf=forcst2a(y0,dr,e)
+% 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.
%
@@ -19,13 +28,10 @@ function yf=forcst2a(y0,dr,e)
global M_ options_
-Sigma_e_ = M_.Sigma_e;
endo_nbr = M_.endo_nbr;
-exo_nbr = M_.exo_nbr;
ykmin_ = M_.maximum_endo_lag;
horizon = size(e,1);
-order = options_.order;
k1 = [ykmin_:-1:1];
k2 = dr.kstate(find(dr.kstate(:,2) <= ykmin_+1),[1 2]);
diff --git a/matlab/gsa/read_data.m b/matlab/gsa/read_data.m
index c5851c0be..0324b8d1c 100644
--- a/matlab/gsa/read_data.m
+++ b/matlab/gsa/read_data.m
@@ -24,7 +24,7 @@ function [gend, data] = read_data()
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
-global options_ bayestopt_
+global options_
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);
end
if options_.prefilter == 1
- bayestopt_.mean_varobs = mean(rawdata,1);
- data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
+ data = transpose(rawdata-ones(gend,1)* mean(rawdata,1));
else
data = transpose(rawdata);
end
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index e418965f3..86e17fba1 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -161,10 +161,12 @@ if info(1) > 0
print_info(info, DynareOptions.noprint, DynareOptions)
end
-if any(abs(DynareResults.steady_state(BayesInfo.mfys))>1e-9) && (DynareOptions.prefilter==1)
- disp(['You are trying to estimate a model with a non zero steady state for the observed endogenous'])
- disp(['variables using demeaned data!'])
- error('You should change something in your mod file...')
+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(['variables using demeaned data!'])
+ error('You should change something in your mod file...')
+ end
end
if ~isequal(DynareOptions.mode_compute,11)
diff --git a/matlab/load_last_mh_history_file.m b/matlab/load_last_mh_history_file.m
index 3e54f9853..7dfcbb0ef 100644
--- a/matlab/load_last_mh_history_file.m
+++ b/matlab/load_last_mh_history_file.m
@@ -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
% assignin statement.
-% Copyright (C) 2013 Dynare Team
+% Copyright (C) 2013-16 Dynare Team
%
% 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.InitialSeeds = 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,'InitialLogLiK');
record = rmfield(record,'Seeds');
@@ -64,6 +65,9 @@ else
if ~isfield(record,'MCMCConcludedSuccessfully')
record.MCMCConcludedSuccessfully = NaN; % This information is forever lost...
end
+ if ~isfield(record,'MAX_nruns')
+ record.MAX_nruns=NaN(size(record.MhDraws,1),1); % This information is forever lost...
+ end
end
if isequal(nargout,0)
diff --git a/matlab/metropolis_hastings_initialization.m b/matlab/metropolis_hastings_initialization.m
index 27e8f4c74..fd6c52654 100644
--- a/matlab/metropolis_hastings_initialization.m
+++ b/matlab/metropolis_hastings_initialization.m
@@ -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_)
-%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.
%
@@ -19,16 +19,16 @@ function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck,
% o oo_ outputs structure
%
% OUTPUTS
-% o ix2 [double] (nblck*npar) vector of starting points for different chains
-% o ilogpo2 [double] (nblck*1) vector of initial posterior values for different chains
+% o ix2 [double] (NumberOfBlocks*npar) vector of starting points 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 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 fline [double] (nblck*1) vector of first draw in each chain (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 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 nblck [scalar] Number of MCM chains requested
-% o nruns [double] (nblck*1) number of draws in each chain
-% o NewFile [scalar] (nblck*1) vector storing the number
+% o NumberOfBlocks [scalar] Number of MCM chains requested
+% o nruns [double] (NumberOfBlocks*1) number of draws in each chain
+% o NewFile [scalar] (NumberOfBlocks*1) vector storing the number
% of the first MH-file to created for each chain when saving additional
% draws
% o MAX_nruns [scalar] maximum number of draws in each MH-file on the harddisk
@@ -59,10 +59,10 @@ ix2 = [];
ilogpo2 = [];
ModelName = [];
MetropolisFolder = [];
-fblck = [];
-fline = [];
+FirstBlock = [];
+FirstLine = [];
npar = [];
-nblck = [];
+NumberOfBlocks = [];
nruns = [];
NewFile = [];
MAX_nruns = [];
@@ -76,15 +76,15 @@ end
MetropolisFolder = CheckPath('metropolis',M_.dname);
BaseName = [MetropolisFolder filesep ModelName];
-nblck = options_.mh_nblck;
-nruns = ones(nblck,1)*options_.mh_replic;
+NumberOfBlocks = options_.mh_nblck;
+nruns = ones(NumberOfBlocks,1)*options_.mh_replic;
npar = length(xparam1);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
d = chol(vv);
if ~options_.load_mh_file && ~options_.mh_recover
% Here we start a new Metropolis-Hastings, previous draws are discarded.
- if nblck > 1
+ if NumberOfBlocks > 1
disp('Estimation::mcmc: Multiple chains mode.')
else
disp('Estimation::mcmc: One Chain mode.')
@@ -108,16 +108,17 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n\n');
fprintf(fidlog,'%% Session 1.\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,' \n');
- % Find initial values for the nblck chains...
- if nblck > 1% Case 1: multiple chains
+ % Find initial values for the NumberOfBlocks chains...
+ if NumberOfBlocks > 1% Case 1: multiple chains
+ set_dynare_seed('default');
fprintf(fidlog,[' Initial values of the parameters:\n']);
disp('Estimation::mcmc: Searching for initial values...')
- ix2 = zeros(nblck,npar);
- ilogpo2 = zeros(nblck,1);
- for j=1:nblck
+ ix2 = zeros(NumberOfBlocks,npar);
+ ilogpo2 = zeros(NumberOfBlocks,1);
+ for j=1:NumberOfBlocks
validate = 0;
init_iter = 0;
trial = 1;
@@ -183,22 +184,23 @@ if ~options_.load_mh_file && ~options_.mh_recover
fprintf(fidlog,' \n');
end
fprintf(fidlog,' \n');
- fblck = 1;
- fline = ones(nblck,1);
- NewFile = ones(nblck,1);
+ FirstBlock = 1;
+ FirstLine = ones(NumberOfBlocks,1);
+ NewFile = ones(NumberOfBlocks,1);
% Delete the mh-history files
delete_mh_history_files(MetropolisFolder, ModelName);
% Create a new record structure
fprintf(['Estimation::mcmc: Write details about the MCMC... ']);
AnticipatedNumberOfFiles = ceil(nruns(1)/MAX_nruns);
AnticipatedNumberOfLinesInTheLastFile = nruns(1) - (AnticipatedNumberOfFiles-1)*MAX_nruns;
- record.Nblck = nblck;
+ record.Nblck = NumberOfBlocks;
record.MhDraws = zeros(1,3);
record.MhDraws(1,1) = nruns(1);
record.MhDraws(1,2) = AnticipatedNumberOfFiles;
record.MhDraws(1,3) = AnticipatedNumberOfLinesInTheLastFile;
- record.AcceptanceRatio = zeros(1,nblck);
- for j=1:nblck
+ record.MAX_nruns=MAX_nruns;
+ 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)
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
@@ -206,8 +208,8 @@ if ~options_.load_mh_file && ~options_.mh_recover
end
record.InitialParameters = ix2;
record.InitialLogPost = ilogpo2;
- record.LastParameters = zeros(nblck,npar);
- record.LastLogPost = zeros(nblck,1);
+ record.LastParameters = zeros(NumberOfBlocks,npar);
+ record.LastLogPost = zeros(NumberOfBlocks,1);
record.LastFileNumber = AnticipatedNumberOfFiles ;
record.LastLineNumber = AnticipatedNumberOfLinesInTheLastFile;
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 lines in the last file: ' int2str(AnticipatedNumberOfLinesInTheLastFile) '.\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']);
for i=1:length(record.InitialSeeds(j).Normal)
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,['%% Session ' int2str(length(record.MhDraws(:,1))+1) '.\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,' \n');
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: 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.' ])
- nblck = past_number_of_blocks;
- options_.mh_nblck = nblck;
+ NumberOfBlocks = past_number_of_blocks;
+ options_.mh_nblck = NumberOfBlocks;
end
% I read the last line of the last mh-file for initialization of the new Metropolis-Hastings simulations:
LastFileNumber = record.LastFileNumber;
LastLineNumber = record.LastLineNumber;
if LastLineNumber < MAX_nruns
- NewFile = ones(nblck,1)*LastFileNumber;
- fline = ones(nblck,1)*(LastLineNumber+1);
+ NewFile = ones(NumberOfBlocks,1)*LastFileNumber;
+ FirstLine = ones(NumberOfBlocks,1)*(LastLineNumber+1);
else
- NewFile = ones(nblck,1)*(LastFileNumber+1);
- fline = ones(nblck,1);
+ NewFile = ones(NumberOfBlocks,1)*(LastFileNumber+1);
+ FirstLine = ones(NumberOfBlocks,1);
end
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
- fblck = 1;
+ FirstBlock = 1;
NumberOfPreviousSimulations = sum(record.MhDraws(:,1),1);
fprintf('Estimation::mcmc: I am writing a new mh-history file... ');
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,2) = AnticipatedNumberOfFiles;
record.MhDraws(end,3) = AnticipatedNumberOfLinesInTheLastFile;
+ record.MAX_nruns = [record.MAX_nruns;MAX_nruns];
record.InitialSeeds = record.LastSeeds;
write_mh_history_file(MetropolisFolder, ModelName, record);
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...
disp('Estimation::mcmc: Recover mode!')
load_last_mh_history_file(MetropolisFolder, ModelName);
- nblck = record.Nblck;% Number of "parallel" mcmc chains.
- options_.mh_nblck = nblck;
+ NumberOfBlocks = record.Nblck;% Number of "parallel" mcmc chains.
+ 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
- OldMh = 0;% The crashed metropolis was the first session.
+ OldMhExists = 0;% The crashed metropolis was the first session.
else
- OldMh = 1;% The crashed metropolis wasn't the first session.
+ OldMhExists = 1;% The crashed metropolis wasn't the first session.
end
% Default initialization:
- if OldMh
+ if OldMhExists
ilogpo2 = record.LastLogPost;
ix2 = record.LastParameters;
else
ilogpo2 = record.InitialLogPost;
ix2 = record.InitialParameters;
end
- % Set NewFile, a nblck*1 vector of integers, and fline (first line), a nblck*1 vector of integers.
- if OldMh
+ % Set NewFile, a NumberOfBlocks*1 vector of integers, and FirstLine (first line), a NumberOfBlocks*1 vector of integers.
+ % 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.
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)
- NewFile = ones(nblck,1)*LastFileNumberInThePreviousMh;
- fline = ones(nblck,1)*(LastLineNumberInThePreviousMh+1);
- else% The last mh files of the previous session are complete.
- NewFile = ones(nblck,1)*(LastFileNumberInThePreviousMh+1);
- fline = ones(nblck,1);
+ %Test if the last mh files of the previous session were not full yet
+ if LastLineNumberInThePreviousMh < MAX_nruns%not full
+ %store starting point if whole chain needs to be redone
+ NewFile = ones(NumberOfBlocks,1)*LastFileNumberInThePreviousMh;
+ FirstLine = ones(NumberOfBlocks,1)*(LastLineNumberInThePreviousMh+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
else
LastLineNumberInThePreviousMh = 0;
LastFileNumberInThePreviousMh = 0;
- NewFile = ones(nblck,1);
- fline = ones(nblck,1);
+ NewFile = ones(NumberOfBlocks,1);
+ FirstLine = ones(NumberOfBlocks,1);
+ LastFileFullIndicator=1;
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 ?
ExpectedNumberOfMhFilesPerBlock = sum(record.MhDraws(:,2),1);
- ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*nblck;
- % I count the total number of saved mh files...
+ ExpectedNumberOfMhFiles = ExpectedNumberOfMhFilesPerBlock*NumberOfBlocks;
+ % How many mh files do we actually have ?
AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
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)
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(' in the estimation command')
error('Estimation::mcmc: mh_recover option not required!')
end
- % I count the number of saved mh files per block.
- NumberOfMhFilesPerBlock = zeros(nblck,1);
- for b = 1:nblck
+ % 2. Something needs to be done; find out what
+ % Count the number of saved mh files per block.
+ NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
+ for b = 1:NumberOfBlocks
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
NumberOfMhFilesPerBlock(b) = length(BlckMhFiles);
end
- % Is there a chain with less mh files than expected ?
- while fblck <= nblck
- if NumberOfMhFilesPerBlock(fblck) < ExpectedNumberOfMhFilesPerBlock
- disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is not complete!'])
+ % Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
+ FirstBlock = 1; %initialize
+ while FirstBlock <= NumberOfBlocks
+ if NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
+ disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is not complete!'])
break
- % The mh_recover session will start from chain fblck.
+ % The mh_recover session will start from chain FirstBlock.
else
- disp(['Estimation::mcmc: Chain ' int2str(fblck) ' is complete!'])
+ disp(['Estimation::mcmc: Chain ' int2str(FirstBlock) ' is complete!'])
end
- fblck = fblck+1;
+ FirstBlock = FirstBlock+1;
end
+
+ %% 3. Overwrite default settings for
% How many mh-files are saved in this block?
- NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(fblck);
- % 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 OldMh
- NumberOfSavedMhFilesInTheCrashedBlck = NumberOfSavedMhFilesInTheCrashedBlck - LastFileNumberInThePreviousMh;
+ NumberOfSavedMhFilesInTheCrashedBlck = NumberOfMhFilesPerBlock(FirstBlock);
+ ExistingDrawsInLastMCFile=0; %initialize: no MCMC draws of current MCMC are in file from last run
+ % Check whether last present file is a file included in the last MCMC run
+ if ~LastFileFullIndicator
+ if NumberOfSavedMhFilesInTheCrashedBlck==NewFile(FirstBlock) %only that last file exists, but no files from current MCMC
+ 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
+ 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
- NumberOfSavedMhFiles = NumberOfSavedMhFilesInTheCrashedBlck+LastFileNumberInThePreviousMh;
+% % 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.
- if NumberOfSavedMhFiles
- load([BaseName '_mh' int2str(NumberOfSavedMhFiles) '_blck' int2str(fblck) '.mat']);
- ilogpo2(fblck) = logpo2(end);
- ix2(fblck,:) = x2(end,:);
+ if NumberOfSavedMhFilesInTheCrashedBlck 0
exo_init = oo_.exo_simul;
else
exo_init = zeros(size(oo_.exo_simul));
end
- make_y_;
+ oo_=make_y_(M_,options_,oo_);
end;
%exo_init
oo_.exo_simul = exo_init;
diff --git a/matlab/perfect-foresight-models/make_ex_.m b/matlab/perfect-foresight-models/make_ex_.m
index 8ad894825..73ab28f43 100644
--- a/matlab/perfect-foresight-models/make_ex_.m
+++ b/matlab/perfect-foresight-models/make_ex_.m
@@ -1,18 +1,20 @@
-function make_ex_()
+function oo_=make_ex_(M_,options_,oo_)
% forms oo_.exo_simul and oo_.exo_det_simul
%
% INPUTS
-% none
+% M_: Dynare model structure
+% options_: Dynare options structure
+% oo_: Dynare results structure
%
% OUTPUTS
-% none
+% oo_: Dynare results structure
%
% ALGORITHM
%
% SPECIAL REQUIREMENTS
%
-% Copyright (C) 1996-2014 Dynare Team
+% Copyright (C) 1996-2016 Dynare Team
%
% 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
% along with Dynare. If not, see .
-global M_ options_ oo_ ex0_
+global ex0_
if isempty(oo_.exo_steady_state)
oo_.exo_steady_state = zeros(M_.exo_nbr,1);
diff --git a/matlab/perfect-foresight-models/make_y_.m b/matlab/perfect-foresight-models/make_y_.m
index 1085b9030..1badcb5e4 100644
--- a/matlab/perfect-foresight-models/make_y_.m
+++ b/matlab/perfect-foresight-models/make_y_.m
@@ -1,18 +1,22 @@
-function make_y_()
-% function make_y_
+function oo_=make_y_(M_,options_,oo_)
+% function oo_=make_y_(M_,options_,oo_)
% forms oo_.endo_simul as guess values for deterministic simulations
%
% INPUTS
-% ...
+% M_: Dynare model structure
+% options_: Dynare options structure
+% oo_: Dynare results structure
+%
% OUTPUTS
-% ...
+% oo_: Dynare results structure
+%
% ALGORITHM
% ...
% SPECIAL REQUIREMENTS
% none
%
-% Copyright (C) 1996-2014 Dynare Team
+% Copyright (C) 1996-2016 Dynare Team
%
% 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
% along with Dynare. If not, see .
-global M_ options_ oo_ ys0_
+global ys0_
if options_.steadystate_flag
[oo_.steady_state,M_.params,check] = ...
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index 081eaf887..3b2825432 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -12,7 +12,7 @@ function perfect_foresight_setup()
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 1996-2014 Dynare Team
+% Copyright (C) 1996-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -53,8 +53,8 @@ end
if ~options_.initval_file
if isempty(options_.datafile)
- make_ex_;
- make_y_;
+ oo_=make_ex_(M_,options_,oo_);
+ oo_=make_y_(M_,options_,oo_);
else
read_data_;
end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index e171c749d..464bfff24 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -1,18 +1,18 @@
function perfect_foresight_solver()
% Computes deterministic simulations
-%
+%
% INPUTS
% None
-%
+%
% OUTPUTS
% none
-%
+%
% ALGORITHM
-%
+%
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 1996-2015 Dynare Team
+% Copyright (C) 1996-2016 Dynare Team
%
% This file is part of Dynare.
%
@@ -45,12 +45,12 @@ if options_.debug
[residual(:,ii)] = model_static(oo_.steady_state, oo_.exo_simul(ii,:),M_.params);
end
problematic_periods=find(any(isinf(residual)) | any(isnan(residual)))-M_.maximum_endo_lag;
- if ~isempty(problematic_periods)
+ if ~isempty(problematic_periods)
period_string=num2str(problematic_periods(1));
for ii=2:length(problematic_periods)
period_string=[period_string, ', ', num2str(problematic_periods(ii))];
end
- fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string);
+ fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string);
fprintf('WARNING: Check for division by 0.\n')
end
end
@@ -62,15 +62,17 @@ oo_ = perfect_foresight_solver_core(M_,options_,oo_);
% If simulation failed try homotopy.
if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
+
skipline()
disp('Simulation of the perfect foresight model failed!')
disp('Switching to a homotopy method...')
+ skipline()
+
if ~M_.maximum_lag
disp('Homotopy not implemented for purely forward models!')
disp('Failed to solve the model!')
disp('Return with empty oo_.endo_simul.')
oo_.endo_simul = [];
- skipline()
return
end
if ~M_.maximum_lead
@@ -78,25 +80,26 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
disp('Failed to solve the model!')
disp('Return with empty oo_.endo_simul.')
oo_.endo_simul = [];
- skipline()
return
end
- skipline()
+
% Disable warnings if homotopy
warning_old_state = warning;
warning off all
% Do not print anything
oldverbositylevel = options_.verbosity;
options_.verbosity = 0;
-
- exosim = oo_.exo_simul;
- exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
-
- 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
- step = .5;
+ % 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);
+
+ % Copy the current paths for the exogenous and endogenous variables.
+ exosim = oo_.exo_simul;
+ endosim = oo_.endo_simul;
+
+ current_weight = 0; % Current weight of target point in convex combination.
+ step = .5; % Set default step size.
success_counter = 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
% But take care of not overwriting the computed part of oo_.endo_simul
oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight);
-
- oo_.endo_simul(:,[initperiods, lastperiods]) = ...
- new_weight*oo_.endo_simul(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
-
+ oo_.endo_simul(:,[initperiods, lastperiods]) = new_weight*endosim(:,[initperiods, lastperiods])+(1-new_weight)*endoinit(:,[initperiods, lastperiods]);
+
+ % Detect Nans or complex numbers in the solution.
path_with_nans = any(any(isnan(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);
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);
end
-
+
+ % Make a copy of the paths.
saved_endo_simul = oo_.endo_simul;
+ % Solve for the paths of the endogenous variables.
[oo_,me] = perfect_foresight_solver_core(M_,options_,oo_);
if oo_.deterministic_simulation.status == 1
@@ -150,6 +156,7 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
end
fprintf('%i \t | %1.5f \t | %s \t | %e\n', iteration, new_weight, 'succeeded', me)
else
+ % If solver failed, then go back.
oo_.endo_simul = saved_endo_simul;
success_counter = 0;
step = step / 2;
@@ -168,11 +175,12 @@ end
if oo_.deterministic_simulation.status == 1
disp('Perfect foresight solution found.')
- skipline()
else
warning('Failed to solve perfect foresight model')
end
+skipline()
+
dyn2vec;
if ~isdates(options_.initial_period) && isnan(options_.initial_period)
@@ -182,4 +190,4 @@ else
end
ts = dseries(transpose(oo_.endo_simul),initial_period,cellstr(M_.endo_names));
-assignin('base', 'Simulated_time_series', ts);
+assignin('base', 'Simulated_time_series', ts);
\ No newline at end of file
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index a2f43343e..aa4114629 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -44,7 +44,8 @@ Y = transpose(dataset.data);
gend = dataset.nobs;
data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state;
-bayestopt_.mean_varobs = dataset_info.descriptive.mean';
+mean_varobs = dataset_info.descriptive.mean;
+
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
@@ -63,6 +64,7 @@ nvobs = length(options_.varobs);
iendo = 1:endo_nbr;
horizon = options_.forecast;
filtered_vars = options_.filtered_vars;
+IdObs = bayestopt_.mfys;
if horizon
i_last_obs = gend+(1-M_.maximum_endo_lag:0);
end
@@ -106,8 +108,6 @@ if horizon
MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
8));
- IdObs = bayestopt_.mfys;
-
end
MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
@@ -161,15 +161,16 @@ localVars.Y=Y;
localVars.data_index=data_index;
localVars.missing_value=missing_value;
localVars.varobs=options_.varobs;
+localVars.mean_varobs=mean_varobs;
localVars.irun=irun;
localVars.endo_nbr=endo_nbr;
localVars.nvn=nvn;
localVars.naK=naK;
localVars.horizon=horizon;
localVars.iendo=iendo;
+localVars.IdObs=IdObs;
if horizon
localVars.i_last_obs=i_last_obs;
- localVars.IdObs=IdObs;
localVars.MAX_nforc1=MAX_nforc1;
localVars.MAX_nforc2=MAX_nforc2;
end
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index 26b793c08..44df174a4 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -57,15 +57,16 @@ Y=myinputs.Y;
data_index=myinputs.data_index;
missing_value=myinputs.missing_value;
varobs=myinputs.varobs;
+mean_varobs=myinputs.mean_varobs;
irun=myinputs.irun;
endo_nbr=myinputs.endo_nbr;
nvn=myinputs.nvn;
naK=myinputs.naK;
horizon=myinputs.horizon;
iendo=myinputs.iendo;
+IdObs=myinputs.IdObs; %index of observables
if horizon
i_last_obs=myinputs.i_last_obs;
- IdObs=myinputs.IdObs;
MAX_nforc1=myinputs.MAX_nforc1;
MAX_nforc2=myinputs.MAX_nforc2;
end
@@ -166,10 +167,10 @@ for b=fpar:B
if run_smoother
[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);
- 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,:)+ ...
repmat(log(SteadyState(dr.order_var)),1,gend);
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,:)+ ...
repmat(SteadyState(dr.order_var),1,gend);
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;
if nvn
stock_error(:,:,irun(3)) = epsilonhat;
end
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,:);
+
+ %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
if horizon
yyyy = alphahat(iendo,i_last_obs);
yf = forcst2a(yyyy,dr,zeros(horizon,exo_nbr));
- if options_.prefilter
- yf(:,IdObs) = yf(:,IdObs)+repmat(bayestopt_.mean_varobs', ...
+ if options_.prefilter
+ % add mean
+ yf(:,IdObs) = yf(:,IdObs)+repmat(mean_varobs, ...
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
- yf(:,IdObs) = yf(:,IdObs)+(gend+[1-maxlag:horizon]')*trend_coeff';
if options_.loglinear
yf = yf+repmat(log(SteadyState'),horizon+maxlag,1);
else
@@ -203,11 +243,17 @@ for b=fpar:B
end
yf1 = forcst2(yyyy,horizon,dr,1);
if options_.prefilter == 1
+ % add mean
yf1(:,IdObs,:) = yf1(:,IdObs,:)+ ...
- repmat(bayestopt_.mean_varobs',[horizon+maxlag,1,1]);
+ repmat(mean_varobs,[horizon+maxlag,1,1]);
+ % add trend, taking into account that last point of sample is still included in forecasts and only cut off later
+ 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]);
end
- yf1(:,IdObs,:) = yf1(:,IdObs,:)+repmat((gend+[1-maxlag:horizon]')* ...
- trend_coeff',[1,1,1]);
if options_.loglinear
yf1 = yf1 + repmat(log(SteadyState'),[horizon+maxlag,1,1]);
else
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index 2c26de0e5..83bbb48aa 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -109,6 +109,7 @@ proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
block_iter=0;
for curr_block = fblck:nblck,
+ LastSeeds=[];
block_iter=block_iter+1;
try
% 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)]);
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
- 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');
fprintf(fidlog,['\n']);
fprintf(fidlog,['%% Mh' int2str(NewFile(curr_block)) 'Blck' int2str(curr_block) ' (' datestr(now,0) ')\n']);
diff --git a/matlab/reversed_extended_path.m b/matlab/reversed_extended_path.m
index 415804f40..f400e01cc 100644
--- a/matlab/reversed_extended_path.m
+++ b/matlab/reversed_extended_path.m
@@ -14,7 +14,7 @@ function innovation_paths = reversed_extended_path(controlled_variable_names, co
%
% SPECIAL REQUIREMENTS
-% Copyright (C) 2010-2011 Dynare Team.
+% Copyright (C) 2010-2016 Dynare Team.
%
% This file is part of Dynare.
%
@@ -56,10 +56,10 @@ options_.order = old_options_order;
options_.periods = 100;
% Set-up oo_.exo_simul.
-make_ex_;
+oo_=make_ex_(M_,options_,oo_);
% Set-up oo_.endo_simul.
-make_y_;
+oo_=make_y_(M_,options_,oo_);
% Get indices of the controlled endogenous variables in endo_simul.
n = length(controlled_variable_names);
diff --git a/matlab/simult_.m b/matlab/simult_.m
index a4ebe4dba..31daa3a6c 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -4,7 +4,7 @@ function y_=simult_(y0,dr,ex_,iorder)
%
% INPUTS
% 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.
% ex_ [double] T*q matrix of innovations.
% 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_(:,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 iorder==1
y_(:,1) = y_(:,1)-dr.ys;
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 750bd4929..0de085f49 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -78,14 +78,16 @@ else
if options_.logged_steady_state %if steady state was previously logged, undo this
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
+ options_.logged_steady_state=0;
end
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
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_.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
if info(1)
options_ = options_old;
diff --git a/matlab/write_smoother_results.m b/matlab/write_smoother_results.m
new file mode 100644
index 000000000..707251399
--- /dev/null
+++ b/matlab/write_smoother_results.m
@@ -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 .
+
+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
\ No newline at end of file
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index c84e3b6b3..cffa58621 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -559,13 +559,6 @@ EstimationStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
it = options_list.num_options.find("dsge_var");
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
mod_file_struct.dsge_var_calibrated = it->second;
@@ -1021,7 +1014,7 @@ ObservationTrendsStatement::ObservationTrendsStatement(const trend_elements_t &t
void
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;
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 481d8eb26..dbf7ae6f2 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -2481,6 +2481,7 @@ calib_smoother_option : o_filtered_vars
| o_prefilter
| o_loglinear
| o_first_obs
+ | o_filter_decomposition
;
extended_path : EXTENDED_PATH ';'
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 42dbdcbb6..00eb9c375 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -798,6 +798,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
<< " save('" << basename << "_results.mat', 'dataset_', '-append');" << endl << "end" << endl
<< "if exist('estimation_info', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'estimation_info', '-append');" << endl << "end" << endl
+ << "if exist('dataset_info', 'var') == 1" << endl
+ << " save('" << basename << "_results.mat', 'dataset_info', '-append');" << endl << "end" << endl
<< "if exist('oo_recursive_', 'var') == 1" << endl
<< " save('" << basename << "_results.mat', 'oo_recursive_', '-append');" << endl << "end" << endl;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 7c0c76824..6fce62bf7 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -7,6 +7,9 @@ MODFILES = \
estimation/fs2000_model_comparison.mod \
estimation/fs2000_fast.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/TaRB/fs2000_tarb.mod \
moments/example1_var_decomp.mod \
@@ -166,6 +169,14 @@ MODFILES = \
ms-sbvar/test_ms_variances_repeated_runs.mod \
ms-dsge/test_ms_dsge.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/algo1.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_missing.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/ds2.mod \
ep/rbc.mod \
@@ -260,6 +276,31 @@ MODFILES = \
differentiate_forward_vars/RBC_differentiate_forward.mod \
TeX/fs2000_corr_ME.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 \
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_uncorr_ME.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/as2007/as2007_steadystate.m \
estimation/fsdat_simul.m \
@@ -469,8 +513,22 @@ EXTRA_DIST = \
smoother2histval/fsdat_simul.m \
optimal_policy/Ramsey/find_c.m \
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/fs2000.common.inc \
+ estimation/MH_recover/fs2000.common.inc \
prior_posterior_function/posterior_function_demo.m
diff --git a/tests/estimation/MH_recover/fs2000.common.inc b/tests/estimation/MH_recover/fs2000.common.inc
new file mode 100644
index 000000000..bbcc293f9
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000.common.inc
@@ -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 .
+ */
+
+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;
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover.mod b/tests/estimation/MH_recover/fs2000_recover.mod
index 52c5664ac..3f421a8ee 100644
--- a/tests/estimation/MH_recover/fs2000_recover.mod
+++ b/tests/estimation/MH_recover/fs2000_recover.mod
@@ -1,139 +1,35 @@
-/*
- * 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.
- */
+//Test mh_recover function for RW-MH
-/*
- * 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 .
- */
+@#include "fs2000.common.inc"
-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_.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'])
+options_.MaxNumberOfBytes=1000*11*8/2;
+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'])
+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(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('fs2000_mh1_blck1.mat');
-temp2=load([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh1_blck1.mat']);
+temp1=load([M_.dname '_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')
end
%check second, affected chain with last unaffected file
-temp1=load('fs2000_mh3_blck2.mat');
-temp2=load([M_.dname filesep 'metropolis' filesep 'fs2000_recover_mh3_blck2.mat']);
+temp1=load([M_.dname '_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 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
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover_2.mod b/tests/estimation/MH_recover/fs2000_recover_2.mod
new file mode 100644
index 000000000..1d503de7c
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000_recover_2.mod
@@ -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
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover_3.mod b/tests/estimation/MH_recover/fs2000_recover_3.mod
new file mode 100644
index 000000000..3f7b920ec
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000_recover_3.mod
@@ -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
\ No newline at end of file
diff --git a/tests/estimation/MH_recover/fs2000_recover_tarb.mod b/tests/estimation/MH_recover/fs2000_recover_tarb.mod
new file mode 100644
index 000000000..9c4fca92a
--- /dev/null
+++ b/tests/estimation/MH_recover/fs2000_recover_tarb.mod
@@ -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
\ No newline at end of file
diff --git a/tests/kalman/lik_init/fs2000_common.inc b/tests/kalman/lik_init/fs2000_common.inc
new file mode 100644
index 000000000..6ab43d50a
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_common.inc
@@ -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 .
+ */
+
+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_);
\ No newline at end of file
diff --git a/tests/kalman/lik_init/fs2000_lik_init_1.mod b/tests/kalman/lik_init/fs2000_lik_init_1.mod
new file mode 100644
index 000000000..693cad39c
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_lik_init_1.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_lik_init_2.mod b/tests/kalman/lik_init/fs2000_lik_init_2.mod
new file mode 100644
index 000000000..2f96618a4
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_lik_init_2.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_lik_init_3.mod b/tests/kalman/lik_init/fs2000_lik_init_3.mod
new file mode 100644
index 000000000..d7cffb5b1
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_lik_init_3.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_lik_init_4.mod b/tests/kalman/lik_init/fs2000_lik_init_4.mod
new file mode 100644
index 000000000..eab036012
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_lik_init_4.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_lik_init_5.mod b/tests/kalman/lik_init/fs2000_lik_init_5.mod
new file mode 100644
index 000000000..9409befa9
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_lik_init_5.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_ns_common.inc b/tests/kalman/lik_init/fs2000_ns_common.inc
new file mode 100644
index 000000000..5bee05834
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_ns_common.inc
@@ -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 .
+ */
+
+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_);
\ No newline at end of file
diff --git a/tests/kalman/lik_init/fs2000_ns_lik_init_2.mod b/tests/kalman/lik_init/fs2000_ns_lik_init_2.mod
new file mode 100644
index 000000000..f2c5f1aba
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_ns_lik_init_2.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_ns_lik_init_3.mod b/tests/kalman/lik_init/fs2000_ns_lik_init_3.mod
new file mode 100644
index 000000000..4e13049da
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_ns_lik_init_3.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman/lik_init/fs2000_ns_lik_init_5.mod b/tests/kalman/lik_init/fs2000_ns_lik_init_5.mod
new file mode 100644
index 000000000..cdbb0273c
--- /dev/null
+++ b/tests/kalman/lik_init/fs2000_ns_lik_init_5.mod
@@ -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 .
+ */
+
+@#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_);
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
new file mode 100644
index 000000000..d2410c55b
--- /dev/null
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000.mod
@@ -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 .
+ */
+
+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--')
\ No newline at end of file
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
new file mode 100644
index 000000000..388d2b7a3
--- /dev/null
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML.mod
@@ -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 .
+ */
+
+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--')
+
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
new file mode 100644
index 000000000..0c4045651
--- /dev/null
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_ML_loglinear.mod
@@ -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 .
+ */
+
+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
\ No newline at end of file
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
new file mode 100644
index 000000000..88150b060
--- /dev/null
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fs2000_loglinear.mod
@@ -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 .
+ */
+
+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--')
\ No newline at end of file
diff --git a/tests/kalman_filter_smoother/compare_results_simulation/fsdat_simul_logged.m b/tests/kalman_filter_smoother/compare_results_simulation/fsdat_simul_logged.m
new file mode 100644
index 000000000..3e442115c
--- /dev/null
+++ b/tests/kalman_filter_smoother/compare_results_simulation/fsdat_simul_logged.m
@@ -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);
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_no_prefilt_first_obs_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_no_prefilt_first_obs_MC.mod
new file mode 100644
index 000000000..640630ff8
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_no_prefilt_first_obs_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_prefilt_first_obs_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_prefilt_first_obs_MC.mod
new file mode 100644
index 000000000..037704634
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglin_prefilt_first_obs_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_no_prefilter_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_no_prefilter_MC.mod
new file mode 100644
index 000000000..e39d4cf03
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_no_prefilter_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_prefilter_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_prefilter_MC.mod
new file mode 100644
index 000000000..b6e0931d6
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_loglinear_prefilter_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_MC.mod
new file mode 100644
index 000000000..7eefa594d
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_first_obs_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_first_obs_MC.mod
new file mode 100644
index 000000000..0038ef3fe
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_no_prefilter_first_obs_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_MC.mod
new file mode 100644
index 000000000..6b8259c2e
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_first_obs_MC.mod b/tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_first_obs_MC.mod
new file mode 100644
index 000000000..191655a99
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/MCMC/Trend_prefilter_first_obs_MC.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter.mod b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter.mod
new file mode 100644
index 000000000..aa659e989
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter_first_obs.mod b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter_first_obs.mod
new file mode 100644
index 000000000..be50874d7
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_no_prefilter_first_obs.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter.mod b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter.mod
new file mode 100644
index 000000000..43eefe0e6
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter_first_obs.mod b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter_first_obs.mod
new file mode 100644
index 000000000..2019dd046
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_loglinear_prefilter_first_obs.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter.mod b/tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter.mod
new file mode 100644
index 000000000..68b963112
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter_first_obs.mod b/tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter_first_obs.mod
new file mode 100644
index 000000000..4b926f1cf
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_no_prefilter_first_obs.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_prefilter.mod b/tests/observation_trends_and_prefiltering/ML/Trend_prefilter.mod
new file mode 100644
index 000000000..6dec6a3dd
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_prefilter.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/ML/Trend_prefilter_first_obs.mod b/tests/observation_trends_and_prefiltering/ML/Trend_prefilter_first_obs.mod
new file mode 100644
index 000000000..d5f45ce2e
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/ML/Trend_prefilter_first_obs.mod
@@ -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"
diff --git a/tests/observation_trends_and_prefiltering/Trend_diagnostics_MCMC_common.inc b/tests/observation_trends_and_prefiltering/Trend_diagnostics_MCMC_common.inc
new file mode 100644
index 000000000..8c6b30c81
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_diagnostics_MCMC_common.inc
@@ -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
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/Trend_diagnostics_ML_common.inc b/tests/observation_trends_and_prefiltering/Trend_diagnostics_ML_common.inc
new file mode 100644
index 000000000..dca528e27
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_diagnostics_ML_common.inc
@@ -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
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/Trend_diagnostics_calib_common.inc b/tests/observation_trends_and_prefiltering/Trend_diagnostics_calib_common.inc
new file mode 100644
index 000000000..2a7af6878
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_diagnostics_calib_common.inc
@@ -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
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/Trend_exp_model_calib_no_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_exp_model_calib_no_prefilter_common.inc
new file mode 100644
index 000000000..9deccb509
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_exp_model_calib_no_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_exp_model_calib_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_exp_model_calib_prefilter_common.inc
new file mode 100644
index 000000000..cb7cd7c6c
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_exp_model_calib_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_exp_model_no_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_exp_model_no_prefilter_common.inc
new file mode 100644
index 000000000..0b00f643e
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_exp_model_no_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_exp_model_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_exp_model_prefilter_common.inc
new file mode 100644
index 000000000..d90cd6c2e
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_exp_model_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_load_data_common.inc b/tests/observation_trends_and_prefiltering/Trend_load_data_common.inc
new file mode 100644
index 000000000..b327d53b3
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_load_data_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_model_calib_no_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_model_calib_no_prefilter_common.inc
new file mode 100644
index 000000000..423266e64
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_model_calib_no_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_model_calib_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_model_calib_prefilter_common.inc
new file mode 100644
index 000000000..e2365c28a
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_model_calib_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_model_no_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_model_no_prefilter_common.inc
new file mode 100644
index 000000000..6f7e8c609
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_model_no_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/Trend_model_prefilter_common.inc b/tests/observation_trends_and_prefiltering/Trend_model_prefilter_common.inc
new file mode 100644
index 000000000..124e6279e
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/Trend_model_prefilter_common.inc
@@ -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;
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefil_f_obs_loglin_cal_smoother.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefil_f_obs_loglin_cal_smoother.mod
new file mode 100644
index 000000000..121b6c3a9
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefil_f_obs_loglin_cal_smoother.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilt_first_obs_cal_smooth.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilt_first_obs_cal_smooth.mod
new file mode 100644
index 000000000..12ca5bcdd
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilt_first_obs_cal_smooth.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_calib_smoother.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_calib_smoother.mod
new file mode 100644
index 000000000..261f20476
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_calib_smoother.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_loglin_calib_smoother.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_loglin_calib_smoother.mod
new file mode 100644
index 000000000..db0f43279
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_no_prefilter_loglin_calib_smoother.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod
new file mode 100644
index 000000000..eb6ae5e1c
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilt_first_obs_cal_smooth.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilt_first_obs_cal_smooth.mod
new file mode 100644
index 000000000..3444502d5
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilt_first_obs_cal_smooth.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_calib_smoother.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_calib_smoother.mod
new file mode 100644
index 000000000..46a2d0310
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_calib_smoother.mod
@@ -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"
\ No newline at end of file
diff --git a/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod
new file mode 100644
index 000000000..6ff184f5b
--- /dev/null
+++ b/tests/observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod
@@ -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"
\ No newline at end of file