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

Merge branch 'use-dseries'

time-shift
Stéphane Adjemian (Charybdis) 2014-06-27 22:38:07 +02:00
parent 87b4254ce5 9c92b20f33
commit a8cb8ccd29
89 changed files with 1562 additions and 884 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
doc/dynare.texi
View File

@ -9733,7 +9733,7 @@ Instantiates an empty @dseries object, with, if defined, an initial date given b
@sp 1
@deftypefn {dseries} dseries (@var{FILENAME})
@deftypefn {dseries} dseries (@var{FILENAME}[, @var{INITIAL_DATE}])
Instantiates and populates a @dseries object with a data file specified by @var{FILENAME}, a string passed as input. Valid file types are @file{.m} file, @file{.mat} file, @file{.csv} file, and @file{.xls} file. A typical @file{.m} file will have the following form:
@ -9746,7 +9746,7 @@ azert = randn(100,1);
yuiop = randn(100,1);
@end example
If a @file{.mat} file is used instead, it should provide the same informations. Note that the @code{INIT__} variable can be either a @dates object or a string which could be used to instantiate the same @dates object.
If a @file{.mat} file is used instead, it should provide the same informations. Note that the @code{INIT__} variable can be either a @dates object or a string which could be used to instantiate the same @dates object. If @code{INIT__} is not provided in the @file{.mat} or @file{.m} file, the initial is by default set equal to @code{dates('1Y')}. If a second input argument is passed to the constructor, @dates object @var{INITIAL_DATE}, the initial date defined in @var{FILENAME} is reset to @var{INITIAL_DATE}. This is typically usefull if @code{INIT__} is not provided in the data file.
@end deftypefn

126
matlab/@dates/dates.m
View File

@ -45,7 +45,7 @@ function dd = dates(varargin) % --*-- Unitary tests --*--
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2013 Dynare Team
% Copyright (C) 2011-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -62,6 +62,7 @@ function dd = dates(varargin) % --*-- Unitary tests --*--
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Initialization.
if nargin>0 && ischar(varargin{1}) && isequal(varargin{1},'initialize')
dd = struct('ndat', 0, 'freq', NaN(0), 'time', NaN(0,2));
dd = class(dd,'dates');
@ -71,61 +72,98 @@ end
dd = evalin('base','emptydatesobject');
switch nargin
case 0
% Returns an empty object
if isequal(nargin, 0)
% Return an empty dates obect
return
case 1
if isa(varargin{1},'dates')
% Returns a copy of the input argument
dd = varargin{1};
elseif isdate(varargin{1})
date = string2date(varargin{1});
dd.ndat = 1;
dd.freq = date.freq;
dd.time = date.time;
elseif isfreq(varargin{1})
% Instantiate an empty dates object (only set frequency)
if ischar(varargin{1})
dd.freq = string2freq(varargin{1});
else
dd.freq = varargin{1};
end
end
if all(cellfun(@isdates, varargin))
% Concatenates dates in a dates object.
dd = horzcat(varargin{:});
return
end
if all(cellfun(@isstringdate,varargin))
% Concatenates dates in a dates object.
tmp = cellfun(@string2date,varargin);
if all([tmp.freq]-tmp(1).freq==0)
dd.freq = tmp(1).freq;
else
error('dates:: Wrong calling sequence of the constructor!')
error('dates::dates: Wrong calling sequence of the constructor! All dates must have common frequency.')
end
otherwise
if isdate(varargin{1})
dd.ndat = nargin;
dd.time = NaN(dd.ndat,2);
date = string2date(varargin{1});
dd.freq = date.freq;
dd.time(1,:) = date.time;
elseif isdates(varargin{1})
dd = horzcat(varargin{:});
return
elseif isfreq(varargin{1})
S.type = '()';
S.subs = varargin;
dd = subsref(dd,S);
return
dd.ndat = length(tmp);
dd.time = transpose(reshape([tmp.time],2,dd.ndat));
return
end
if isequal(nargin,1) && isfreq(varargin{1})
% Instantiate an empty dates object (only set frequency)
if ischar(varargin{1})
dd.freq = string2freq(varargin{1});
else
error('dates::dates: Wrong calling sequence!')
dd.freq = varargin{1};
end
for i=2:dd.ndat
if isdate(varargin{i})
date = string2date(varargin{i});
if isequal(date.freq,dd.freq)
dd.time(i,:) = date.time;
return
end
if isequal(nargin,3) && isfreq(varargin{1})
if ischar(varargin{1})
dd.freq = string2freq(varargin{1});
else
dd.freq = varargin{1};
end
if (isnumeric(varargin{2}) && isvector(varargin{2}) && isint(varargin{2}))
if isnumeric(varargin{3}) && isvector(varargin{3}) && isint(varargin{3})
if all(varargin{3}>=1) && all(varargin{3}<=dd.freq)
dd.time = [varargin{2}(:), varargin{3}(:)];
dd.ndat = size(dd.time,1);
else
error(sprintf('dates::dates: Check that all the inputs have the same frequency (see input number %i)!',i))
error(sprintf('dates::dates: Wrong calling sequence of the constructor! Third input must contain integers between 1 and %i.',dd.freq))
end
else
error(sprintf('dates::dates: Input %i has to be a string date!',i))
error('dates::dates: Wrong calling sequence of the constructor! Third input must be a vector of integers.')
end
else
error('dates::dates: Wrong calling sequence of the constructor! Second input must be a vector of integers.')
end
return
end
if isequal(nargin,2) && isfreq(varargin{1})
if ischar(varargin{1})
dd.freq = string2freq(varargin{1});
else
dd.freq = varargin{1};
end
if isequal(dd.freq, 1)
if (isnumeric(varargin{2}) && isvector(varargin{2}) && isint(varargin{2}))
dd.time = [varargin{2}, ones(length(varargin{2}),1)];
dd.ndat = size(dd.time,1);
return
else
error('dates::dates: Wrong calling sequence of the constructor! Second input must be a vector of integers.')
end
else
if isequal(size(varargin{2},2), 2)
if all(isint(varargin{2}(:,1))) && all(isint(varargin{2}(:,1)))
if all(varargin{2}(:,2)>=1) && all(varargin{2}(:,2)<=dd.freq)
dd.time = [varargin{2}(:,1), varargin{2}(:,2)];
dd.ndat = size(dd.time,1);
else
error(sprintf('dates::dates: Wrong calling sequence of the constructor! Second column of the last input must contain integers between 1 and %i.',dd.freq))
end
else
error('dates::dates: Wrong calling sequence! Second input argument must be an array of integers.')
end
else
error('dates::dates: Wrong calling sequence!')
end
end
return
end
error('dates::dates: Wrong calling sequence!')
%@test:1
%$ % Define some dates
%$ B1 = '1945Q3';

10
matlab/@dseries/dseries.m
View File

@ -10,7 +10,7 @@ function ts = dseries(varargin) % --*-- Unitary tests --*--
%! @sp 2
%! If @code{nargin==0} then an empty dseries object is created. The object can be populated with data subsequently using the overloaded subsref method.
%! @sp 2
%! If @code{nargin==1} and if the input argument is a @ref{dynDate} object, then a dseries object without data is created. This object can be populated with the overload subsref method.
%! If @code{nargin==1} and if the input argument is a @ref{dates} object, then a dseries object without data is created. This object can be populated with the overload subsref method.
%! @sp 2
%! If @code{nargin==1} and if the input argument is a string for the name of a csv, m or mat file containing data, then a dseries object is created from these data.
%! @sp 2
@ -54,7 +54,7 @@ function ts = dseries(varargin) % --*-- Unitary tests --*--
%! frequency is unspecified. @var{freq} is equal to 4 if data are on a quaterly basis. @var{freq} is equal to
%! 12 if data are on a monthly basis. @var{freq} is equal to 52 if data are on a weekly basis.
%! @item init
%! @ref{dynDate} object, initial date of the dataset.
%! @ref{dates} object, initial date of the dataset.
%! @end table
%! @end deftypefn
%@eod:
@ -153,6 +153,12 @@ switch nargin
ts.dates = dates(1,1):dates(1,1)+(ts.nobs-1);
end
case {2,3,4}
if isequal(nargin,2) && ischar(varargin{1}) && isdates(varargin{2})
% Instantiate dseries object with a data file and force the initial date to be as given by the second input argument.
ds = dseries(varargin{1});
ts = dseries(ds.data, varargin{2}, ds.name, ds.tex);
return
end
a = varargin{1};
b = varargin{2};
if nargin<4

2
matlab/@dseries/ydiff.m
View File

@ -151,7 +151,7 @@ end
%$
%$ try
%$ data = transpose(1:100);
%$ ts = dseries(data,1950,{'A1'},{'A_1'});
%$ ts = dseries(data,'1950Y',{'A1'},{'A_1'});
%$ ts = ts.ydiff;
%$ t(1) = 1;
%$ catch

2
matlab/@dynTimeIndex/subsasgn.m
View File

@ -17,4 +17,4 @@ function val = subsasgn(val, idx, rhs)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
error('dynTimeIndex::subsasgn: Members of dynDate class are private!')
error('dynTimeIndex::subsasgn: Members of dynTimeIndex class are private!')

2
matlab/DsgeSmoother.m
View File

@ -59,7 +59,7 @@ R = [];
P = [];
PK = [];
decomp = [];
nobs = size(options_.varobs,1);
nobs = length(options_.varobs);
smpl = size(Y,2);
M_ = set_all_parameters(xparam1,estim_params_,M_);

6
matlab/GetPosteriorParametersStatistics.m
View File

@ -180,17 +180,17 @@ if nvn
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
else
try
name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
[post_mean,hpd_interval,post_var] = Extractoo(oo_,name,type);
catch
Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
posterior_moments(Draws,1,options_.mh_conf_sig);
name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
oo_ = Filloo(oo_,name,type,post_mean,hpd_interval,post_median,post_var,post_deciles,density);
end
end

2
matlab/PlotPosteriorDistributions.m
View File

@ -89,7 +89,7 @@ for i=1:npar
eval(['pmod = oo_.posterior_mode.shocks_std.' name ';'])
end
elseif i <= nvx+nvn
name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i-nvx,1),:));
name = options_.varobs{estim_params_.nvn_observable_correspondence(i-nvx,1)};
eval(['x1 = oo_.posterior_density.measurement_errors_std.' name '(:,1);'])
eval(['f1 = oo_.posterior_density.measurement_errors_std.' name '(:,2);'])
eval(['oo_.prior_density.mearsurement_errors_std.' name '(:,1) = x2;'])

8
matlab/PosteriorFilterSmootherAndForecast.m
View File

@ -48,7 +48,7 @@ MaxNumberOfPlotPerFigure = 4;% The square root must be an integer!
MaxNumberOfBytes=options_.MaxNumberOfBytes;
endo_nbr=M_.endo_nbr;
exo_nbr=M_.exo_nbr;
nvobs = size(options_.varobs,1);
nvobs = length(options_.varobs);
nn = sqrt(MaxNumberOfPlotPerFigure);
iendo = 1:endo_nbr;
i_last_obs = gend+(1-M_.maximum_endo_lag:0);
@ -74,9 +74,9 @@ B = 200;
MAX_nruns = min(B,ceil(options_.MaxNumberOfBytes/(npar+2)/8));
MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)*gend)/8));
if naK
MAX_naK = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ...
MAX_naK = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)* ...
length(options_.filter_step_ahead)*gend)/8));
end
if horizon
@ -272,4 +272,4 @@ if options_.forecast
pm3(endo_nbr,horizon+maxlag,ifil6,B,'Forecasted variables (total)',...
M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
'names2','smooth',MetropolisFolder,'_forc_total')
end
end

76
matlab/PosteriorIRF.m
View File

@ -1,11 +1,11 @@
function PosteriorIRF(type)
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
% o type [char] string specifying the joint density of the
% deep parameters ('prior','posterior').
%
% OUTPUTS
% deep parameters ('prior','posterior').
%
% OUTPUTS
% None (oo_ and plots).
%
% SPECIAL REQUIREMENTS
@ -34,15 +34,16 @@ function PosteriorIRF(type)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ estim_params_ oo_ M_ bayestopt_ dataset_
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
% Set the number of periods
if isempty(options_.irf) || ~options_.irf
if isempty(options_.irf) || ~options_.irf
options_.irf = 40;
end
% Set varlist if necessary
varlist = options_.varlist;
if isempty(varlist)
varlist = options_.varobs;
varlist = char(options_.varobs);
end
options_.varlist = varlist;
nvar = size(varlist,1);
@ -58,7 +59,7 @@ end
% Get index of shocks for requested IRFs
irf_shocks_indx = getIrfShocksIndx();
% Set various parameters & Check or create directories
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
@ -68,8 +69,8 @@ np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np; clear('nvx','nvn','ncx','ncn','np');
nvobs = dataset_.info.nvobs;
gend = dataset_.info.ntobs;
nvobs = dataset_.vobs;
gend = dataset_.nobs;
MaxNumberOfPlotPerFigure = 9;
nn = sqrt(MaxNumberOfPlotPerFigure);
MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8);
@ -111,14 +112,14 @@ else% type = 'prior'
B = options_.prior_draws;
options_.B = B;
end
try
try
delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
catch
catch
disp('No _IRFs (dsge) files to be deleted!')
end
try
try
delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat'])
catch
catch
disp('No _IRFs (bvar-dsge) files to be deleted!')
end
irun = 0;
@ -144,10 +145,6 @@ if MAX_nirfs_dsgevar
else
stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B);
end
[mYY,mXY,mYX,mXX,Ydata,Xdata] = ...
var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,...
options_.dsge_varlag,-1,options_.datafile,options_.varobs,...
options_.xls_sheet,options_.xls_range);
NumberOfLags = options_.dsge_varlag;
NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
if options_.noconstant
@ -174,7 +171,7 @@ localVars.npar = npar;
localVars.type=type;
if strcmpi(type,'posterior')
while b<B
while b<B
b = b + 1;
x(b,:) = GetOneDraw(type);
end
@ -192,7 +189,7 @@ if options_.dsge_var
localVars.NumberOfLags = options_.dsge_varlag;
localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
end
localVars.nvar=nvar;
localVars.IndxVariables=IndxVariables;
localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
@ -225,14 +222,15 @@ else
localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
globalVars = struct('M_',M_, ...
'options_', options_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_, ...
'dataset_',dataset_);
'dataset_',dataset_, ...
'dataset_info',dataset_info);
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
@ -244,13 +242,13 @@ else
for j=1:length(fout),
nosaddle = nosaddle + fout(j).nosaddle;
end
end
% END first parallel section!
if nosaddle
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
end
ReshapeMatFiles('irf_dsge',type)
@ -323,7 +321,7 @@ if MAX_nirfs_dsgevar
MedianIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);
HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);
fprintf('Estimation::mcmc: Posterior (bvar-dsge) IRFs...\n');
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
@ -341,7 +339,7 @@ if MAX_nirfs_dsgevar
end
kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
end
clear STOCK_IRF_BVARDSGE;
clear STOCK_IRF_BVARDSGE;
for i = irf_shocks_indx
for j = 1:nvar
name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
@ -384,7 +382,7 @@ localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure;
if options_.dsge_var
localVars.HPDIRFdsgevar=HPDIRFdsgevar;
localVars.MeanIRFdsgevar = MeanIRFdsgevar;
end
end
% The files .TeX are genereted in sequential way always!
@ -394,14 +392,14 @@ if options_.TeX
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
for i=irf_shocks_indx
NAMES = [];
TEXNAMES = [];
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
name = deblank(varlist(j,:));
texname = deblank(varlist_TeX(j,:));
@ -413,7 +411,7 @@ if options_.TeX
TEXNAMES = char(TEXNAMES,['$' texname '$']);
end
end
end
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
@ -430,10 +428,10 @@ if options_.TeX
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
% The others file format are generated in parallel by PosteriorIRF_core2!
@ -453,7 +451,7 @@ if ~isoctave
else
globalVars = struct('M_',M_, ...
'options_', options_);
[fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
end
end
@ -463,4 +461,4 @@ end
% END parallel code!
fprintf('Estimation::mcmc: Posterior IRFs, done!\n');
fprintf('Estimation::mcmc: Posterior IRFs, done!\n');

50
matlab/PosteriorIRF_core1.m
View File

@ -41,7 +41,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ estim_params_ oo_ M_ bayestopt_ dataset_
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
if nargin<4,
whoiam=0;
@ -189,24 +189,24 @@ while fpar<B
end
if MAX_nirfs_dsgevar
IRUN = IRUN+1;
[fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(dataset_.info.ntobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.info.ntobs*(dsge_prior_weight+1)));
DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1)));
explosive_var = 1;
while explosive_var
% draw from the marginal posterior of SIGMA
SIGMAu_draw = rand_inverse_wishart(dataset_.info.nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
SIGMAu_draw = rand_inverse_wishart(dataset_.vobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
SIGMA_inv_upper_chol);
% draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.info.nvobs, PHI, ...
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.vobs, PHI, ...
chol(SIGMAu_draw)', chol(iXX)');
Companion_matrix(1:dataset_.info.nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
Companion_matrix(1:dataset_.vobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
% Check for stationarity
explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
end
% Get the mean
mu = zeros(1,dataset_.info.nvobs);
mu = zeros(1,dataset_.vobs);
% Get rotation
if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
@ -216,24 +216,24 @@ while fpar<B
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (options_.irf,dataset_.info.nvobs*M_.exo_nbr);
tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA;
irfs = zeros (options_.irf,dataset_.vobs*M_.exo_nbr);
tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)';
for t = 2:options_.irf
PHIpower = Companion_matrix*PHIpower;
tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA;
tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
end
tmp_dsgevar = kron(ones(options_.irf,1),mu);
for j = 1:(dataset_.info.nvobs*M_.exo_nbr)
for j = 1:(dataset_.vobs*M_.exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j));
end
end
if IRUN < MAX_nirfs_dsgevar
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr);
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr);
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
eval(['save ' instr]);
@ -242,7 +242,7 @@ while fpar<B
end
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,dataset_.info.nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
stock_irf_dsgevar = zeros(options_.irf,dataset_.vobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end
end
if irun == MAX_nirfs_dsge || irun == B || fpar == B
@ -279,20 +279,6 @@ while fpar<B
ifil2 = ifil2 + 1;
irun2 = 0;
end
% if isoctave
% if (whoiam==0),
% printf(['Posterior IRF %3.f%% done\r'],(fpar/B*100));
% end
% elseif ~whoiam,
% waitbar(fpar/B,h);
% end
% if whoiam,
% if ~isoctave
% fprintf('Done! \n');
% end
% waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(B) ' done.'];
% fMessageStatus((fpar-fpar0)/(B-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
% end
dyn_waitbar((fpar-fpar0)/(B-fpar0),h);
end
@ -310,9 +296,5 @@ end
myoutput.OutputFileName = [OutputFileName_dsge;
OutputFileName_param;
OutputFileName_bvardsge];
myoutput.nosaddle = nosaddle;

23
matlab/PosteriorIRF_core2.m
View File

@ -5,7 +5,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Perform in parallel execution a portion of the PosteriorIRF.m code.
% See also the comment in random_walk_metropolis_hastings_core.m funtion.
%
% INPUTS
% INPUTS
% See the comment in random_walk_metropolis_hastings_core.m funtion.
%
% OUTPUTS
@ -13,8 +13,8 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Contained:
% OutputFileName (i.e. the figures without the file .txt).
%
% ALGORITHM
% Portion of PosteriorIRF.m function code. Specifically the last 'for' cycle.
% ALGORITHM
% Portion of PosteriorIRF.m function code. Specifically the last 'for' cycle.
%
% SPECIAL REQUIREMENTS.
% None.
@ -36,7 +36,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ M_
global options_ M_
if nargin<4,
whoiam=0;
@ -87,7 +87,7 @@ OutputFileName={};
subplotnum = 0;
for i=fpar:npar,
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
subplotnum = subplotnum+1;
@ -96,7 +96,7 @@ for i=fpar:npar,
elseif subplotnum == 1 && ~options_.relative_irf
hh = dyn_figure(options_,'Name',['Orthogonalized shock to ' tit(i,:)]);
end
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
if ~MAX_nirfs_dsgevar
@ -108,12 +108,12 @@ for i=fpar:npar,
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min(HPDIRF(:,1,j,i)));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
hold off
else
else
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
@ -136,9 +136,9 @@ for i=fpar:npar,
else
if options_.debug
fprintf('POSTERIOR_IRF: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n',deblank(varlist(j,:)),tit(i,:),options_.impulse_responses.plot_threshold)
end
end
end
if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0)
figunumber = figunumber+1;
dyn_saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)],options_);
@ -154,9 +154,8 @@ for i=fpar:npar,
% fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
dyn_waitbar((i-fpar+1)/(npar-fpar+1),[],waitbarString);
end
end% loop over exo_var
end% loop over exo_var
myoutput.OutputFileName = OutputFileName;

4
matlab/ReshapeMatFiles.m
View File

@ -70,7 +70,7 @@ switch type
TYPEarray = 4;
case 'irf_bvardsge'
CAPtype = 'IRF_BVARDSGE';
TYPEsize = [ options_.irf , size(options_.varobs,1) , M_.exo_nbr ];
TYPEsize = [ options_.irf , length(options_.varobs) , M_.exo_nbr ];
TYPEarray = 4;
case 'smooth'
CAPtype = 'SMOOTH';
@ -82,7 +82,7 @@ switch type
TYPEarray = 3;
case 'error'
CAPtype = 'ERROR';
TYPEsize = [ size(options_.varobs,1) , options_.nobs ];
TYPEsize = [ length(options_.varobs) , options_.nobs ];
TYPEarray = 3;
case 'innov'
CAPtype = 'INNOV';

10
matlab/bvar_forecast.m
View File

@ -122,7 +122,7 @@ for i = 1:ny
dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
options_.varobs(i, :));
options_.varobs{i});
end
dyn_saveas(dyn_graph.fh,[OutputDirectoryName '/' M_.fname '_BVAR_forecast_',num2str(nlags)],options_)
@ -146,8 +146,8 @@ if ~isempty(forecast_data.realized_val)
fprintf('RMSE of BVAR(%d):\n', nlags);
for i = 1:size(options_.varobs, 1)
fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i));
for i = 1:length(options_.varobs)
fprintf('%s: %10.4f\n', options_.varobs{i}, rmse(i));
end
end
@ -162,8 +162,8 @@ if ~isdir(DirectoryName)
end
save([ DirectoryName '/simulations.mat'], 'sims_no_shock', 'sims_with_shocks');
for i = 1:size(options_.varobs, 1)
name = options_.varobs(i, :);
for i = 1:length(options_.varobs)
name = options_.varobs{i};
sims = squeeze(sims_with_shocks(:,i,:));
eval(['oo_.bvar.forecast.with_shocks.Mean.' name ' = mean(sims, 2);']);

4
matlab/bvar_irf.m
View File

@ -138,9 +138,9 @@ save([ DirectoryName '/simulations.mat'], 'sampled_irfs');
% Save results in oo_
for i=1:ny
shock_name = options_.varobs(i, :);
shock_name = options_.varobs{i};
for j=1:ny
variable_name = options_.varobs(j, :);
variable_name = options_.varobs{j};
eval(['oo_.bvar.irf.Mean.' variable_name '.' shock_name ' = posterior_mean_irfs(' int2str(j) ',' int2str(i) ',:);'])
eval(['oo_.bvar.irf.Median.' variable_name '.' shock_name ' = posterior_median_irfs(' int2str(j) ',' int2str(i) ',:);'])
eval(['oo_.bvar.irf.Var.' variable_name '.' shock_name ' = posterior_variance_irfs(' int2str(j) ',' int2str(i) ',:);'])

12
matlab/check_dsge_var_model.m
View File

@ -20,25 +20,25 @@ function check_dsge_var_model(Model, EstimatedParameters, BayesInfo)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if EstimatedParameters.nvn
error('DsgeVarLikelihood:: Measurement errors are not allowed!')
error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!')
end
if EstimatedParameters.ncn
error('DsgeVarLikelihood:: Measurement errors are not allowed!')
error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!')
end
if any(vec(Model.H))
error('DsgeVarLikelihood:: Measurement errors are not allowed!')
error('Estimation::DsgeVarLikelihood: Measurement errors are not allowed!')
end
if EstimatedParameters.ncx
error('DsgeVarLikelihood:: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.')
error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.')
end
if Model.exo_nbr>1 && any(vec(tril(Model.Sigma_e,-1)))
error('DsgeVarLikelihood:: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.')
error('Estimation::DsgeVarLikelihood: Structural innovations cannot be correlated using Dynare''s interface! Introduce the correlations in the model block instead.')
end
if isequal(BayesInfo.with_trend,1)
error('DsgeVarLikelihood :: Linear trend is not yet implemented!')
error('Estimation::DsgeVarLikelihood: Linear trend is not yet implemented!')
end

4
matlab/check_file_extension.m
View File

@ -19,6 +19,10 @@ function b = check_file_extension(file,type)
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
% Clean-up path
file = strrep(file, '../', '');
file = strrep(file, './', '');
remain = file;
while ~isempty(remain)
[ext, remain] = strtok(remain,'.');

15
matlab/check_list_of_variables.m
View File

@ -46,7 +46,7 @@ if options_.dsge_var && options_.bayesian_irf
msg = 1;
end
end
if size(varlist,1)~=size(options_.varobs)
if size(varlist,1)~=length(options_.varobs)
msg = 1;
end
if msg
@ -55,7 +55,7 @@ if options_.dsge_var && options_.bayesian_irf
skipline()
end
end
varlist = options_.varobs;
varlist = char(options_.varobs);
return
end
@ -128,7 +128,7 @@ elseif isempty(varlist) && isempty(options_.endo_vars_for_moment_computations_in
if choice==1
varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
elseif choice==2
varlist = options_.varobs;
varlist = char(options_.varobs);
elseif choice==3
varlist = NaN;
else
@ -138,15 +138,14 @@ elseif isempty(varlist) && isempty(options_.endo_vars_for_moment_computations_in
end
end
end
if isnan(varlist)
edit([M_.fname '.mod'])
end
skipline()
end
if isnan(varlist)
edit([M_.fname '.mod'])
end
skipline()
end
function format_text(remain, max_number_of_words_per_line)
index = 0;
line_of_text = [];

4
matlab/compute_moments_varendo.m
View File

@ -36,14 +36,14 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
if strcmpi(type,'posterior')
posterior = 1;
if nargin==4
var_list_ = options_.varobs;
var_list_ = char(options_.varobs);
end
elseif strcmpi(type,'prior')
posterior = 0;
if nargin==4
var_list_ = options_.prior_analysis_endo_var_list;
if isempty(var_list_)
options_.prior_analysis_var_list = options_.varobs;
options_.prior_analysis_var_list = char(options_.varobs);
end
end
else

6
matlab/display_estimation_results_table.m
View File

@ -105,7 +105,7 @@ if nvx
disp(tit1)
ip = nvx+1;
for i=1:nvn
name = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:));
name = options_.varobs{estim_params_.nvn_observable_correspondence(i,1)};
if strcmp(field_name,'posterior')
fprintf('%-*s %7.3f %8.4f %7.4f %4s %6.4f \n', ...
header_width,name,bayestopt_.p1(ip), ...
@ -276,7 +276,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
ip = nvx+1;
for i=1:nvn
idx = strmatch(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:),M_.endo_names);
idx = strmatch(options_.varobs{estim_params_.nvn_observable_correspondence(i,1)},M_.endo_names);
fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
deblank(M_.endo_names_tex(idx,:)), ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
@ -469,7 +469,7 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
ip = nvx+1;
for i=1:nvn
idx = strmatch(options_.varobs(estim_params_.nvn_observable_correspondence(i,1),:),M_.endo_names);
idx = strmatch(options_.varobs{estim_params_.nvn_observable_correspondence(i,1)},M_.endo_names);
fprintf(fidTeX,'$%s$ & %8.4f & %7.4f & %7.4f \\\\ \n',...
deblank(M_.endo_names_tex(idx,:)), ...
xparam1(ip),...

36
matlab/dsge_likelihood.m
View File

@ -1,4 +1,4 @@
function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
function [fval,DLIK,Hess,exit_flag,ys,trend_coeff,info,Model,DynareOptions,BayesInfo,DynareResults] = dsge_likelihood(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults,derivatives_info)
% Evaluates the posterior kernel of a dsge model using the specified
% kalman_algo; the resulting posterior includes the 2*pi constant of the
% likelihood function
@ -295,7 +295,7 @@ BayesInfo.mf = BayesInfo.mf1;
% Define the constant vector of the measurement equation.
if DynareOptions.noconstant
constant = zeros(DynareDataset.info.nvobs,1);
constant = zeros(DynareDataset.vobs,1);
else
if DynareOptions.loglinear
constant = log(SteadyState(BayesInfo.mfys));
@ -306,28 +306,28 @@ end
% Define the deterministic linear trend of the measurement equation.
if BayesInfo.with_trend
trend_coeff = zeros(DynareDataset.info.nvobs,1);
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.info.ntobs)+trend_coeff*[1:DynareDataset.info.ntobs];
trend = repmat(constant,1,DynareDataset.nobs)+trend_coeff*[1:DynareDataset.nobs];
else
trend = repmat(constant,1,DynareDataset.info.ntobs);
trend = repmat(constant,1,DynareDataset.nobs);
end
% Get needed informations for kalman filter routines.
start = DynareOptions.presample+1;
Z = BayesInfo.mf; % old mf
no_missing_data_flag = ~DynareDataset.missing.state;
mm = length(T); % old np
pp = DynareDataset.info.nvobs;
Z = BayesInfo.mf;
no_missing_data_flag = ~DatasetInfo.missing.state;
mm = length(T);
pp = DynareDataset.vobs;
rr = length(Q);
kalman_tol = DynareOptions.kalman_tol;
riccati_tol = DynareOptions.riccati_tol;
Y = DynareDataset.data-trend;
Y = transpose(DynareDataset.data)-trend;
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
@ -406,7 +406,7 @@ switch DynareOptions.lik_init
kalman_tol, riccati_tol, DynareOptions.presample, ...
T,R,Q,H,Z,mm,pp,rr);
else
[dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
[dLIK,dlik,a,Pstar] = missing_observations_kalman_filter_d(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
@ -441,9 +441,9 @@ switch DynareOptions.lik_init
% no need to test again for correlation elements
correlated_errors_have_been_checked = 1;
[dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,...
DynareDataset.missing.number_of_observations,...
DynareDataset.missing.no_more_missing_observations, ...
[dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DatasetInfo.missing.aindex,...
DatasetInfo.missing.number_of_observations,...
DatasetInfo.missing.no_more_missing_observations, ...
Y, 1, size(Y,2), ...
zeros(mmm,1), Pinf, Pstar, ...
kalman_tol, riccati_tol, DynareOptions.presample, ...
@ -523,7 +523,7 @@ if analytic_derivation,
DLIK = [];
AHess = [];
iv = DynareResults.dr.restrict_var_list;
if nargin<8 || isempty(derivatives_info)
if nargin<9 || isempty(derivatives_info)
[A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
if ~isempty(EstimatedParameters.var_exo)
indexo=EstimatedParameters.var_exo(:,1);
@ -655,10 +655,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
else
if 0 %DynareOptions.block
[err, LIK,lik] = block_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,...
[err, LIK,lik] = block_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,...
T,R,Q,H,Pstar,Y,start,Z,kalman_tol,riccati_tol, Model.nz_state_var, Model.n_diag, Model.nobs_non_statevar);
else
[LIK,lik] = missing_observations_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
[LIK,lik] = missing_observations_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a, Pstar, ...
kalman_tol, DynareOptions.riccati_tol, ...
DynareOptions.presample, ...
@ -733,7 +733,7 @@ if (kalman_algo==2) || (kalman_algo==4)
end
end
[LIK, lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
[LIK, lik] = univariate_kalman_filter(DatasetInfo.missing.aindex,DatasetInfo.missing.number_of_observations,DatasetInfo.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
a,Pstar, ...
DynareOptions.kalman_tol, ...
DynareOptions.riccati_tol, ...

44
matlab/dsge_var_likelihood.m
View File

@ -1,4 +1,4 @@
function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function [fval,grad,hess,exit_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(xparam1,DynareDataset,DynareInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Evaluates the posterior kernel of the bvar-dsge model.
%
% INPUTS
@ -55,7 +55,11 @@ end
nx = EstimatedParameters.nvx + EstimatedParameters.np;
% Get the number of observed variables in the VAR model.
NumberOfObservedVariables = DynareDataset.info.nvobs;
NumberOfObservedVariables = DynareDataset.vobs;
% Get the number of observations.
NumberOfObservations = DynareDataset.nobs;
% Get the number of lags in the VAR model.
NumberOfLags = DynareOptions.dsge_varlag;
@ -110,8 +114,8 @@ Model.Sigma_e = Q;
dsge_prior_weight = Model.params(dsge_prior_weight_idx);
% Is the dsge prior proper?
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs;
fval = objective_function_penalty_base+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/NumberOfObservations;
fval = objective_function_penalty_base+abs(NumberOfObservations*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
exit_flag = 0;
info = 51;
return
@ -197,9 +201,9 @@ assignin('base','GXX',GXX);
assignin('base','GYX',GYX);
if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is finite.
tmp0 = dsge_prior_weight*DynareDataset.info.ntobs*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
tmp1 = dsge_prior_weight*DynareDataset.info.ntobs*GYX + mYX;
tmp2 = inv(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX);
tmp0 = dsge_prior_weight*NumberOfObservations*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
tmp1 = dsge_prior_weight*NumberOfObservations*GYX + mYX;
tmp2 = inv(dsge_prior_weight*NumberOfObservations*GXX+mXX);
SIGMAu = tmp0 - tmp1*tmp2*tmp1'; clear('tmp0');
[SIGMAu_is_positive_definite, penalty] = ispd(SIGMAu);
if ~SIGMAu_is_positive_definite
@ -208,28 +212,28 @@ if ~isinf(dsge_prior_weight)% Evaluation of the likelihood of the dsge-var model
exit_flag = 0;
return;
end
SIGMAu = SIGMAu / (DynareDataset.info.ntobs*(1+dsge_prior_weight));
SIGMAu = SIGMAu / (NumberOfObservations*(1+dsge_prior_weight));
PHI = tmp2*tmp1'; clear('tmp1');
prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*DynareDataset.info.ntobs- ...
prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*NumberOfObservations- ...
NumberOfObservedVariables*NumberOfLags ...
+1-(1:NumberOfObservedVariables)')));
prodlng2 = sum(gammaln(.5*(dsge_prior_weight*DynareDataset.info.ntobs- ...
prodlng2 = sum(gammaln(.5*(dsge_prior_weight*NumberOfObservations- ...
NumberOfObservedVariables*NumberOfLags ...
+1-(1:NumberOfObservedVariables)')));
lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX+mXX)) ...
+ .5*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters)*log(det((dsge_prior_weight+1)*DynareDataset.info.ntobs*SIGMAu)) ...
- .5*NumberOfObservedVariables*log(det(dsge_prior_weight*DynareDataset.info.ntobs*GXX)) ...
- .5*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters)*log(det(dsge_prior_weight*DynareDataset.info.ntobs*(GYY-GYX*inv(GXX)*GYX'))) ...
+ .5*NumberOfObservedVariables*DynareDataset.info.ntobs*log(2*pi) ...
- .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*DynareDataset.info.ntobs-NumberOfParameters) ...
+ .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*DynareDataset.info.ntobs-NumberOfParameters) ...
lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX+mXX)) ...
+ .5*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters)*log(det((dsge_prior_weight+1)*NumberOfObservations*SIGMAu)) ...
- .5*NumberOfObservedVariables*log(det(dsge_prior_weight*NumberOfObservations*GXX)) ...
- .5*(dsge_prior_weight*NumberOfObservations-NumberOfParameters)*log(det(dsge_prior_weight*NumberOfObservations*(GYY-GYX*inv(GXX)*GYX'))) ...
+ .5*NumberOfObservedVariables*NumberOfObservations*log(2*pi) ...
- .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*NumberOfObservations-NumberOfParameters) ...
+ .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*NumberOfObservations-NumberOfParameters) ...
- prodlng1 + prodlng2;
else% Evaluation of the likelihood of the dsge-var model when the dsge prior weight is infinite.
iGXX = inv(GXX);
SIGMAu = GYY - GYX*iGXX*transpose(GYX);
PHI = iGXX*transpose(GYX);
lik = DynareDataset.info.ntobs * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ...
trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/DynareDataset.info.ntobs));
lik = NumberOfObservations * ( log(det(SIGMAu)) + NumberOfObservedVariables*log(2*pi) + ...
trace(inv(SIGMAu)*(mYY - transpose(mYX*PHI) - mYX*PHI + transpose(PHI)*mXX*PHI)/NumberOfObservations));
lik = .5*lik;% Minus likelihood
end
@ -254,7 +258,7 @@ if (nargout==9)
iGXX = inv(GXX);
prior.SIGMAstar = GYY - GYX*iGXX*GYX';
prior.PHIstar = iGXX*transpose(GYX);
prior.ArtificialSampleSize = fix(dsge_prior_weight*DynareDataset.info.ntobs);
prior.ArtificialSampleSize = fix(dsge_prior_weight*NumberOfObservations);
prior.DF = prior.ArtificialSampleSize - NumberOfParameters - NumberOfObservedVariables;
prior.iGXX = iGXX;
end

6
matlab/dsgevar_posterior_density.m
View File

@ -1,4 +1,4 @@
function bvar = dsgevar_posterior_density(deep,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function bvar = dsgevar_posterior_density(deep,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% This function characterizes the posterior distribution of a bvar with
% a dsge prior (as in Del Negro and Schorfheide 2003) for a given value
% of the deep parameters (structural parameters + the size of the
@ -38,7 +38,7 @@ dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
bvar.NumberOfLags = options_.varlag;
bvar.NumberOfVariables = size(options_.varobs,1);
bvar.NumberOfVariables = length(options_.varobs);
bvar.Constant = 'no';
bvar.NumberOfEstimatedParameters = bvar.NumberOfLags*bvar.NumberOfVariables;
if ~options_.noconstant
@ -47,7 +47,7 @@ if ~options_.noconstant
bvar.NumberOfVariables;
end
[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval,cost_flag,info,PHI,SIGMAu,iXX,prior] = dsge_var_likelihood(deep',DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
% Conditionnal posterior density of the lagged matrices (given Sigma) ->
% Matric-variate normal distribution.

4
matlab/dyn_forecast.m
View File

@ -81,8 +81,8 @@ switch task
order_var = oo_.dr.order_var;
i_var_obs = [];
trend_coeffs = [];
for i=1:size(var_obs,1)
tmp = strmatch(var_obs(i,:),endo_names(i_var,:),'exact');
for i=1:length(var_obs)
tmp = strmatch(var_obs{i},endo_names(i_var,:),'exact');
if ~isempty(tmp)
i_var_obs = [ i_var_obs; tmp];
trend_coeffs = [trend_coeffs; oo_.Smoother.TrendCoeffs(i)];

4
matlab/dynare_estimation.m
View File

@ -105,7 +105,7 @@ if nnobs > 1 && horizon > 0
end
endo_names = M_.endo_names;
n_varobs = size(options_.varobs,1);
n_varobs = length(options_.varobs);
if isempty(var_list)
var_list = endo_names;
@ -125,7 +125,7 @@ if nnobs > 1 && horizon > 0
IdObs = zeros(n_varobs,1);
for j=1:n_varobs
iobs = strmatch(options_.varobs(j,:),var_list,'exact');
iobs = strmatch(options_.varobs{j},var_list,'exact');
if ~isempty(iobs)
IdObs(j,1) = iobs;
end

117
matlab/dynare_estimation_1.m
View File

@ -29,7 +29,7 @@ function dynare_estimation_1(var_list_,dname)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ oo_ estim_params_ bayestopt_ dataset_
global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
% Set particle filter flag.
if options_.order > 1
@ -78,10 +78,20 @@ else
objective_function = str2func('dsge_var_likelihood');
end
[dataset_,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
[dataset_,dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_] = ...
dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
if options_.dsge_var
check_dsge_var_model(M_, estim_params_, bayestopt_);
if dataset_info.missing.state
error('Estimation::DsgeVarLikelihood: I cannot estimate a DSGE-VAR model with missing observations!')
end
if options_.noconstant
var_sample_moments(options_.dsge_varlag, -1, dataset_);
else
% The steady state is non zero ==> a constant in the VAR is needed!
var_sample_moments(options_.dsge_varlag, 0, dataset_);
end
end
%check for calibrated covariances before updating parameters
@ -114,14 +124,14 @@ if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(esti
end
data = dataset_.data;
rawdata = dataset_.rawdata;
data_index = dataset_.missing.aindex;
missing_value = dataset_.missing.state;
rawdata = dataset_.data;
data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state;
% Set number of observations
gend = options_.nobs;
gend = dataset_.nobs;
% Set the number of observed variables.
n_varobs = size(options_.varobs,1);
n_varobs = length(options_.varobs);
% Get the number of parameters to be estimated.
nvx = estim_params_.nvx; % Variance of the structural innovations (number of parameters).
nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters).
@ -154,29 +164,10 @@ if ~isempty(estim_params_)
M_ = set_all_parameters(xparam1,estim_params_,M_);
end
% compute sample moments if needed (bvar-dsge)
if options_.dsge_var
if dataset_.missing.state
error('I cannot estimate a DSGE-VAR model with missing observations!')
end
if options_.noconstant
evalin('base',...
['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
'var_sample_moments(options_.first_obs,' ...
'options_.first_obs+options_.nobs-1,options_.dsge_varlag,-1,' ...
'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);'])
else% The steady state is non zero ==> a constant in the VAR is needed!
evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
'var_sample_moments(options_.first_obs,' ...
'options_.first_obs+options_.nobs-1,options_.dsge_varlag,0,' ...
'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);'])
end
end
%% perform initial estimation checks;
try
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,M_,estim_params_,options_,bayestopt_,oo_);
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,oo_);
catch % if check fails, provide info on using calibration if present
e = lasterror();
if full_calibration_detected %calibrated model present and no explicit starting values
@ -246,7 +237,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7);
end
[xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ...
fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 2
error('ESTIMATION: mode_compute=2 option (Lester Ingber''s Adaptive Simulated Annealing) is no longer available')
case 3
@ -264,10 +255,10 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
optim_options = optimset(optim_options,'GradObj','on');
end
if ~isoctave
[xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,dataset_info_,options_,M_,estim_params_,bayestopt_,oo_);
else
% Under Octave, use a wrapper, since fminunc() does not have a 4th arg
func = @(x) objective_function(x, dataset_,options_,M_,estim_params_,bayestopt_,oo_);
func = @(x) objective_function(x, dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval,exitflag] = fminunc(func,xparam1,optim_options);
end
case 4
@ -306,7 +297,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
% Call csminwell.
[fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
% Disp value at the mode.
disp(sprintf('Objective function at mode: %f',fval))
case 5
@ -334,7 +325,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
else
nit=1000;
end
[xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,oo_);
if options_.analytic_derivation,
options_.analytic_derivation = ana_deriv;
end
@ -384,7 +375,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
end
% Evaluate the objective function.
fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
fval = feval(objective_function,xparam1,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,oo_);
OldMode = fval;
if ~exist('MeanPar','var')
MeanPar = xparam1;
@ -416,8 +407,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
flag = 'LastCall';
end
[xparam1,PostVar,Scale,PostMean] = ...
gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.mh_jscale = Scale;
mouvement = max(max(abs(PostVar-OldPostVar)));
skipline()
@ -436,11 +427,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
[xparam1,PostVar,Scale,PostMean] = ...
gmhmaxlik(objective_function,xparam1,[lb ub],...
gmhmaxlikOptions,Scale,flag,PostMean,PostVar,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
gmhmaxlikOptions,Scale,flag,PostMean,PostVar,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.mh_jscale = Scale;
mouvement = max(max(abs(PostVar-OldPostVar)));
fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
fval = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
skipline()
disp('========================================================== ')
disp([' Change in the covariance matrix = ' num2str(mouvement) '.'])
@ -473,7 +464,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
if isfield(options_,'optim_opt')
eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
end
[xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 8
% Dynare implementation of the simplex algorithm.
simplexOptions = options_.simplex;
@ -498,7 +489,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
end
end
[xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 9
% Set defaults
H0 = 1e-4*ones(nx,1);
@ -523,7 +514,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
warning('off','CMAES:NonfinitenessRange');
warning('off','CMAES:InitialSigma');
[x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
xparam1=BESTEVER.x;
disp(sprintf('\n Objective function at mode: %f',fval))
case 10
@ -554,16 +545,16 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
simpsaOptionsList = options2cell(simpsaOptions);
simpsaOptions = simpsaset(simpsaOptionsList{:});
[xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 11
options_.cova_compute = 0 ;
[xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ;
[xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_) ;
case 101
myoptions=soptions;
myoptions(2)=1e-6; %accuracy of argument
myoptions(3)=1e-6; %accuracy of function (see Solvopt p.29)
myoptions(5)= 1.0;
[xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
case 102
%simulating annealing
% LB=zeros(size(xparam1))-20;
@ -600,10 +591,10 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
disp(['c vector ' num2str(c')]);
[xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(objective_function,xparam1,maxy,rt_,epsilon,ns,nt ...
,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
otherwise
if ischar(options_.mode_compute)
[xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
else
error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!'])
end
@ -614,11 +605,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 2;
[junk1, junk2, hh] = feval(objective_function,xparam1, ...
dataset_,options_,M_,estim_params_,bayestopt_,oo_);
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.analytic_derivation = ana_deriv;
else
hh = reshape(hessian(objective_function,xparam1, ...
options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
end
end
end
@ -698,7 +689,7 @@ end
if options_.mode_check.status == 1 && ~options_.mh_posterior_mode_estimation
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 0;
mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
mode_check(objective_function,xparam1,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.analytic_derivation = ana_deriv;
end
@ -725,14 +716,14 @@ if ~options_.cova_compute
end
if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
%% display results table and store parameter estimates and standard errors in results
% display results table and store parameter estimates and standard errors in results
oo_=display_estimation_results_table(xparam1,stdh,M_,options_,estim_params_,bayestopt_,oo_,pnames,'Posterior','posterior');
%% Laplace approximation to the marginal log density:
% Laplace approximation to the marginal log density:
if options_.cova_compute
estim_params_nbr = size(xparam1,1);
scale_factor = -sum(log10(diag(invhess)));
log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess));
likelihood = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
skipline()
disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation))
@ -779,7 +770,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
ana_deriv = options_.analytic_derivation;
options_.analytic_derivation = 0;
if options_.cova_compute
feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
else
error('I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.')
end
@ -827,7 +818,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if error_flag
error('Estimation::mcmc: I cannot compute the posterior statistics!')
end
prior_posterior_statistics('posterior',dataset_);
prior_posterior_statistics('posterior',dataset_,dataset_info);
end
xparam = get_posterior_parameters('mean');
M_ = set_all_parameters(xparam,estim_params_,M_);
@ -842,7 +833,7 @@ 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_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state);
[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;
@ -949,8 +940,11 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
%% Smooth observational errors...
%%
if options_.prefilter == 1
yf = atT(bayestopt_.mf,:)+repmat(dataset_.descriptive.mean',1,gend);
yf = atT(bayestopt_.mf,:)+repmat(dataset_info.descriptive.mean',1,gend);
elseif options_.loglinear == 1
gend
bayestopt_.mfys
ys
yf = atT(bayestopt_.mf,:)+repmat(log(ys(bayestopt_.mfys)),1,gend)+...
trend_coeff*[1:gend];
else
@ -965,8 +959,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.' deblank(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);
@ -995,7 +988,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
hold on
plot(1:gend,measurement_error(index(k),:),marker_string{2,1},'linewidth',1)
hold off
name = deblank(options_.varobs(index(k),:));
name = options_.varobs{index(k)};
if gend>1
xlim([1 gend])
end
@ -1009,7 +1002,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
set(gca,'XTickLabel',options_.XTickLabel)
end
if options_.TeX
idx = strmatch(options_.varobs(index(k),:),M_.endo_names,'exact');
idx = strmatch(options_.varobs{index(k)},M_.endo_names,'exact');
texname = M_.endo_names_tex(idx,:);
if isempty(TeXNAMES)
TeXNAMES = ['$ ' deblank(texname) ' $'];
@ -1070,7 +1063,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
hold on
plot(1:gend,rawdata(:,k),marker_string{2,1},'linewidth',1)
hold off
name = deblank(options_.varobs(k,:));
name = options_.varobs{k};
if isempty(NAMES)
NAMES = name;
else
@ -1084,7 +1077,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
xlim([1 gend])
end
if options_.TeX
idx = strmatch(options_.varobs(k,:),M_.endo_names,'exact');
idx = strmatch(options_.varobs{k},M_.endo_names,'exact');
texname = M_.endo_names_tex(idx,:);
if isempty(TeXNAMES)
TeXNAMES = ['$ ' deblank(texname) ' $'];

124
matlab/dynare_estimation_init.m
View File

@ -1,4 +1,4 @@
function [dataset_, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
function [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, gsa_flag, M_, options_, oo_, estim_params_, bayestopt_)
% function dynare_estimation_init(var_list_, gsa_flag)
% preforms initialization tasks before estimation or
@ -41,47 +41,61 @@ hh = [];
if isempty(gsa_flag)
gsa_flag = 0;
else% Decide if a DSGE or DSGE-VAR has to be estimated.
else
% Decide if a DSGE or DSGE-VAR has to be estimated.
if ~isempty(strmatch('dsge_prior_weight',M_.param_names))
options_.dsge_var = 1;
end
% Get the list of the endogenous variables for which posterior statistics wil be computed
var_list_ = check_list_of_variables(options_, M_, var_list_);
options_.varlist = var_list_;
end
% Get the indices of the observed variables in M_.endo_names.
options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1);
for i = 1:size(M_.endo_names,1)
tmp = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
if ~isempty(tmp)
if length(tmp)>1
skipline()
error(['Multiple declarations of ' deblank(M_.endo_names(i,:)) ' as an observed variable is not allowed!'])
end
options_.lgyidx2varobs(i) = tmp;
% Set the number of observed variables.
options_.number_of_observed_variables = length(options_.varobs);
% Test if observed variables are declared.
if ~options_.number_of_observed_variables
error('VAROBS is missing!')
end
% Check that each declared observed variable is also an endogenous variable.
for i = 1:options_.number_of_observed_variables
id = strmatch(options_.varobs{i}, M_.endo_names, 'exact');
if isempty(id)
error(['Unknown variable (' options_.varobs{i} ')!'])
end
end
% Check that a variable is not declared as observed more than once.
if ~isequal(options_.varobs,unique(options_.varobs))
for i = 1:options_.number_of_observed_variables
if length(strmatch(options_.varobs{i},options_.varobs))>1
error(['A variable cannot be declared as observed more than once (' options_.varobs{i} ')!'])
end
end
end
% Check the perturbation order (nonlinear filters with third order perturbation, or higher order, are not yet implemented).
if options_.order>2
error(['I cannot estimate a model with a ' int2str(options_.order) ' order approximation of the model!'])
end
% Set options_.lik_init equal to 3 if diffuse filter is used or
% kalman_algo refers to a diffuse filter algorithm.
if (options_.diffuse_filter==1) || (options_.kalman_algo > 2)
if options_.lik_init == 2
error(['options diffuse_filter, lik_init and/or kalman_algo have ' ...
'contradictory settings'])
% Set options_.lik_init equal to 3 if diffuse filter is used or kalman_algo refers to a diffuse filter algorithm.
if isequal(options_.diffuse_filter,1) || (options_.kalman_algo>2)
if isequal(options_.lik_init,2)
error(['options diffuse_filter, lik_init and/or kalman_algo have contradictory settings'])
else
options_.lik_init = 3;
end
end
% If options_.lik_init == 1
% set by default options_.qz_criterium to 1-1e-6
% and check options_.qz_criterium < 1-eps if options_.lik_init == 1
% Else set by default options_.qz_criterium to 1+1e-6
if options_.lik_init == 1
% set by default options_.qz_criterium to 1-1e-6
% and check options_.qz_criterium < 1-eps if options_.lik_init == 1
% Else
% set by default options_.qz_criterium to 1+1e-6
if isequal(options_.lik_init,1)
if isempty(options_.qz_criterium)
options_.qz_criterium = 1-1e-6;
elseif options_.qz_criterium > 1-eps
@ -111,9 +125,6 @@ else
M_.dname = dname;
end
% Set the number of observed variables.
n_varobs = size(options_.varobs,1);
% Set priors over the estimated parameters.
if ~isempty(estim_params_)
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
@ -297,7 +308,10 @@ end
% storing prior parameters in results
oo_.prior.mean = bayestopt_.p1;
oo_.prior.mode = bayestopt_.p5;
oo_.prior.variance = diag(bayestopt_.p2.^2);
oo_.prior.hyperparameters.first = bayestopt_.p6;
oo_.prior.hyperparameters.second = bayestopt_.p7;
% Is there a linear trend in the measurement equation?
if ~isfield(options_,'trend_coeffs') % No!
@ -307,7 +321,7 @@ else% Yes!
bayestopt_.trend_coeff = {};
trend_coeffs = options_.trend_coeffs;
nt = length(trend_coeffs);
for i=1:n_varobs
for i=1:options_.number_of_observed_variables
if i > length(trend_coeffs)
bayestopt_.trend_coeff{i} = '0';
else
@ -326,17 +340,12 @@ nstatic = M_.nstatic; % Number of static variables.
npred = M_.nspred; % Number of predetermined variables.
nspred = M_.nspred; % Number of predetermined variables in the state equation.
% Test if observed variables are declared.
if isempty(options_.varobs)
error('VAROBS is missing')
end
% Setting resticted state space (observed + predetermined variables)
var_obs_index = [];
k1 = [];
for i=1:n_varobs
var_obs_index = [var_obs_index; strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')];
k1 = [k1; strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
for i=1:options_.number_of_observed_variables
var_obs_index = [var_obs_index; strmatch(options_.varobs{i},M_.endo_names(dr.order_var,:),'exact')];
k1 = [k1; strmatch(options_.varobs{i},M_.endo_names, 'exact')];
end
% Define union of observed and state variables
@ -359,10 +368,8 @@ k3 = [];
k3p = [];
if options_.selected_variables_only
for i=1:size(var_list_,1)
k3 = [k3; strmatch(var_list_(i,:),M_.endo_names(dr.order_var,:), ...
'exact')];
k3p = [k3; strmatch(var_list_(i,:),M_.endo_names, ...
'exact')];
k3 = [k3; strmatch(var_list_(i,:),M_.endo_names(dr.order_var,:), 'exact')];
k3p = [k3; strmatch(var_list_(i,:),M_.endo_names, 'exact')];
end
else
k3 = (1:M_.endo_nbr)';
@ -383,14 +390,12 @@ if options_.block == 1
bayestopt_.mf2 = var_obs_index;
bayestopt_.mfys = k1;
oo_.dr.restrict_columns = [size(i_posA,1)+(1:size(M_.state_var,2))];
[k2, i_posA, i_posB] = union(k3p, M_.state_var', 'rows');
bayestopt_.smoother_var_list = [k3p(i_posA); M_.state_var(sort(i_posB))'];
[junk,junk,bayestopt_.smoother_saved_var_list] = intersect(k3p,bayestopt_.smoother_var_list(:));
[junk,ic] = intersect(bayestopt_.smoother_var_list,M_.state_var);
bayestopt_.smoother_restrict_columns = ic;
[junk,bayestopt_.smoother_mf] = ismember(k1, ...
bayestopt_.smoother_var_list);
[junk,bayestopt_.smoother_mf] = ismember(k1, bayestopt_.smoother_var_list);
else
k2 = union(var_obs_index,[M_.nstatic+1:M_.nstatic+M_.nspred]', 'rows');
% Set restrict_state to postion of observed + state variables in expanded state vector.
@ -404,13 +409,11 @@ else
bayestopt_.mfys = k1;
[junk,ic] = intersect(k2,nstatic+(1:npred)');
oo_.dr.restrict_columns = [ic; length(k2)+(1:nspred-npred)'];
bayestopt_.smoother_var_list = union(k2,k3);
[junk,junk,bayestopt_.smoother_saved_var_list] = intersect(k3,bayestopt_.smoother_var_list(:));
[junk,ic] = intersect(bayestopt_.smoother_var_list,nstatic+(1:npred)');
bayestopt_.smoother_restrict_columns = ic;
[junk,bayestopt_.smoother_mf] = ismember(var_obs_index, ...
bayestopt_.smoother_var_list);
[junk,bayestopt_.smoother_mf] = ismember(var_obs_index, bayestopt_.smoother_var_list);
end;
if options_.analytic_derivation,
@ -440,37 +443,12 @@ if options_.analytic_derivation,
end
end
% Test if the data file is declared.
if isempty(options_.datafile)
if gsa_flag
dataset_ = [];
% rawdata = [];
% data_info = [];
return
else
error('datafile option is missing')
end
end
% If jscale isn't specified for an estimated parameter, use global option options_.jscale, set to 0.2, by default.
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
% Load and transform data.
transformation = [];
if options_.loglinear && ~options_.logdata
transformation = @log;
end
xls.sheet = options_.xls_sheet;
xls.range = options_.xls_range;
if ~isfield(options_,'nobs')
options_.nobs = [];
end
dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls);
options_.nobs = dataset_.info.ntobs;
% Build the dataset
[dataset_, dataset_info] = makedataset(options_, options_.dsge_var*options_.dsge_varlag, gsa_flag);
% setting steadystate_check_flag option
if options_.diffuse_filter
@ -479,6 +457,7 @@ else
steadystate_check_flag = 1;
end
% If the steady state of the observed variables is non zero, set noconstant equal 0 ()
M = M_;
nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
@ -502,5 +481,4 @@ else
disp('variables in the model block, or drop the prefiltering.')
error('The option "prefilter" is inconsistent with the non-zero mean measurement equations in the model.')
end
end
end

25
matlab/dynare_identification.m
View File

@ -128,11 +128,11 @@ options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0;
options_.plot_priors = 0;
options_.smoother=1;
[dataset_,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
[dataset_,dataset_info,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
options_ident.analytic_derivation_mode = options_.analytic_derivation_mode;
if isempty(dataset_),
dataset_.info.ntobs = periods;
dataset_.info.nvobs = rows(options_.varobs);
dataset_.info.nvobs = length(options_.varobs);
dataset_.info.varobs = options_.varobs;
dataset_.rawdata = [];
dataset_.missing.state = 0;
@ -145,17 +145,8 @@ if isempty(dataset_),
dataset_.missing.no_more_missing_observations = 1;
dataset_.descriptive.mean = [];
dataset_.data = [];
% data_info.gend = periods;
% data_info.data = [];
% data_info.data_index = [];
% data_info.number_of_observations = periods*size(options_.varobs,1);
% data_info.no_more_missing_observations = 0;
% data_info.missing_value = 0;
end
% results = prior_sampler(0,M_,bayestopt_,options_,oo_);
if prior_exist
if any(bayestopt_.pshape > 0)
if options_ident.prior_range
@ -278,7 +269,7 @@ if iload <=0,
disp('Testing current parameter values')
end
[idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
identification_analysis(params,indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1);
if info(1)~=0,
skipline()
disp('----------- ')
@ -321,7 +312,7 @@ if iload <=0,
kk=kk+1;
params = prior_draw();
[idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex,1);
end
end
if info(1)
@ -375,7 +366,7 @@ if iload <=0,
params = prior_draw();
end
[dum1, ideJ, ideH, ideGP, dum2 , info] = ...
identification_analysis(params,indx,indexo,options_MC,dataset_, prior_exist, name_tex,0);
identification_analysis(params,indx,indexo,options_MC,dataset_, dataset_info, prior_exist, name_tex,0);
if iteration==0 && info(1)==0,
MAX_tau = min(SampleSize,ceil(MaxNumberOfBytes/(size(ideH.siH,1)*nparam)/8));
stoH = zeros([size(ideH.siH,1),nparam,MAX_tau]);
@ -546,7 +537,7 @@ if SampleSize > 1,
disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
if ~iload,
[idehess_max, idemoments_max, idemodel_max, idelre_max, derivatives_info_max] = ...
identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex,1);
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_max', 'idemoments_max','idemodel_max', 'idelre_max', 'jmax', '-append');
end
disp_identification(pdraws(jmax,:), idemodel_max, idemoments_max, name,1);
@ -561,7 +552,7 @@ if SampleSize > 1,
disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
if ~iload,
[idehess_min, idemoments_min, idemodel_min, idelre_min, derivatives_info_min] = ...
identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1);
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_min', 'idemoments_min','idemodel_min', 'idelre_min', 'jmin', '-append');
end
disp_identification(pdraws(jmin,:), idemodel_min, idemoments_min, name,1);
@ -576,7 +567,7 @@ if SampleSize > 1,
disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
if ~iload,
[idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j)] = ...
identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1);
end
disp_identification(pdraws(jcrit(j),:), idemodel_(j), idemoments_(j), name,1);
close all,

24
matlab/dynare_sensitivity.m
View File

@ -84,7 +84,7 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse,
options_.mode_compute = 0;
options_.filtered_vars = 1;
options_.plot_priors = 0;
[dataset_,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
[dataset_,dataset_info,xparam1,hh, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
% computes a first linear solution to set up various variables
else
if isempty(options_.qz_criterium)
@ -154,7 +154,7 @@ options_gsa = set_default_option(options_gsa,'namlagendo',[]);
options_gsa = set_default_option(options_gsa,'namexo',[]);
% RMSE mapping
options_gsa = set_default_option(options_gsa,'lik_only',0);
options_gsa = set_default_option(options_gsa,'var_rmse',options_.varobs);
options_gsa = set_default_option(options_gsa,'var_rmse',char(options_.varobs));
options_gsa = set_default_option(options_gsa,'pfilt_rmse',0.1);
options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1);
options_gsa = set_default_option(options_gsa,'alpha_rmse',0.001);
@ -343,7 +343,7 @@ if options_gsa.rmse,
end
end
prior_posterior_statistics('gsa',dataset_);
prior_posterior_statistics('gsa',dataset_, dataset_info);
if options_.bayesian_irf
PosteriorIRF('gsa');
end
@ -366,7 +366,7 @@ if options_gsa.rmse,
end
clear a;
% filt_mc_(OutputDirectoryName,data_info);
filt_mc_(OutputDirectoryName,options_gsa,dataset_);
filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info);
end
options_.opt_gsa = options_gsa;
@ -405,9 +405,9 @@ if options_gsa.glue,
Obs.data = data;
Obs.time = [1:gend];
Obs.num = gend;
for j=1:size(options_.varobs,1)
Obs.name{j} = deblank(options_.varobs(j,:));
vj=deblank(options_.varobs(j,:));
for j=1:length(options_.varobs)
Obs.name{j} = options_.varobs{j};
vj = options_.varobs{j};
jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
js = strmatch(vj,lgy_,'exact');
@ -445,7 +445,7 @@ if options_gsa.glue,
ismoo(j)=jxj;
end
jsmoo = size(options_.varobs,1);
jsmoo = length(options_.varobs);
for j=1:M_.endo_nbr,
if ~ismember(j,ismoo),
jsmoo=jsmoo+1;
@ -470,10 +470,10 @@ if options_gsa.glue,
Exo(j).name = deblank(tit(j,:));
end
if ~options_gsa.ppost
Lik(size(options_.varobs,1)+1).name = 'logpo';
Lik(size(options_.varobs,1)+1).ini = 'yes';
Lik(size(options_.varobs,1)+1).isam = 1;
Lik(size(options_.varobs,1)+1).data = -logpo2;
Lik(length(options_.varobs)+1).name = 'logpo';
Lik(length(options_.varobs)+1).ini = 'yes';
Lik(length(options_.varobs)+1).isam = 1;
Lik(length(options_.varobs)+1).data = -logpo2;
end
Sam.name = bayestopt_.name;
Sam.dim = [size(x) 0];

21
matlab/evaluate_likelihood.m
View File

@ -37,7 +37,7 @@ function [llik,parameters] = evaluate_likelihood(parameters)
global options_ M_ bayestopt_ oo_ estim_params_
persistent dataset
persistent dataset dataset_info
if nargin==0
parameters = 'posterior mode';
@ -67,22 +67,9 @@ if ischar(parameters)
end
if isempty(dataset)
% Load and transform data.
transformation = [];
if options_.loglinear && ~options_.logdata
transformation = @log;
end
xls.sheet = options_.xls_sheet;
xls.range = options_.xls_range;
if ~isfield(options_,'nobs')
options_.nobs = [];
end
dataset = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls);
[dataset, dataset_info] = makedataset(options_);
end
llik = -dsge_likelihood(parameters,dataset,options_,M_,estim_params_,bayestopt_,oo_);
llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
ldens = evaluate_prior(parameters);
llik = llik - ldens;
llik = llik - ldens;

6
matlab/evaluate_smoother.m
View File

@ -44,14 +44,14 @@ function oo_=evaluate_smoother(parameters,var_list)
global options_ M_ bayestopt_ oo_ estim_params_ % estim_params_ may be emty
persistent dataset_
persistent dataset_ dataset_info
if ischar(parameters) && strcmp(parameters,'calibration')
options_.smoother=1;
end
if isempty(dataset_) || isempty(bayestopt_)
[dataset_,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list, M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_);
[dataset_,dataset_info,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list, M_.fname, [], M_, options_, oo_, estim_params_, bayestopt_);
end
if nargin==0
@ -88,7 +88,7 @@ if ischar(parameters)
end
[atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = ...
DsgeSmoother(parameters,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state);
DsgeSmoother(parameters,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state);
oo_.Smoother.SteadyState = ys;
oo_.Smoother.TrendCoeffs = trend_coeff;

11
matlab/forecast_graphs.m
View File

@ -24,10 +24,6 @@ nr = 3;
exo_nbr = M_.exo_nbr;
endo_names = M_.endo_names;
fname = M_.fname;
% $$$ varobs = options_.varobs;
% $$$ y = oo_.SmoothedVariables;
% $$$ ys = oo_.dr.ys;
% $$$ gend = size(y,2);
yf = oo_.forecast.Mean;
hpdinf = oo_.forecast.HPDinf;
hpdsup = oo_.forecast.HPDsup;
@ -44,13 +40,6 @@ for i = 1:size(var_list)
end
nvar = length(i_var);
% $$$ % build trend for smoothed variables if necessary
% $$$ trend = zeros(size(varobs,1),10);
% $$$ if isfield(oo_.Smoother,'TrendCoeffs')
% $$$ trend_coeffs = oo_.Smoother.TrendCoeffs;
% $$$ trend = trend_coeffs*(gend-9:gend);
% $$$ end
% create subdirectory <fname>/graphs if id doesn't exist
if ~exist(fname, 'dir')
mkdir('.',fname);

42
matlab/get_file_extension.m Normal file
View File

@ -0,0 +1,42 @@
function ext = get_file_extension(file)
% returns the extension of a file.
%
% INPUTS
% o file string, name of the file
%
% OUTPUTS
% o ext string, extension.
%
% REMARKS
% If the provided file name has no extension, the routine will return an empty array.
% Copyright (C) 2013 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Clean-up path
file = strrep(file, '../', '');
file = strrep(file, './', '');
remain = file;
while ~isempty(remain)
[ext, remain] = strtok(remain,'.');
end
if strcmp(ext,file)
ext = [];
end

6
matlab/get_posterior_parameters.m
View File

@ -50,7 +50,7 @@ end
for i=1:nvn
k1 = estim_params_.nvn_observable_correspondence(i,1);
name1 = deblank(options_.varobs(k1,:));
name1 = options_.varobs{k1};
xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_std.' name1]);
m = m+1;
end
@ -69,8 +69,8 @@ end
for i=1:ncn
k1 = estim_params_.corrn_observable_correspondence(i,1);
k2 = estim_params_.corrn_observable_correspondence(i,2);
name1 = deblank(options_.varobs(k1,:));
name2 = deblank(options_.varobs(k2,:));
name1 = options_.varobs{k1};
name2 = options_.varobs{k2};
xparam(m) = eval(['oo_.posterior_' type '.measurement_errors_corr.' name1 '_' name2]);
m = m+1;
end

2
matlab/get_the_name.m
View File

@ -73,7 +73,7 @@ if k <= nvx
texnam = ['$ SE_{' tname '} $'];
end
elseif k <= (nvx+nvn)
vname = deblank(options_.varobs(estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1),:));
vname = options_.varobs{estim_params_.nvn_observable_correspondence(k-estim_params_.nvx,1)};
nam=['SE_EOBS_',vname];
if TeX
tname = deblank(M_.endo_names_tex(estim_params_.var_endo(k-estim_params_.nvx,1),:));

2
matlab/get_variables_list.m
View File

@ -35,7 +35,7 @@ function [ivar,vartan,options_] = get_variables_list(options_,M_)
varlist = options_.varlist;
if isempty(varlist)
varlist = options_.varobs;
varlist = char(options_.varobs);
options_.varlist = varlist;
end
nvar = rows(varlist);

14
matlab/global_initialization.m
View File

@ -33,6 +33,7 @@ global oo_ M_ options_ estim_params_ bayestopt_ estimation_info ex0_ ys0_
estim_params_ = [];
bayestopt_ = [];
options_.datafile = '';
options_.dataset = [];
options_.verbosity = 1;
options_.terminal_condition = 0;
options_.rplottype = 0;
@ -349,11 +350,13 @@ estimation_info.structural_innovation_corr_prior_index = {};
estimation_info.structural_innovation_corr_options_index = {};
estimation_info.structural_innovation_corr.range_index = {};
options_.initial_period = dates(1,1);
options_.dataset.firstobs = options_.initial_period;
options_.dataset.lastobs = NaN;
options_.dataset.file = [];
options_.dataset.series = [];
options_.dataset.firstobs = dates();
options_.dataset.lastobs = dates();
options_.dataset.nobs = NaN;
options_.dataset.xls_sheet = NaN;
options_.dataset.xls_range = NaN;
options_.dataset.xls_sheet = [];
options_.dataset.xls_range = [];
options_.Harvey_scale_factor = 10;
options_.MaxNumberOfBytes = 1e6;
options_.MaximumNumberOfMegaBytes = 111;
@ -364,7 +367,8 @@ options_.bayesian_th_moments = 0;
options_.diffuse_filter = 0;
options_.filter_step_ahead = [];
options_.filtered_vars = 0;
options_.first_obs = 1;
options_.first_obs = NaN;
options_.nobs = NaN;
options_.kalman_algo = 0;
options_.fast_kalman = 0;
options_.kalman_tol = 1e-10;

10
matlab/gsa/filt_mc_.m
View File

@ -1,4 +1,4 @@
function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_)
function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% function [rmse_MC, ixx] = filt_mc_(OutDir)
% inputs (from opt_gsa structure)
% vvarvecm = options_gsa_.var_rmse;
@ -120,10 +120,10 @@ if ~loadSA,
ys_mean=steady_(M_,options_,oo_);
end
% eval(options_.datafile)
Y = dataset_.data;
gend = dataset_.info.ntobs;
data_index = dataset_.missing.aindex;
missing_value = dataset_.missing.state;
Y = transpose(dataset_.data);
gend = dataset_.nobs;
data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state;
for jx=1:gend, data_indx(jx,data_index{jx})=true; end
%stock_gend=data_info.gend;
%stock_data = data_info.data;

11
matlab/identification_analysis.m
View File

@ -1,4 +1,4 @@
function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init)
function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex, init)
% function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init)
% given the parameter vector params, wraps all identification analyses
%
@ -134,17 +134,18 @@ if info(1)==0,
options_.noprint = 1;
options_.order = 1;
options_.SpectralDensity.trigger = 0;
options_.periods = data_info.info.ntobs+100;
options_.periods = dataset_.nobs+100;
if options_.kalman_algo > 2,
options_.kalman_algo = 1;
end
analytic_derivation = options_.analytic_derivation;
options_.analytic_derivation = -2;
info = stoch_simul(options_.varobs);
data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);
info = stoch_simul(char(options_.varobs));
dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end),dataset_.dates(1),dataset_.names,dataset_.tex);
%data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);
% datax=data;
derivatives_info.no_DLIK=1;
[fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
[fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
% fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
options_.analytic_derivation = analytic_derivation;
AHess=-AHess;

28
matlab/imcforecast.m
View File

@ -102,30 +102,13 @@ if estimated_model
error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
end
end
set_parameters(xparam);
% Load and transform data.
transformation = [];
if options_.loglinear && ~options_.logdata
transformation = @log;
end
xls.sheet = options_.xls_sheet;
xls.range = options_.xls_range;
if ~isfield(options_,'nobs')
options_.nobs = [];
end
dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls);
data = dataset_.data;
data_index = dataset_.missing.aindex;
gend = options_.nobs;
missing_value = dataset_.missing.state;
[dataset_,dataset_info] = makedataset(options_);
data = transpose(dataset_.data);
data_index = dataset_info.missing.aindex;
gend = dataset_.nobs;
missing_value = dataset_info.missing.state;
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data,data_index,missing_value);
trend = repmat(ys,1,options_cond_fcst.periods+1);
for i=1:M_.endo_nbr
j = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact');
@ -134,7 +117,6 @@ if estimated_model
end
end
trend = trend(oo_.dr.order_var,:);
InitState(:,1) = atT(:,end);
else
InitState(:,1) = zeros(M_.endo_nbr,1);

13
matlab/initial_estimation_checks.m
View File

@ -1,11 +1,12 @@
function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,Model,EstimatedParameters,DynareOptions,BayesInfo,DynareResults)
function DynareResults = initial_estimation_checks(objective_function,xparam1,DynareDataset,DatasetInfo,Model,EstimatedParameters,DynareOptions,BayesInfo,DynareResults)
% function initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% Checks data (complex values, ML evaluation, initial values, BK conditions,..)
%
% INPUTS
% objective_function [function handle] of the objective function
% xparam1: [vector] of parameters to be estimated
% DynareDataset: [structure] storing the dataset information
% DynareDataset: [dseries] object storing the dataset
% DataSetInfo: [structure] storing informations about the sample.
% Model: [structure] decribing the model
% EstimatedParameters [structure] characterizing parameters to be estimated
% DynareOptions [structure] describing the options
@ -35,11 +36,11 @@ function DynareResults = initial_estimation_checks(objective_function,xparam1,Dy
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if DynareDataset.info.nvobs>Model.exo_nbr+EstimatedParameters.nvn
if DynareDataset.vobs>Model.exo_nbr+EstimatedParameters.nvn
error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables!'])
end
if DynareDataset.info.nvobs>length(find(diag(Model.Sigma_e)))+EstimatedParameters.nvn
if DynareDataset.vobs>length(find(diag(Model.Sigma_e)))+EstimatedParameters.nvn
error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
end
@ -72,7 +73,7 @@ end
% Evaluate the likelihood.
ana_deriv = DynareOptions.analytic_derivation;
DynareOptions.analytic_derivation=0;
[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
DynareOptions.analytic_derivation=ana_deriv;
if DynareOptions.dsge_var || strcmp(func2str(objective_function),'non_linear_dsge_likelihood')

2
matlab/load_m_file_data.m
View File

@ -36,7 +36,7 @@ function [freq,init,data,varlist,tex] = load_m_file_data(file)
if isoctave
run(file);
else
[basename, ext] = strtok(file,'.');
basename = file(1:end-2);
run(basename);
end

8
matlab/metropolis_hastings_initialization.m
View File

@ -1,7 +1,7 @@
function [ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
%function [ ix2, ilogpo2, ModelName, MhDirectoryName, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] =
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,options_,M_,estim_params_,bayestopt_,oo_)
% metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds, dataset_,dataset_info,,options_,M_,estim_params_,bayestopt_,oo_)
% Metropolis-Hastings initialization.
%
% INPUTS
@ -110,7 +110,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2))
ix2(j,:) = candidate;
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
if ~isfinite(ilogpo2(j)) % if returned log-density is
% Inf or Nan (penalized value)
validate = 0;
@ -152,7 +152,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
candidate = transpose(xparam1(:));%
if all(candidate(:) > mh_bounds(:,1)) && all(candidate(:) < mh_bounds(:,2))
ix2 = candidate;
ilogpo2 = - feval(TargetFun,ix2',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
ilogpo2 = - feval(TargetFun,ix2',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
disp('Estimation::mcmc: Initialization at the posterior mode.')
skipline()
fprintf(fidlog,[' Blck ' int2str(1) 'params:\n']);

6
matlab/mode_check.m
View File

@ -1,4 +1,4 @@
function mode_check(fun,x,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function mode_check(fun,x,hessian,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
% Checks the estimated ML mode or Posterior mode.
%@info:
@ -62,7 +62,7 @@ if ~isempty(hessian);
[ s_min, k ] = min(diag(hessian));
end
fval = feval(fun,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
fval = feval(fun,x,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
if ~isempty(hessian);
skipline()
@ -138,7 +138,7 @@ for plt = 1:nbplt,
end
for i=1:length(z)
xx(kk) = z(i);
[fval, junk1, junk2, exit_flag] = feval(fun,xx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval, junk1, junk2, exit_flag] = feval(fun,xx,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
if exit_flag
y(i,1) = fval;
else

24
matlab/mr_gstep.m
View File

@ -1,9 +1,17 @@
function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin)
% function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin)
%
% Gibbs type step in optimisation
%
% varargin{1} --> DynareDataset
% varargin{2} --> DatasetInfo
% varargin{3} --> DynareOptions
% varargin{4} --> Model
% varargin{5} --> EstimatedParameters
% varargin{6} --> BayesInfo
% varargin{1} --> DynareResults
% Copyright (C) 2006-2012 Dynare Team
% Copyright (C) 2006-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -22,7 +30,7 @@ function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,DynareDataset,DynareOptions,Mod
n=size(x,1);
if isempty(h1),
h1=DynareOptions.gradient_epsilon*ones(n,1);
h1=varargin{3}.gradient_epsilon*ones(n,1);
end
@ -31,7 +39,7 @@ if isempty(htol0)
else
htol = htol0;
end
f0=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
f0=feval(func0,x,varargin{:});
xh1=x;
f1=zeros(size(f0,1),n);
@ -45,10 +53,10 @@ while i<n
hcheck=0;
dx=[];
xh1(i)=x(i)+h1(i);
fx = feval(func0,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
fx = feval(func0,xh1,varargin{:});
f1(:,i)=fx;
xh1(i)=x(i)-h1(i);
fx = feval(func0,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
fx = feval(func0,xh1,varargin{:});
f_1(:,i)=fx;
if hcheck && htol<1
htol=min(1,max(min(abs(dx))*2,htol*10));
@ -61,9 +69,9 @@ while i<n
gg(i)=(f1(i)'-f_1(i)')./(2.*h1(i));
hh(i) = 1/max(1.e-9,abs( (f1(i)+f_1(i)-2*f0)./(h1(i)*h1(i)) ));
if gg(i)*(hh(i)*gg(i))/2 > htol
[f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),varargin{:});
ig(i)=1;
fprintf(['Done for param %s = %8.4f\n'],BayesInfo.name{i},x(i))
fprintf(['Done for param %s = %8.4f\n'],varargin{6}.name{i},x(i))
end
xh1=x;
end

50
matlab/mr_hessian.m
View File

@ -1,10 +1,10 @@
function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin)
% [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin)
%
% numerical gradient and Hessian, with 'automatic' check of numerical
% error
%
% adapted from Michel Juillard original rutine hessian.m
% adapted from Michel Juillard original routine hessian.m
%
% func = function handle. The function must give two outputs:
% - the log-likelihood AND the single contributions at times t=1,...,T
@ -22,9 +22,17 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hf
% htol0 = 'precision' of increment of function values for numerical
% derivatives
%
% varargin: other parameters of func
% varargin{1} --> DynareDataset
% varargin{2} --> DatasetInfo
% varargin{3} --> DynareOptions
% varargin{4} --> Model
% varargin{5} --> EstimatedParameters
% varargin{6} --> BayesInfo
% varargin{1} --> DynareResults
% Copyright (C) 2004-2012 Dynare Team
% Copyright (C) 2004-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -45,16 +53,16 @@ persistent h1 htol
n=size(x,1);
if init
gstep_=DynareOptions.gstep;
gstep_=varargin{3}.gstep;
htol = 1.e-4;
h1=DynareOptions.gradient_epsilon*ones(n,1);
h1=varargin{3}.gradient_epsilon*ones(n,1);
return
end
[f0, ff0]=feval(func,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
h2=BayesInfo.ub-BayesInfo.lb;
hmax=BayesInfo.ub-x;
hmax=min(hmax,x-BayesInfo.lb);
[f0, ff0]=feval(func,x,varargin{:});
h2=varargin{6}.ub-varargin{6}.lb;
hmax=varargin{6}.ub-x;
hmax=min(hmax,x-varargin{6}.lb);
if isempty(ff0),
outer_product_gradient=0;
else
@ -83,7 +91,7 @@ while i<n
hcheck=0;
xh1(i)=x(i)+h1(i);
try
[fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fx, ffx]=feval(func,xh1,varargin{:});
catch
fx=1.e8;
end
@ -104,7 +112,7 @@ while i<n
h1(i) = max(h1(i),1.e-10);
xh1(i)=x(i)+h1(i);
try
[fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fx, ffx]=feval(func,xh1,varargin{:});
catch
fx=1.e8;
end
@ -113,21 +121,21 @@ while i<n
h1(i)= htol/abs(dx(it))*h1(i);
xh1(i)=x(i)+h1(i);
try
[fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fx, ffx]=feval(func,xh1,varargin{:});
catch
fx=1.e8;
end
while (fx-f0)==0
h1(i)= h1(i)*2;
xh1(i)=x(i)+h1(i);
[fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fx, ffx]=feval(func,xh1,varargin{:});
ic=1;
end
end
it=it+1;
dx(it)=(fx-f0);
h0(it)=h1(i);
if (h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i))% || (icount==10 && abs(dx(it))>(3*htol)),
if (h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i))
ic=1;
hcheck=1;
end
@ -141,7 +149,7 @@ while i<n
end
end
xh1(i)=x(i)-h1(i);
[fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fx, ffx]=feval(func,xh1,varargin{:});
f_1(:,i)=fx;
if outer_product_gradient,
if any(isnan(ffx)) || isempty(ffx),
@ -181,8 +189,8 @@ if outer_product_gradient,
xh1(j)=x(j)+h_1(j);
xh_1(i)=x(i)-h1(i);
xh_1(j)=x(j)-h_1(j);
temp1 = feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
temp2 = feval(func,xh_1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
temp1 = feval(func,xh1,varargin{:});
temp2 = feval(func,xh_1,varargin{:});
hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
xh1(i)=x(i);
xh1(j)=x(j);
@ -203,7 +211,7 @@ if outer_product_gradient,
end
end
end
gga=ggh.*kron(ones(size(ff1)),2.*h1'); % re-scaled gradient
hh_mat=gga'*gga; % rescaled outer product hessian
hh_mat0=ggh'*ggh; % outer product hessian
@ -222,7 +230,7 @@ if outer_product_gradient,
sd=sqrt(diag(ihh)); %standard errors
sdh=sqrt(1./diag(hh)); %diagonal standard errors
for j=1:length(sd)
sd0(j,1)=min(BayesInfo.p2(j), sd(j)); %prior std
sd0(j,1)=min(varargin{6}.p2(j), sd(j)); %prior std
sd0(j,1)=10^(0.5*(log10(sd0(j,1))+log10(sdh(j,1))));
end
ihh=ihh./(sd*sd').*(sd0*sd0'); %inverse outer product with modified std's
@ -239,7 +247,7 @@ if outer_product_gradient,
if hflag<2
hessian_mat=hh_mat0(:);
end
if any(isnan(hessian_mat))
hh_mat0=eye(length(hh_mat0));
ihh=hh_mat0;

2
matlab/ms-sbvar/ms_sbvar_setup.m
View File

@ -122,7 +122,7 @@ pctindx = [];
% Select the variable to use and rearrange columns if desired
%vlist = [3 1 2];
%options_.ms.vlist = [1 2 3];
options_.ms.vlist = 1:size(options_.varobs,1);
options_.ms.vlist = 1:length(options_.varobs);
vlist1=options_.ms.vlist;
%==========================================================================

2
matlab/ms-sbvar/ms_write_markov_file.m
View File

@ -28,7 +28,7 @@ function ms_write_markov_file(fname, options)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
n_chains = length(options.ms.ms_chain);
nvars = size(options.varobs,1);
nvars = length(options.varobs);
fh = fopen(fname,'w');
%/******************************************************************************/

8
matlab/ms-sbvar/msstart_setup.m
View File

@ -48,16 +48,16 @@ end
%1 CBO output gap -- log(x_t)-log(x_t potential)
%2 GDP deflator -- (P_t/P_{t-1})^4-1.0
%2 FFR/100.
options_.ms.vlist = [1:size(options_.varobs,1)]; % 1: U; 4: PCE inflation.
options_.ms.varlist=cellstr(options_.varobs);
options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,options_.varobs)); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already).
options_.ms.vlist = [1:length(options_.varobs)]; % 1: U; 4: PCE inflation.
options_.ms.varlist=cellstr(options_.varobs');
options_.ms.log_var = sort(varlist_indices(options_.ms.vlistlog,char(options_.varobs))); % subset of "options_.ms.vlist. Variables in log level so that differences are in **monthly** growth, unlike R and U which are in annual percent (divided by 100 already).
options_.ms.percent_var =setdiff(options_.ms.vlist,options_.ms.log_var);
%options_.ms.restriction_fname='ftd_upperchol3v'; %Only used by msstart2.m.
ylab = options_.ms.varlist;
xlab = options_.ms.varlist;
%----------------
nvar = size(options_.varobs,1); % number of endogenous variables
nvar = length(options_.varobs); % number of endogenous variables
nlogeno = length(options_.ms.log_var); % number of endogenous variables in options_.ms.log_var
npereno = length(options_.ms.percent_var); % number of endogenous variables in options_.ms.percent_var
if (nvar~=(nlogeno+npereno))

75
matlab/newrat.m
View File

@ -1,4 +1,4 @@
function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, varargin)
% [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, varargin)
%
% Optimiser with outer product gradient and with sequences of univariate steps
@ -21,9 +21,16 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ft
% with correlation structure as from outer product gradient,
% flagg = 2, full numerical Hessian
%
% varargin = list of parameters for func0
% varargin{1} --> DynareDataset
% varargin{2} --> DatasetInfo
% varargin{3} --> DynareOptions
% varargin{4} --> Model
% varargin{5} --> EstimatedParameters
% varargin{6} --> BayesInfo
% varargin{1} --> DynareResults
% Copyright (C) 2004-2013 Dynare Team
% Copyright (C) 2004-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -52,19 +59,19 @@ ftol=ftol0;
gtol=1.e-3;
htol=htol_base;
htol0=htol_base;
gibbstol=length(BayesInfo.pshape)/50; %25;
gibbstol=length(varargin{6}.pshape)/50; %25;
% func0 = str2func([func2str(func0),'_hh']);
% func0 = func0;
[fval0,gg,hh]=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval0,gg,hh]=feval(func0,x,varargin{:});
fval=fval0;
% initialize mr_gstep and mr_hessian
outer_product_gradient=1;
if isempty(hh)
mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
mr_hessian(1,x,[],[],[],varargin{:});
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,varargin{:});
if isempty(dum),
outer_product_gradient=0;
igg = 1e-4*eye(nx);
@ -111,9 +118,9 @@ while norm(gg)>gtol && check==0 && jit<nit
objective_function_penalty_base = fval0(icount);
disp([' '])
disp(['Iteration ',num2str(icount)])
[fval,x0,fc,retcode] = csminit1(func0,xparam1,fval0(icount),gg,0,H,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval,x0,fc,retcode] = csminit1(func0,xparam1,fval0(icount),gg,0,H,varargin{:});
if igrad
[fval1,x01,fc,retcode1] = csminit1(func0,x0,fval,gg,0,inx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval1,x01,fc,retcode1] = csminit1(func0,x0,fval,gg,0,inx,varargin{:});
if (fval-fval1)>1
disp('Gradient step!!')
else
@ -122,28 +129,26 @@ while norm(gg)>gtol && check==0 && jit<nit
fval=fval1;
x0=x01;
end
% if icount==1 || (icount>1 && (fval0(icount-1)-fval0(icount))>1) || ((fval0(icount)-fval)<1.e-2*(gg'*(H*gg))/2 && igibbs),
if length(find(ig))<nx
ggx=ggx*0;
ggx(find(ig))=gg(find(ig));
if analytic_derivation,
hhx=hh;
else
hhx = reshape(dum,nx,nx);
end
iggx=eye(length(gg));
iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
[fvala,x0,fc,retcode] = csminit1(func0,x0,fval,ggx,0,iggx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
if length(find(ig))<nx
ggx=ggx*0;
ggx(find(ig))=gg(find(ig));
if analytic_derivation,
hhx=hh;
else
hhx = reshape(dum,nx,nx);
end
[fvala, x0, ig] = mr_gstep(h1,x0,func0,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
nig=[nig ig];
disp('Sequence of univariate steps!!')
fval=fvala;
% end
iggx=eye(length(gg));
iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
[fvala,x0,fc,retcode] = csminit1(func0,x0,fval,ggx,0,iggx,varargin{:});
end
[fvala, x0, ig] = mr_gstep(h1,x0,func0,htol,varargin{:});
nig=[nig ig];
disp('Sequence of univariate steps!!')
fval=fvala;
if (fval0(icount)-fval)<ftol && flagit==0
disp('Try diagonal Hessian')
ihh=diag(1./(diag(hhg)));
[fval2,x0,fc,retcode2] = csminit1(func0,x0,fval,gg,0,ihh,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval2,x0,fc,retcode2] = csminit1(func0,x0,fval,gg,0,ihh,varargin{:});
if (fval-fval2)>=ftol
disp('Diagonal Hessian successful')
end
@ -152,7 +157,7 @@ while norm(gg)>gtol && check==0 && jit<nit
if (fval0(icount)-fval)<ftol && flagit==0
disp('Try gradient direction')
ihh0=inx.*1.e-4;
[fval3,x0,fc,retcode3] = csminit1(func0,x0,fval,gg,0,ihh0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fval3,x0,fc,retcode3] = csminit1(func0,x0,fval,gg,0,ihh0,varargin{:});
if (fval-fval3)>=ftol
disp('Gradient direction successful')
end
@ -165,14 +170,14 @@ while norm(gg)>gtol && check==0 && jit<nit
disp('No further improvement is possible!')
check=1;
if analytic_derivation,
[fvalx,gg,hh]=feval(func0,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fvalx,gg,hh]=feval(func0,xparam1,varargin{:});
hhg=hh;
H = inv(hh);
H = inv(hh);
else
if flagit==2
hh=hh0;
elseif flagg>0
[dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,flagg,ftol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,flagg,ftol0,varargin{:});
if flagg==2
hh = reshape(dum,nx,nx);
ee=eig(hh);
@ -208,7 +213,7 @@ while norm(gg)>gtol && check==0 && jit<nit
catch
save m1.mat x fval0 nig
end
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,flagit,htol,varargin{:});
if isempty(dum),
outer_product_gradient=0;
end
@ -237,7 +242,7 @@ while norm(gg)>gtol && check==0 && jit<nit
H = igg;
end
elseif analytic_derivation,
[fvalx,gg,hh]=feval(func0,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
[fvalx,gg,hh]=feval(func0,xparam1,varargin{:});
hhg=hh;
H = inv(hh);
end
@ -276,6 +281,4 @@ if check==1,
disp(['Estimation successful.'])
end
return
return

1
matlab/osr.m
View File

@ -41,7 +41,6 @@ function osr_res = osr(var_list,params,i_var,W)
global M_ options_ oo_
options_.order = 1;
options_ = set_default_option(options_,'simul',0);
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;

5
matlab/particle/solve_model_for_online_filter.m
View File

@ -201,9 +201,10 @@ end
% Get the off-diagonal elements of the covariance matrix for the measurement errors. Test if H is positive definite.
if EstimatedParameters.ncn
corrn_observable_correspondence = EstimatedParameters.corrn_observable_correspondence;
for i=1:EstimatedParameters.ncn
k1 = DynareOptions.lgyidx2varobs(EstimatedParameters.corrn(i,1));
k2 = DynareOptions.lgyidx2varobs(EstimatedParameters.corrn(i,2));
k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2);
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
H(k2,k1) = H(k1,k2);
end

25
matlab/prior_posterior_statistics.m
View File

@ -1,4 +1,4 @@
function prior_posterior_statistics(type,dataset)
function prior_posterior_statistics(type,dataset,dataset_info)
% function prior_posterior_statistics(type,dataset)
% Computes Monte Carlo filter smoother and forecasts
@ -40,11 +40,11 @@ global options_ estim_params_ oo_ M_ bayestopt_
localVars=[];
Y = dataset.data;
gend = dataset.info.ntobs;
data_index = dataset.missing.aindex;
missing_value = dataset.missing.state;
bayestopt_.mean_varobs = dataset.descriptive.mean';
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';
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
@ -58,7 +58,7 @@ naK = length(options_.filter_step_ahead);
MaxNumberOfBytes=options_.MaxNumberOfBytes;
endo_nbr=M_.endo_nbr;
exo_nbr=M_.exo_nbr;
nvobs = size(options_.varobs,1);
nvobs = length(options_.varobs);
iendo = 1:endo_nbr;
horizon = options_.forecast;
filtered_vars = options_.filtered_vars;
@ -97,7 +97,7 @@ MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
if naK
MAX_naK = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ...
MAX_naK = min(B,ceil(MaxNumberOfBytes/(length(options_.varobs)* ...
length(options_.filter_step_ahead)*gend)/8));
end
@ -106,7 +106,7 @@ if horizon
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)));
@ -232,7 +232,7 @@ else
'options_', options_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_);
'oo_', oo_);
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
@ -240,11 +240,10 @@ else
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
end
[fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'prior_posterior_statistics_core', localVars,globalVars, options_.parallel_info);
end
ifil = fout(end).ifil;
% ??????????
stock_gend=gend;
stock_data=Y;
save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data');
@ -301,4 +300,4 @@ if ~isnumeric(options_.parallel),
if leaveSlaveOpen == 0,
closeSlave(options_.parallel,options_.parallel_info.RemoteTmpFolder),
end
end
end

24
matlab/prior_posterior_statistics_core.m
View File

@ -133,14 +133,14 @@ end
if run_smoother
stock_smooth=NaN(endo_nbr,gend,MAX_nsmoo);
stock_update=NaN(endo_nbr,gend,MAX_nsmoo);
stock_innov=NaN(M_.exo_nbr,gend,MAX_ninno);
stock_innov=NaN(M_.exo_nbr,gend,MAX_ninno);
if horizon
stock_forcst_mean= NaN(endo_nbr,horizon+maxlag,MAX_nforc1);
stock_forcst_point = NaN(endo_nbr,horizon+maxlag,MAX_nforc2);
end
end
if nvn
stock_error = NaN(size(varobs,1),gend,MAX_nerro);
stock_error = NaN(length(varobs),gend,MAX_nerro);
end
if naK
stock_filter_step_ahead =NaN(length(options_.filter_step_ahead),endo_nbr,gend+max(options_.filter_step_ahead),MAX_naK);
@ -150,10 +150,6 @@ stock_logpo = NaN(MAX_nruns,1);
stock_ys = NaN(MAX_nruns,endo_nbr);
for b=fpar:B
% [deep, logpo] = GetOneDraw(type);
% set_all_parameters(deep);
% dr = resol(oo_.steady_state,0);
if strcmpi(type,'prior')
[deep, logpo] = GetOneDraw(type);
@ -304,22 +300,7 @@ for b=fpar:B
end
irun(7) = 1;
end
% if moments_varendo && (irun(8) > MAX_momentsno || b == B)
% stock = stock_moments(1:irun(8)-1);
% ifil(8) = ifil(8) + 1;
% save([DirectoryName '/' M_.fname '_moments' int2str(ifil(8)) '.mat'],'stock');
% if RemoteFlag==1,
% OutputFileName_moments = [OutputFileName_moments; {[DirectoryName filesep], [M_.fname '_moments' int2str(ifil(8)) '.mat']}];
% end
% irun(8) = 1;
% end
% DirectoryName=TempPath;
dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
end
myoutput.ifil=ifil;
@ -332,7 +313,6 @@ if RemoteFlag==1,
OutputFileName_param;
OutputFileName_forc_mean;
OutputFileName_forc_point];
% OutputFileName_moments];
end
dyn_waitbar_close(h);

5
matlab/random_walk_metropolis_hastings.m
View File

@ -1,4 +1,4 @@
function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
function random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_)
%function record=random_walk_metropolis_hastings(TargetFun,ProposalFun,xparam1,vv,mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_)
% Random walk Metropolis-Hastings algorithm.
%
@ -60,7 +60,7 @@ objective_function_penalty_base = Inf;
% Initialization of the random walk metropolis-hastings chains.
[ ix2, ilogpo2, ModelName, MetropolisFolder, fblck, fline, npar, nblck, nruns, NewFile, MAX_nruns, d ] = ...
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
metropolis_hastings_initialization(TargetFun, xparam1, vv, mh_bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
InitSizeArray = min([repmat(MAX_nruns,nblck,1) fline+nruns-1],[],2);
@ -101,6 +101,7 @@ localVars = struct('TargetFun', TargetFun, ...
'InitSizeArray',InitSizeArray, ...
'record', record, ...
'dataset_', dataset_, ...
'dataset_info', dataset_info, ...
'options_', options_, ...
'M_',M_, ...
'bayestopt_', bayestopt_, ...

3
matlab/random_walk_metropolis_hastings_core.m
View File

@ -88,6 +88,7 @@ d=myinputs.d;
InitSizeArray=myinputs.InitSizeArray;
record=myinputs.record;
dataset_ = myinputs.dataset_;
dataset_info = myinputs.dataset_info;
bayestopt_ = myinputs.bayestopt_;
estim_params_ = myinputs.estim_params_;
options_ = myinputs.options_;
@ -164,7 +165,7 @@ for b = fblck:nblck,
par = feval(ProposalFun, ix2(b,:), proposal_covariance_Cholesky_decomposition, n);
if all( par(:) > mh_bounds(:,1) ) && all( par(:) < mh_bounds(:,2) )
try
logpost = - feval(TargetFun, par(:),dataset_,options_,M_,estim_params_,bayestopt_,oo_);
logpost = - feval(TargetFun, par(:),dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
catch
logpost = -inf;
end

11
matlab/read_variables.m
View File

@ -42,7 +42,7 @@ if ~isempty(directory)
end
dyn_size_01 = size(dyn_data_01,1);
var_size_01 = size(var_names_01,1);
var_size_01 = length(var_names_01);
% Auto-detect extension if not provided
if isempty(extension)
@ -71,7 +71,7 @@ switch (extension)
case '.m'
eval(basename);
for dyn_i_01=1:var_size_01
dyn_tmp_01 = eval(var_names_01(dyn_i_01,:));
dyn_tmp_01 = eval(var_names_01{dyn_i_01});
if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0
cd(old_pwd)
error('data size is too large')
@ -81,7 +81,7 @@ switch (extension)
case '.mat'
s = load(basename);
for dyn_i_01=1:var_size_01
dyn_tmp_01 = s.(deblank(var_names_01(dyn_i_01,:)));
dyn_tmp_01 = s.(var_names_01{dyn_i_01});
if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0
cd(old_pwd)
error('data size is too large')
@ -106,9 +106,8 @@ switch (extension)
end
case '.csv'
[freq,init,data,varlist] = load_csv_file_data(fullname);
%var_names_01 = deblank(var_names_01);
for dyn_i_01=1:var_size_01
iv = strmatch(strtrim(var_names_01(dyn_i_01,:)),varlist,'exact');
iv = strmatch(var_names_01{dyn_i_01},varlist,'exact');
if ~isempty(iv)
dyn_tmp_01 = [data(:,iv)]';
if length(dyn_tmp_01) > dyn_size_01 && dyn_size_01 > 0
@ -118,7 +117,7 @@ switch (extension)
dyn_data_01(:,dyn_i_01) = dyn_tmp_01;
else
cd(old_pwd)
error([strtrim(var_names_01(dyn_i_01,:)) ' not found in ' fullname])
error([var_names_01{dyn_i_01} ' not found in ' fullname])
end
end
otherwise

11
matlab/set_default_option.m
View File

@ -33,6 +33,17 @@ function options=set_default_option(options,field,default)
if ~isfield(options,field)
options.(field) = default;
return
end
if isempty(options.(field))
options.(field) = default;
return
end
if isnan(options.(field))
options.(field) = default;
return
end
% 06/07/03 MJ added ; to eval expression

12
matlab/set_prior.m
View File

@ -76,12 +76,12 @@ end
if nvn
estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable
if isequal(M_.H,0) %if no previously set measurement error, initialize H
nvarobs = size(options_.varobs,1);
nvarobs = length(options_.varobs);
M_.H = zeros(nvarobs,nvarobs);
M_.Correlation_matrix_ME = eye(nvarobs);
end
for i=1:nvn
obsi_ = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact');
obsi_ = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),options_.varobs,'exact');
if isempty(obsi_)
error(['The variable ' deblank(M_.endo_names(estim_params_.var_endo(i,1),:)) ' has to be declared as observable since you assume a measurement error on it.'])
end
@ -96,7 +96,7 @@ if nvn
bayestopt_.p3 = [ bayestopt_.p3; estim_params_.var_endo(:,8)]; %take generalized distribution into account
bayestopt_.p4 = [ bayestopt_.p4; estim_params_.var_endo(:,9)]; %take generalized distribution into account
bayestopt_.jscale = [ bayestopt_.jscale; estim_params_.var_endo(:,10)];
bayestopt_.name = [ bayestopt_.name; cellstr(options_.varobs(estim_params_.nvn_observable_correspondence,:))];
bayestopt_.name = [ bayestopt_.name; transpose(options_.varobs(estim_params_.nvn_observable_correspondence))];
end
if ncx
xparam1 = [xparam1; estim_params_.corrx(:,3)];
@ -115,7 +115,7 @@ end
if ncn
estim_params_.corrn_observable_correspondence=NaN(ncn,2);
if isequal(M_.H,0)
nvarobs = size(options_.varobs,1);
nvarobs = length(options_.varobs);
M_.H = zeros(nvarobs,nvarobs);
M_.Correlation_matrix_ME = eye(nvarobs);
end
@ -134,8 +134,8 @@ if ncn
for i=1:ncn
k1 = estim_params_.corrn(i,1);
k2 = estim_params_.corrn(i,2);
obsi1 = strmatch(deblank(M_.endo_names(k1,:)),deblank(options_.varobs),'exact'); %find correspondence to varobs to construct H in set_all_paramters
obsi2 = strmatch(deblank(M_.endo_names(k2,:)),deblank(options_.varobs),'exact');
obsi1 = strmatch(deblank(M_.endo_names(k1,:)),options_.varobs,'exact'); %find correspondence to varobs to construct H in set_all_paramters
obsi2 = strmatch(deblank(M_.endo_names(k2,:)),options_.varobs,'exact');
estim_params_.corrn_observable_correspondence(i,:)=[obsi1,obsi2]; %save correspondence
end
end

19
matlab/simplex_optimization_routine.m
View File

@ -11,8 +11,15 @@ function [x,fval,exitflag] = simplex_optimization_routine(objective_function,x,o
% o objective_function [string] Name of the objective function to be minimized.
% o x [double] n*1 vector, starting guess of the optimization routine.
% o options [structure] Options of this implementation of the simplex algorithm.
% o varargin [cell of structures] Structures to be passed to the objective function: dataset_,
% options_, M_, estim_params_, bayestopt_, and oo_.
% o varargin [cell of structures] Structures to be passed to the objective function.
%
% varargin{1} --> DynareDataset
% varargin{2} --> DatasetInfo
% varargin{3} --> DynareOptions
% varargin{4} --> Model
% varargin{5} --> EstimatedParameters
% varargin{6} --> BayesInfo
% varargin{1} --> DynareResults
%
% OUTPUTS
% o x [double] n*1 vector, estimate of the optimal inputs.
@ -187,7 +194,7 @@ else
disp(['Current parameter values: '])
fprintf(1,'%s: \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \n','Names','Best point', 'Worst point', 'Mean values', 'Min values', 'Max values');
for i=1:number_of_variables
fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{5}.name{i},v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2));
fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{6}.name{i},v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2));
end
skipline()
end
@ -399,7 +406,7 @@ while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex
disp(['Current parameter values: '])
fprintf(1,'%s: \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \n','Names','Best point', 'Worst point', 'Mean values', 'Min values', 'Max values');
for i=1:number_of_variables
fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{5}.name{i}, v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2));
fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{6}.name{i}, v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2));
end
skipline()
end
@ -425,7 +432,7 @@ while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex
disp(['values for the control variables. '])
disp(['New value of delta (size of the new simplex) is: '])
for i=1:number_of_variables
fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{5}.name{i}, delta(i));
fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{6}.name{i}, delta(i));
end
end
% Reset counters
@ -472,7 +479,7 @@ while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex
disp(['values for the control variables. '])
disp(['New value of delta (size of the new simplex) is: '])
for i=1:number_of_variables
fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{5}.name{i}, delta(i));
fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{6}.name{i}, delta(i));
end
end
% Reset counters

2
matlab/simulated_moment_uncertainty.m
View File

@ -30,7 +30,7 @@ for j=1:replic;
options_.noprint = 1;
options_.order = 1;
options_.periods = periods;
info = stoch_simul(options_.varobs);
info = stoch_simul(char(options_.varobs));
dum=[oo_.mean; dyn_vech(oo_.var)];
sd = sqrt(diag(oo_.var));
for i=1:options_.ar;

2
matlab/stoch_simul.m
View File

@ -115,7 +115,7 @@ if ~options_.noprint
skipline()
if isfield(options_,'varobs')&& ~isempty(options_.varobs)
PCL_varobs=options_.varobs;
PCL_varobs=char(options_.varobs);
disp('OBSERVED VARIABLES')
else
PCL_varobs=M_.endo_names;

2
matlab/subset.m
View File

@ -20,8 +20,6 @@ function jndx = subset()
global options_ estim_params_ M_
ExcludedParamNames = options_.ExcludedParams;
VarObs = options_.varobs;
VarExo = M_.exo_names;
info = options_.ParamSubSet;
nvx = estim_params_.nvx;

59
matlab/utilities/dataset/compute_corr.m
View File

@ -1,59 +0,0 @@
function dataset_ = compute_corr(dataset_)
% Computes the correlation matrix of the sample (possibly with missing observations).
%@info:
%! @deftypefn {Function File} {@var{dataset_} =} compute_corr(@var{dataset_})
%! @anchor{compute_corr}
%! This function computes covariance matrix of the sample (possibly with missing observations).
%!
%! @strong{Inputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @end table
%!
%! @strong{Outputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @end table
%!
%! @strong{This function is called by:}
%! @ref{descriptive_statistics}.
%!
%! @strong{This function calls:}
%! @ref{ndim}, @ref{compute_cova}.
%!
%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.corr} is a
%! @tex{n\times n} vector (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}).
%!
%! @strong{Remark 2.} If @code{dataset_.descriptive.cova} does not exist, the covariance matrix is computed prior to the
%! computation of the correlation matrix.
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr
if ~isfield(dataset_.descriptive,'cova')
dataset_ = compute_cova(dataset_);
end
normalization_matrix = diag(1./sqrt(diag(dataset_.descriptive.cova)));
dataset_.descriptive.corr = normalization_matrix*dataset_.descriptive.cova*normalization_matrix;

67
matlab/utilities/dataset/compute_cova.m
View File

@ -1,67 +0,0 @@
function dataset_ = compute_cova(dataset_)
% Computes the covariance matrix of the sample (possibly with missing observations).
%@info:
%! @deftypefn {Function File} {@var{dataset_} =} compute_corr(@var{dataset_})
%! @anchor{compute_corr}
%! This function computes covariance matrix of the sample (possibly with missing observations).
%!
%! @strong{Inputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @end table
%!
%! @strong{Outputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @end table
%!
%! @strong{This function is called by:}
%! @ref{descriptive_statistics}.
%!
%! @strong{This function calls:}
%! @ref{ndim}, @ref{demean}, @ref{nandemean}.
%!
%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.cova} is a
%! @tex{n\times n} vector (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}).
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr
dataset_.descriptive.cova = zeros(dataset_.nvobs);
data = transpose(dataset_.data);
for i=1:dataset_.info.nvobs
for j=i:dataset_.info.nvobs
if dataset_.missing.state
dataset_.descriptive.cova(i,j) = nanmean(nandemean(data(:,i)).*nandemean(data(:,j)));
else
dataset_.descriptive.cova(i,j) = mean(demean(data(:,i)).*demean(data(:,j)));
end
if j>i
dataset_.descriptive.cova(j,i) = dataset_.descriptive.cova(i,j);
end
end
end

14
matlab/utilities/dataset/describe_missing_data.m
View File

@ -26,9 +26,9 @@ function [i,n,s,j] = describe_missing_data(data)
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2008-2012 Dynare Team
%
% Copyright (C) 2008-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
@ -85,7 +85,7 @@ end
%@test:1
%$ % Define a data set.
%$ % Define a data set.
%$ A = [ 1 1 ; ...
%$ 1 NaN ; ...
%$ NaN 1 ; ...
@ -99,7 +99,7 @@ end
%$ 1 1 ];
%$
%$ % Define expected results.
%$ eB = cell(1,11);
%$ eB = cell(1,11);
%$ eB(1) = { transpose(1:2) };
%$ eB(2) = { 1 };
%$ eB(3) = { 2 };
@ -107,7 +107,7 @@ end
%$ eB(5) = { [] };
%$ eB(6) = { 1 };
%$ eB(7) = { 1 };
%$ eB(8) = { transpose(1:2) };
%$ eB(8) = { transpose(1:2) };
%$ eB(9) = { transpose(1:2) };
%$ eB(10) = { transpose(1:2) };
%$ eB(11) = { transpose(1:2) };
@ -117,7 +117,7 @@ end
%$ eE(1) = { [1; 2; 4; transpose(6:11)] };
%$ eE(2) = { [1; 3; 4; transpose(8:11)] };
%$
%$ % Call the tested routine.
%$ % Call the tested routine.
%$ [B,C,D,E] = describe_missing_data(transpose(A));
%$
%$ % Check the results.

14
matlab/utilities/dataset/initialize_dataset.m
View File

@ -1,4 +1,4 @@
function dataset_ = initialize_dataset(datafile,varobs,first,nobs,transformation,prefilter,xls)
function dataset_ = initialize_dataset(datafile,varobs,first,nobs,logged_data_flag,prefilter,xls)
% Initializes a structure describing the dataset.
% Copyright (C) 2011-2013 Dynare Team
@ -18,8 +18,6 @@ function dataset_ = initialize_dataset(datafile,varobs,first,nobs,transformation
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr
if isempty(datafile)
error('Estimation::initialize_dataset: You have to declare a dataset file!')
end
@ -43,7 +41,7 @@ end
% Fill the dataset structure
dataset_.info.ntobs = nobs;
dataset_.info.nvobs = rows(varobs);
dataset_.info.nvobs = length(varobs);
dataset_.info.varobs = varobs;
% Test the number of variables in the database.
@ -61,9 +59,9 @@ if size(rawdata,1)~=dataset_.info.ntobs
end
rawdata = rawdata(first:(first+dataset_.info.ntobs-1),:);
% Take the log (or anything else) of the variables if needed
if isempty(transformation)
dataset_.rawdata = rawdata;
% Take the log of the variables if needed
if logged_data_flag
dataset_.rawdata = log(rawdata);
else
if isequal(transformation,@log)
if ~isreal(rawdata)
@ -101,4 +99,4 @@ if prefilter == 1
dataset_.data = transpose(bsxfun(@minus,dataset_.rawdata,dataset_.descriptive.mean));
else
dataset_.data = transpose(dataset_.rawdata);
end
end

241
matlab/utilities/dataset/makedataset.m Normal file
View File

@ -0,0 +1,241 @@
function [DynareDataset, DatasetInfo] = makedataset(DynareOptions, initialconditions, gsa_flag)
% Initialize a dataset as a dseries object.
%
%
% INPUTS
% ======
%
% DynareOptions [struct] Structure of options built by Dynare's preprocessor.
%
%
% OUTPUTS
% =======
%
% DynareDataset [dseries] The dataset.
% DatasetInfo [struct] Various informations about the dataset (descriptive statistics and missing observations).
%
% EXAMPLE
% =======
%
% [dataset_, dataset_info] = makedataset(options_) ;
%
%
% See also dynare_estimation_init
if nargin<3
gsa_flag = 0;
end
if nargin<2 || isempty(initialconditions)
% If a the sample is to be used for the estimation of a VAR or DSGE-VAR model
% the second argument must be a strictly positive integer (the number of lags).
initialconditions = 0;
end
if isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file) && isempty(DynareOptions.dataset.series)
if gsa_flag
DynareDataset = dseries();
DatasetInfo = [];
return
else
error('makedataset: datafile option is missing!')
end
end
if isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.file)
datafile = DynareOptions.dataset.file;
newdatainterface = 1;
elseif isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.series)
try
dseriesobjectforuserdataset = evalin('base', DynareOptions.dataset.series);
catch
error(sprintf('makedataset: %s is unknown!', DynareOptions.dataset.series))
end
if ~isdseries(dseriesobjectforuserdataset)
error(sprintf('makedataset: %s has to be a dseries object!', DynareOptions.dataset.series))
end
datafile = [];
newdatainterface = 1;
elseif ~isempty(DynareOptions.datafile) && isempty(DynareOptions.dataset.file)
datafile = DynareOptions.datafile;
newdatainterface = 0;
elseif isempty(DynareOptions.datafile) && ~isempty(DynareOptions.dataset.file)
error('You cannot use simultaneously the data command and the datafile option (in the estimation command)!')
else
error('You have to specify the datafile!')
end
% Check extension.
if ~isempty(datafile)
allowed_extensions = {'m','mat','csv','xls','xlsx'};
datafile_extension = get_file_extension(datafile);
if isempty(datafile_extension)
available_extensions = {}; j = 1;
for i=1:length(allowed_extensions)
if exist([datafile '.' allowed_extensions{i}])
available_extensions(j) = {allowed_extensions{i}};
j = j+1;
end
end
if isempty(available_extensions)
error(['I can''t find a datafile (with allowed extension)!'])
end
if length(available_extensions)>1
error(sprintf(['You did not specify an extension for the datafile, but more than one candidate ' ...
'are available in the designed folder!\nPlease, add an extension to the datafile ' ...
'(m, mat, csv, xls or xlsx are legal extensions).']));
end
datafile = [datafile '.' available_extensions{1}];
end
end
% Load the data in a dseries object.
if ~isempty(datafile)
DynareDataset = dseries(datafile);
else
DynareDataset = dseriesobjectforuserdataset;
clear('dseriesobjectforuserdataset');
end
% Select a subset of the variables.
DynareDataset = DynareDataset{DynareOptions.varobs{:}};
% Apply log function if needed.
if DynareOptions.loglinear && ~DynareOptions.logdata
DynareDataset = DynareDataset.log();
end
% Test if an initial period (different from its default value) is explicitely defined in the datafile.
if isequal(DynareDataset.init, dates(1,1))
dataset_default_initial_period = 1;
else
dataset_default_initial_period = 0;
end
% Test if an initial period (different from its default value) is explicitely defined in the mod file with the set_time command.
if isequal(DynareOptions.initial_period, dates(1,1))
set_time_default_initial_period = 1;
else
set_time_default_initial_period = 0;
end
if ~set_time_default_initial_period && dataset_default_initial_period
% Overwrite the initial period in dataset (it was set to default).
% Note that the updates of freq and time members are auto-magically
% done by dseries::subsasgn overloaded method.
DynareDataset.init = DynareOptions.initial_period;
end
if set_time_default_initial_period && ~dataset_default_initial_period
% Overwrite the global initial period defined by set_time (it was set to default).
DynareOptions.initial_period = DynareDataset.init;
end
if ~set_time_default_initial_period && ~dataset_default_initial_period
% Check if dataset.init and options_.initial_period are identical.
if DynareOptions.initial_period>DynareDataset.init
%if ~isequal(DynareOptions.initial_period, DynareDataset.init)
error('The date as defined by the set_time command is not consistent with the initial period in the database!')
end
end
% Set firstobs, lastobs and nobs
if newdatainterface
if isempty(DynareOptions.dataset.firstobs)
% first_obs option was not used in the data command.
firstobs = DynareDataset.init;
else
firstobs = DynareOptions.dataset.firstobs;
end
if isnan(DynareOptions.dataset.nobs)
% nobs option was not used in the data command.
if isempty(DynareOptions.dataset.lastobs)
% last_obs option was not used in the data command.
nobs = DynareDataset.nobs;
lastobs = DynareDataset.dates(end);
else
lastobs = DynareOptions.dataset.lastobs;
nobs = lastobs-firstobs+1;
end
else
nobs = DynareOptions.dataset.nobs;
if isempty(DynareOptions.dataset.lastobs)
% last_obs option was not used in the data command.
lastobs = firstobs+(nobs-1);
else
% last_obs and nobs were used in the data command. Check that they are consistent (with firstobs).
if ~isequal(lastobs,firstobs+(nobs-1))
error(sprintf('Options last_obs (%s), first_obs (%s) and nobs (%s) are not consistent!',char(lastobs),char(firstobs),num2str(nobs)));
end
end
end
else
if isnan(DynareOptions.first_obs)
firstobs = DynareDataset.init;
else
firstobs = DynareDataset.dates(DynareOptions.first_obs);
end
if isnan(DynareOptions.nobs)
lastobs = DynareDataset.dates(end);
nobs = lastobs-firstobs+1;
else
nobs = DynareOptions.nobs;
lastobs = firstobs+(nobs-1);
end
end
% Add initial conditions if needed
FIRSTOBS = firstobs-initialconditions;
% Check that firstobs belongs to DynareDataset.dates
if firstobs<DynareDataset.init
error(sprintf('first_obs (%s) cannot be less than the first date in the dataset (%s)!',char(firstobs),char(DynareDataset.init)))
end
% Check that FIRSTOBS belongs to DynareDataset.dates
if initialconditions && FIRSTOBS<DynareDataset.init
error(sprintf('first_obs (%s) - %i cannot be less than the first date in the dataset (%s)!\nReduce the number of lags in the VAR model or increase the value of first_obs.', char(firstobs), initialconditions, char(DynareDataset.init)));
end
% Check that lastobs belongs to DynareDataset.dates...
if newdatainterface
if lastobs>DynareDataset.dates(end)
error(sprintf('last_obs (%s) cannot be greater than the last date in the dataset (%s)!',char(lastobs),char(DynareDataset.dates(end))))
end
else
% ... or check that nobs is smaller than the number of observations in DynareDataset.
if nobs>DynareDataset.nobs
error(sprintf('nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(nobs)))
end
end
% Select a subsample.
DynareDataset = DynareDataset(FIRSTOBS:lastobs);
% Initialize DatasetInfo structure.
DatasetInfo = struct('missing', struct('state', NaN, 'aindex', [], 'vindex', [], 'number_of_observations', NaN, 'no_more_missing_observations', NaN), ...
'descriptive', struct('mean', [], 'covariance', [], 'correlation', [], 'autocovariance', []));
% Fill DatasetInfo.missing if some observations are missing
DatasetInfo.missing.state = isanynan(DynareDataset.data);
if DatasetInfo.missing.state
[DatasetInfo.missing.aindex, DatasetInfo.missing.number_of_observations, DatasetInfo.missing.no_more_missing_observations, DatasetInfo.missing.vindex] = ...
describe_missing_data(DynareDataset.data);
else
DatasetInfo.missing.aindex = num2cell(transpose(repmat(1:DynareDataset.vobs,DynareDataset.nobs,1)),1);
DatasetInfo.missing.no_more_missing_observations = 1;
end
% Compute the empirical mean of the observed variables.
DatasetInfo.descriptive.mean = nanmean(DynareDataset.data);
% Compute the empirical covariance matrix of the observed variables.
DatasetInfo.descriptive.covariance = nancovariance(DynareDataset.data);
% Compute the empirical correlation matrix of the observed variables.
normalization_matrix = diag(1./sqrt(diag(DatasetInfo.descriptive.covariance)));
DatasetInfo.descriptive.correlation = normalization_matrix*DatasetInfo.descriptive.covariance*normalization_matrix;
% Compute autocorrelation function.
DatasetInfo.descriptive.autocovariance = nanautocovariance(DynareDataset.data, DynareOptions.ar);

63
matlab/utilities/dataset/compute_acov.m → matlab/utilities/dataset/nanautocovariance.m
View File

@ -1,4 +1,4 @@
function dataset_ = compute_acov(dataset_)
function autocov = nanautocovariance(data,order)
% Computes the (multivariate) auto-covariance function of the sample (possibly with missing observations).
%@info:
@ -8,36 +8,36 @@ function dataset_ = compute_acov(dataset_)
%!
%! @strong{Inputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @item nlag
%! Integer scalar. The maximum number of lags of the autocovariance function.
%! @item data
%! T*N array of real numbers.
%! @item order
%! Integer scalar. The maximum number of lags of the autocovariance function.
%! @end table
%!
%! @strong{Outputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @item autocov
%! A N*N*order array of real numbers.
%! @end table
%!
%! @strong{This function is called by:}
%!
%! @strong{This function is called by:}
%! @ref{descriptive_statistics}.
%!
%!
%! @strong{This function calls:}
%! @ref{ndim}, @ref{compute_cova}, @ref{demean}, @ref{nandemean}.
%!
%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.acov} is a
%! @tex{n\times n\times p} array (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs},
%! @ref{ndim}, @ref{nancovariance}, @ref{demean}, @ref{nandemean}.
%!
%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.acov} is a
%! @tex{n\times n\times p} array (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs},
%! and @tex{n} is the maximum number of lags given by the second input @code{nlag}).
%!
%! @strong{Remark 2.} If @code{dataset_.descriptive.cova} does not exist, the covariance matrix is computed prior to the
%!
%! @strong{Remark 2.} If @code{dataset_.descriptive.cova} does not exist, the covariance matrix is computed prior to the
%! computation of the auto-covariance function.
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2012 Dynare Team
%
% Copyright (C) 2011-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
@ -53,22 +53,23 @@ function dataset_ = compute_acov(dataset_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr
n = size(data,2);
missing = isanynan(data);
if ~isfield(dataset_.descriptive,'cova')
dataset_ = compute_cova(dataset_);
end
dataset_.descriptive.acov = zeros(dataset_.nvobs,dataset_.nvobs,nlag);
autocov = zeros(n, n, order);
data = transpose(dataset_.data);
for lag=1:nlag
for i=1:dataset_.info.nvobs
for j=1:dataset_.info.nvobs
if dataset_.missing.state
dataset_.descriptive.acov(i,j,lag) = nanmean(nandemean(data(lag+1:end,i)).*nandemean(data(1:end-lag,j)));
for lag=1:order
if missing
data = nandemean(data);
else
data = demean(data);
end
for i=1:n
for j=1:n
if missing
autocov(i,j,lag) = nanmean(data((lag+1):end,i).*data(1:end-lag,j));
else
dataset_.descriptive.acov(i,j,lag) = mean(demean(data(lag+1:end,i)).*demean(data(1:end-lag,j)));
autocov(i,j,lag) = mean(data((lag+1):end,i).*data(1:end-lag,j));
end
end
end

101
matlab/utilities/dataset/nancovariance.m Normal file
View File

@ -0,0 +1,101 @@
function CovarianceMatrix = nancovariance(data)
% Computes the covariance matrix of a sample (possibly with missing observations).
%@info:
%! @deftypefn {Function File} {@var{CovarianceMatrix} =} compute_corr(@var{data})
%! @anchor{compute_cova}
%! This function computes covariance matrix of a sample defined by a dseries object (possibly with missing observations).
%!
%! @strong{Inputs}
%! @table @var
%! @item data
%! a T*N array of real numbers.
%! @end table
%!
%! @strong{Outputs}
%! @table @var
%! @item CovarianceMatrix
%! Array of real numbers.
%! @end table
%!
%! @strong{This function is called by:}
%! @ref{descriptive_statistics}.
%!
%! @strong{This function calls:}
%! @ref{ndim}, @ref{demean}, @ref{nandemean}.
%!
%! @strong{Remark 1.} On exit, a new field is appended to the structure: @code{dataset_.descriptive.cova} is a
%! @tex{n\times n} vector (where @tex{n} is the number of observed variables as defined by @code{dataset_.info.nvobs}).
%!
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Initialize the output.
CovarianceMatrix = zeros(size(data,2));
if isanynan(data)
data = bsxfun(@minus,data,nanmean(data));
for i=1:size(data,2)
for j=i:size(data,2)
CovarianceMatrix(i,j) = nanmean(data(:,i).*data(:,j));
if j>i
CovarianceMatrix(j,i) = CovarianceMatrix(i,j);
end
end
end
else
data = bsxfun(@minus,data,mean(data));
CovarianceMatrix = (transpose(data)*data)/size(data,1);
end
%@test:1
%$
%$ % Define a dataset.
%$ data1 = randn(10000000,2);
%$
%$ % Same dataset with missing observations.
%$ data2 = data1;
%$ data2(45,1) = NaN;
%$ data2(57,2) = NaN;
%$ data2(367,:) = NaN(1,2);
%$
%$ t = zeros(2,1);
%$
%$ % Call the tested routine.
%$ try
%$ c1 = nancovariance(data1);
%$ t(1) = 1;
%$ catch
%$ t(1) = 0;
%$ end
%$ try
%$ c2 = nancovariance(data2);
%$ t(2) = 1;
%$ catch
%$ t(2) = 0;
%$ end
%$
%$ if t(1) && t(2)
%$ t(3) = max(max(abs(c1-c2)))<1e-4;
%$ end
%$
%$ % Check the results.
%$ T = all(t);
%@eof:1

25
matlab/utilities/dataset/nanmoments.m Normal file
View File

@ -0,0 +1,25 @@
function m = nanmoments(data, n)
% Compute centered marginal moments of order n (possibly with missing observations).
% Copyright (C) 2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if isanynan(data)
m = transpose(nanmean(bsxfun(@power,nandemean(data),n)));
else
m = transpose(mean(bsxfun(@power,demean(data),n)));
end

26
matlab/utilities/dataset/compute_stdv.m → matlab/utilities/dataset/nanvariance.m
View File

@ -1,21 +1,21 @@
function dataset_ = compute_stdv(dataset_)
function variances = nanvariance(data)
% Compute the standard deviation for each observed variable (possibly with missing observations).
%@info:
%! @deftypefn {Function File} {@var{dataset_} =} compute_stdv(@var{dataset_})
%! @anchor{compute_stdv}
%! This function computes the standard deviation of the observed variables (possibly with missing observations).
%! @deftypefn {Function File} {@var{variances} =} nanvariance(@var{data})
%! @anchor{nanvariance}
%! This function computes the variances of the observed variables (possibly with missing observations).
%!
%! @strong{Inputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @item datas
%! A T*N array of real numbers.
%! @end table
%!
%! @strong{Outputs}
%! @table @var
%! @item dataset_
%! Dynare structure describing the dataset, built by @ref{initialize_dataset}
%! @item variances
%! A N*1 vector of real numbers
%! @end table
%!
%! @strong{This function is called by:}
@ -30,7 +30,7 @@ function dataset_ = compute_stdv(dataset_)
%! @end deftypefn
%@eod:
% Copyright (C) 2011-2012 Dynare Team
% Copyright (C) 2011-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -47,10 +47,8 @@ function dataset_ = compute_stdv(dataset_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr
if dataset_.missing.state
dataset_.descriptive.stdv = sqrt(nanmean(bsxfun(@power,nandemean(transpose(dataset_.data)),2)));
if isanynan(data)
variances = transpose(nanmean(bsxfun(@power,nandemean(data),2)));
else
dataset_.descriptive.stdv = sqrt(mean(bsxfun(@power,demean(transpose(dataset_.data)),2)));
variances = transpose(mean(bsxfun(@power,demean(data),2)));
end

13
matlab/utilities/dates/isdate.m
View File

@ -25,15 +25,16 @@ function b = isdate(str) % --*-- Unitary tests --*--
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%if length(str)>1
b = isquaterly(str) || isyearly(str) || ismonthly(str) || isweekly(str);
%else
%b = 0;
%end
if isnumeric(str) && isscalar(str)
b = 1;
return
end
b = isstringdate(str);
%@test:1
%$
%$ date_1 = '1950M2';
%$ date_1 = 1950;
%$ date_2 = '1950m2';
%$ date_3 = '-1950m2';
%$ date_4 = '1950m52';

13
matlab/utilities/dates/isfreq.m
View File

@ -25,4 +25,15 @@ function B = isfreq(A)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
B = (ischar(A) && ismember(upper(A),{'Y','A','Q','M','W'})) || ismember(A,[1 4 12 52]);
B = 0;
if ischar(A)
if isequal(length(A),1) && ismember(upper(A),{'Y','A','Q','M','W'})
B = 1;
return
end
end
if isnumeric(A) && isequal(length(A),1) && ismember(A,[1 4 12 52])
B = 1;
end

54
matlab/utilities/dates/isstringdate.m Normal file
View File

@ -0,0 +1,54 @@
function b = isstringdate(str) % --*-- Unitary tests --*--
% Tests if the input string can be interpreted as a date.
%
% INPUTS
% o str string.
%
% OUTPUTS
% o b integer scalar, equal to 1 if str can be interpreted as a date or 0 otherwise.
% Copyright (C) 2013 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if ischar(str)
b = isquaterly(str) || isyearly(str) || ismonthly(str) || isweekly(str);
else
b = 0;
end
%@test:1
%$
%$ date_1 = '1950M2';
%$ date_2 = '1950m2';
%$ date_3 = '-1950m2';
%$ date_4 = '1950m52';
%$ date_5 = ' 1950';
%$ date_6 = '1950Y';
%$ date_7 = '-1950a';
%$ date_8 = '1950m ';
%$
%$ t(1) = dyn_assert(isstringdate(date_1),1);
%$ t(2) = dyn_assert(isstringdate(date_2),1);
%$ t(3) = dyn_assert(isstringdate(date_3),1);
%$ t(4) = dyn_assert(isstringdate(date_4),0);
%$ t(5) = dyn_assert(isstringdate(date_5),0);
%$ t(6) = dyn_assert(isstringdate(date_6),1);
%$ t(7) = dyn_assert(isstringdate(date_7),1);
%$ t(8) = dyn_assert(isstringdate(date_8),0);
%$ T = all(t);
%@eof:1

24
matlab/utilities/dseries/from.m
View File

@ -1,4 +1,4 @@
function from(varargin)
function from(varargin) % --*-- Unitary tests --*--
% Copyright (C) 2014 Dynare Team
%
@ -219,10 +219,10 @@ for i=1:number_of_variables
end
if i>1
if ismember(var.name,variable_names)
error('dseries::from: All the dseries objects should contain variables with different names!')
else
variable_names(i) = {var.name{1}};
% Locally change variable name.
var = var.rename(var.name{1},get_random_string(20));
end
variable_names(i) = {var.name{1}};
else
variable_names(i) = {var.name{1}};
end
@ -421,4 +421,18 @@ function i = isassignedvariable(var,expr)
return
end
end
i = 0;
i = 0;
%@test:1
%$ try
%$ y = dseries(zeros(400,1),dates('1950Q1')) ;
%$ v = dseries(randn(400,1),dates('1950Q1')) ;
%$ u = dseries(randn(400,1),dates('1950Q1')) ;
%$ from 1950Q2 to 2049Q4 do y(t) = (1+.01*u(t))*y(t-1) + v(t)
%$ t(1) = 1;
%$ catch
%$ t(1) = 0;
%$ end
%$
%$ T = all(t);
%@eof:1

21
matlab/utilities/general/isanynan.m Normal file
View File

@ -0,0 +1,21 @@
function yes = isanynan(array)
% Return one if the array contains at least one NaN, 0 otherwise.
% Copyright (C) 2011-2014 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
yes = any(isnan(array(:)));

94
matlab/var_sample_moments.m
View File

@ -1,29 +1,28 @@
function [YtY,XtY,YtX,XtX,Y,X] = ...
var_sample_moments(FirstObservation,LastObservation,qlag,var_trend_order,datafile,varobs,xls_sheet,xls_range)
function var_sample_moments(nlag, var_trend_order, dataset_)%datafile,varobs,xls_sheet,xls_range)
% Computes the sample moments of a VAR model.
%
% The VAR(p) model is defined by:
%
% y_t = \sum_{k=1}^p y_{t-k} A_k + z_t C + e_t for t = 1,...,T
% y_t = \sum_{k=1}^p y_{t-k} A_k + z_t C + e_t for t = 1,...,T
%
% where y_t is a 1*m vector of observed endogenous variables, p is the
% number of lags, A_k is an m*m real matrix, z_t is a 1*q vector of
% exogenous (deterministic) variables, C is a q*m real matrix and
% e_t is a vector of exogenous stochastic shocks. T is the number
% of observations. The deterministic exogenous variables are assumed to
% be a polynomial trend of order q = "var_trend_order".
% of observations. The deterministic exogenous variables are assumed to
% be a polynomial trend of order q = "var_trend_order".
%
% We define:
% We define:
%
% <> Y = (y_1',y_2',...,y_T')' a T*m matrix,
%
% <> x_t = (y_{t-1},y_{t-2},...,y_{t-p},z_t) a 1*(mp+q) row vector,
% <> x_t = (y_{t-1},y_{t-2},...,y_{t-p},z_t) a 1*(mp+q) row vector,
%
% <> X = (x_1',x_2',...,x_T')' a T*(mp+q) matrix,
% <> X = (x_1',x_2',...,x_T')' a T*(mp+q) matrix,
%
% <> E = (e_1',e_2',...,e_T')' a T*m matrix and
%
% <> A = (A_1',A_2',...,A_p',C')' an (mp+q)*m matrix of coefficients.
% <> A = (A_1',A_2',...,A_p',C')' an (mp+q)*m matrix of coefficients.
%
% So that we can equivalently write the VAR(p) model using the following
% matrix representation:
@ -31,18 +30,17 @@ function [YtY,XtY,YtX,XtX,Y,X] = ...
% Y = X * A +E
%
%
% INPUTS
% o FirstObservation [integer] First observation.
% o LastObservation [integer] Last observation.
% o qlag [integer] Number of lags in the VAR model.
% o var_trend_order [integer] Order of the polynomial exogenous trend:
% INPUTS
% o nlag [integer] Number of lags in the VAR model.
% o var_trend_order [integer] Order of the polynomial exogenous trend:
% = -1 no constant and no linear trend,
% = 0 constant and no linear trend,
% = 1 constant and linear trend.
% o dataset_ [dseries] The sample.
%
% OUTPUTS
% OUTPUTS
% o YtY [double] Y'*Y an m*m matrix.
% o XtY [double] X'*Y an (mp+q)*m matrix.
% o XtY [double] X'*Y an (mp+q)*m matrix.
% o YtX [double] Y'*X an m*(mp+q) matrix.
% o XtX [double] X'*X an (mp+q)*(mp+q) matrix.
% o Y [double] Y a T*m matrix.
@ -50,8 +48,11 @@ function [YtY,XtY,YtX,XtX,Y,X] = ...
%
% SPECIAL REQUIREMENTS
% None.
%
% REMARKS
% Outputs are saved in the base workspace (not returned by the function).
% Copyright (C) 2007-2009 Dynare Team
% Copyright (C) 2007-2014 Dynare Team
%
% This file is part of Dynare.
%
@ -68,42 +69,35 @@ function [YtY,XtY,YtX,XtX,Y,X] = ...
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
X = [];
Y = [];
YtY = [];
YtX = [];
XtY = [];
XtX = [];
LastObservation = dataset_.dates(end);
FirstObservation = dataset_.dates(1)+nlag;
data = read_variables(datafile,varobs,[],xls_sheet,xls_range);
NumberOfObservations = LastObservation-FirstObservation+1;
NumberOfVariables = dataset_.vobs;
if qlag > FirstObservation
error('VarSampleMoments :: not enough data to initialize! Try to increase FirstObservation.')
return
end
NumberOfObservations = LastObservation-FirstObservation+1;% This is T.
NumberOfVariables = size(varobs,1);% This is m.
if var_trend_order == -1% No constant no linear trend case.
X = zeros(NumberOfObservations,NumberOfVariables*qlag);
elseif var_trend_order == 0% Constant and no linear trend case.
X = ones(NumberOfObservations,NumberOfVariables*qlag+1);
indx = NumberOfVariables*qlag+1;
elseif var_trend_order == 1;% Constant and linear trend case.
X = ones(NumberOfObservations,NumberOfVariables*qlag+2);
indx = NumberOfVariables*qlag+1:NumberOfVariables*qlag+2;
if isequal(var_trend_order,-1)
% No constant no linear trend case.
X = zeros(NumberOfObservations,NumberOfVariables*nlag);
elseif isequal(var_trend_order,0)
% Constant and no linear trend case.
X = ones(NumberOfObservations,NumberOfVariables*nlag+1);
indx = NumberOfVariables*nlag+1;
elseif isequal(var_trend_order,1)
% Constant and linear trend case.
X = ones(NumberOfObservations,NumberOfVariables*nlag+2);
indx = NumberOfVariables*nlag+1:NumberOfVariables*nlag+2;
else
error('var_sample_moments :: trend must be equal to -1,0 or 1!')
error('Estimation::var_sample_moments: trend must be equal to -1,0 or 1!')
return
end
% I build matrices Y and X
Y = data(FirstObservation:LastObservation,:);
% I build matrices Y and X
Y = dataset_(FirstObservation:LastObservation).data;
for t=1:NumberOfObservations
line = t + FirstObservation-1;
for lag = 1:qlag
X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = data(line-lag,:);
currentdate = FirstObservation+(t-1);
for lag = 1:nlag
X(t,(lag-1)*NumberOfVariables+1:lag*NumberOfVariables) = dataset_(currentdate-lag).data;
end
end
@ -111,7 +105,9 @@ if (var_trend_order == 1)
X(:,end) = transpose(1:NumberOfObservations)
end
YtY = Y'*Y;
YtX = Y'*X;
XtY = X'*Y;
XtX = X'*X;
assignin('base', 'mYY', Y'*Y);
assignin('base', 'mYX', Y'*X);
assignin('base', 'mXY', X'*Y);
assignin('base', 'mXX', X'*X);
assignin('base', 'Ydata', Y);
assignin('base', 'Xdata', X);

16
preprocessor/ComputingTasks.cc
View File

@ -1876,9 +1876,17 @@ EstimationDataStatement::checkPass(ModFileStructure &mod_file_struct, WarningCon
exit(EXIT_FAILURE);
}
if (options_list.string_options.find("file") == options_list.string_options.end())
if ((options_list.string_options.find("file") == options_list.string_options.end()) &&
(options_list.string_options.find("series") == options_list.string_options.end()))
{
cerr << "ERROR: The file option must be passed to the data statement." << endl;
cerr << "ERROR: The file or series option must be passed to the data statement." << endl;
exit(EXIT_FAILURE);
}
if ((options_list.string_options.find("file") != options_list.string_options.end()) &&
(options_list.string_options.find("series") != options_list.string_options.end()))
{
cerr << "ERROR: The file and series options cannot be used simultaneously in the data statement." << endl;
exit(EXIT_FAILURE);
}
}
@ -1887,8 +1895,8 @@ void
EstimationDataStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output, "options_.dataset");
if (options_list.date_options.find("first_obs") == options_list.date_options.end())
output << "options_.dataset.firstobs = options_.initial_period;" << endl;
//if (options_list.date_options.find("first_obs") == options_list.date_options.end())
// output << "options_.dataset.first_obs = options_.initial_period;" << endl;
}
SubsamplesStatement::SubsamplesStatement(const string &name1_arg,

22
preprocessor/DynareBison.yy
View File

@ -90,7 +90,7 @@ class ParsingDriver;
%token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED
%token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
%token CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED
%token DATAFILE FILE DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME
%token <string_val> FLOAT_NUMBER DATES
@ -833,7 +833,7 @@ restriction_expression : expression {driver.check_restriction_expression_constan
restriction_expression_1 : restriction_elem_expression
| restriction_expression_1 restriction_elem_expression
;
;
restriction_elem_expression : COEFF '(' symbol COMMA INT_NUMBER ')'
{ driver.add_positive_restriction_element($3,$5);}
@ -1366,9 +1366,10 @@ data_options_list : data_options_list COMMA data_options
;
data_options : o_file
| o_new_estimation_data_first_obs
| o_last_obs
| o_new_estimation_data_nobs
| o_series
| o_data_first_obs
| o_data_last_obs
| o_data_nobs
| o_xls_sheet
| o_xls_range
;
@ -2493,6 +2494,7 @@ o_simul_seed : SIMUL_SEED EQUAL INT_NUMBER { driver.error("'simul_seed' option i
o_qz_criterium : QZ_CRITERIUM EQUAL non_negative_number { driver.option_num("qz_criterium", $3); };
o_qz_zero_threshold : QZ_ZERO_THRESHOLD EQUAL non_negative_number { driver.option_num("qz_zero_threshold", $3); };
o_file : FILE EQUAL filename { driver.option_str("file", $3); };
o_series : SERIES EQUAL symbol { driver.option_str("series", $3); };
o_datafile : DATAFILE EQUAL filename { driver.option_str("datafile", $3); };
o_nobs : NOBS EQUAL vec_int
{ driver.option_vec_int("nobs", $3); }
@ -2505,12 +2507,9 @@ o_conditional_variance_decomposition : CONDITIONAL_VARIANCE_DECOMPOSITION EQUAL
{ driver.option_vec_int("conditional_variance_decomposition", $3); }
;
o_first_obs : FIRST_OBS EQUAL INT_NUMBER { driver.option_num("first_obs", $3); };
o_new_estimation_data_first_obs : FIRST_OBS EQUAL date_expr
{ driver.option_date("first_obs", $3); }
;
o_last_obs : LAST_OBS EQUAL date_expr
{ driver.option_date("last_obs", $3); }
;
o_data_first_obs : FIRST_OBS EQUAL date_expr { driver.option_date("firstobs", $3); } ;
o_data_last_obs : LAST_OBS EQUAL date_expr { driver.option_date("lastobs", $3); } ;
o_data_nobs : NOBS EQUAL INT_NUMBER { driver.option_num("nobs", $3); };
o_shift : SHIFT EQUAL signed_number { driver.option_num("shift", $3); };
o_shape : SHAPE EQUAL prior_distribution { driver.prior_shape = $3; };
o_mode : MODE EQUAL signed_number { driver.option_num("mode", $3); };
@ -2523,7 +2522,6 @@ o_bounds : BOUNDS EQUAL vec_value_w_inf { driver.option_num("bounds", $3); };
o_domain : DOMAINN EQUAL vec_value { driver.option_num("domain", $3); };
o_interval : INTERVAL EQUAL vec_value { driver.option_num("interval", $3); };
o_variance : VARIANCE EQUAL expression { driver.set_prior_variance($3); }
o_new_estimation_data_nobs : NOBS EQUAL INT_NUMBER { driver.option_num("nobs", $3); };
o_prefilter : PREFILTER EQUAL INT_NUMBER { driver.option_num("prefilter", $3); };
o_presample : PRESAMPLE EQUAL INT_NUMBER { driver.option_num("presample", $3); };
o_lik_algo : LIK_ALGO EQUAL INT_NUMBER { driver.option_num("lik_algo", $3); };

1
preprocessor/DynareFlex.ll
View File

@ -604,6 +604,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
<DYNARE_STATEMENT>simul_replic {return token::SIMUL_REPLIC;}
<DYNARE_STATEMENT>xls_sheet {return token::XLS_SHEET;}
<DYNARE_STATEMENT>xls_range {return token::XLS_RANGE;}
<DYNARE_STATEMENT>series {return token::SERIES;}
<DYNARE_STATEMENT>mh_recover {return token::MH_RECOVER;}
<DYNARE_STATEMENT>planner_discount {return token::PLANNER_DISCOUNT;}
<DYNARE_STATEMENT>calibration {return token::CALIBRATION;}

15
preprocessor/Statement.cc
View File

@ -135,7 +135,20 @@ OptionsList::writeOutput(ostream &output) const
void
OptionsList::writeOutput(ostream &output, const string &option_group) const
{
output << option_group << " = struct();" << endl;
// Initialize option_group as an empty struct iff the field does not exist!
unsigned idx = option_group.find_last_of(".");
if (idx<UINT_MAX)
{
output << option_group << endl;
output << idx << endl;
output << "if ~isfield(" << option_group.substr(0,idx) << ",'" << option_group.substr(idx+1) << "')" << endl;
output << " " << option_group << " = struct();" << endl;
output << "end" << endl;
}
else
{
output << option_group << " = struct();" << endl;
}
for (num_options_t::const_iterator it = num_options.begin();
it != num_options.end(); it++)

12
preprocessor/SymbolTable.cc
View File

@ -270,14 +270,14 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
if (observedVariablesNbr() > 0)
{
output << "options_.varobs = [];" << endl;
int ic = 1;
output << "options_.varobs = cell(1);" << endl;
for (vector<int>::const_iterator it = varobs.begin();
it != varobs.end(); it++)
if (it == varobs.begin())
output << "options_.varobs = '" << getName(*it) << "';" << endl;
else
output << "options_.varobs = char(options_.varobs, '" << getName(*it) << "');" << endl;
{
output << "options_.varobs(" << ic << ") = {'" << getName(*it) << "'};" << endl;
ic++;
}
output << "options_.varobs_id = [ ";
for (vector<int>::const_iterator it = varobs.begin();
it != varobs.end(); it++)

1
tests/Makefile.am
View File

@ -80,6 +80,7 @@ MODFILES = \
fs2000/fs2000e.mod \
fs2000/fs2000_cmaes.mod \
fs2000/fs2000_calib.mod \
fs2000/fs2000_calib_dseries.mod \
fs2000/fs2000_analytic_derivation.mod \
fs2000/fs2000_missing_data.mod \
fs2000/fs2000_sd.mod \

73
tests/fs2000/fs2000_calib_dseries.mod Normal file
View File

@ -0,0 +1,73 @@
// See fs2000.mod in the examples/ directory for details on the model
set_time(1950Q1);
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;
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;
check;
varobs gp_obs gy_obs;
fsdata = dseries('fsdat_simul.m',1950Q1);
fsdata = fsdata{'g@p,y@_obs'}(1960Q1:1980Q4);
fsdata
data(series=fsdata);//, first_obs=1960Q1, last_obs=1980Q4);
calib_smoother(filtered_vars, filter_step_ahead = [3:4]) m P c e W R k d n l y dA;

79
tests/fs2000/fs2000_data.mod Normal file
View File

@ -0,0 +1,79 @@
// See fs2000.mod in the examples/ directory for details on the model
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;
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;
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_.solve_tolf = 1e-12;
set_time(1970Q3); // Interpreted as the first date available in the sample loaded below.
data(file='fsdat_simul.m',first_obs=1971Q1, nobs=40);
estimation(order=1,nobs=192,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8);

77
tests/fs2000/fs2000_dseries_a.mod Normal file
View File

@ -0,0 +1,77 @@
set_time(1950Q3);
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;
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;
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_.solve_tolf = 1e-12;
data(file=fsdat_simul_dseries,first_obs=1950Q3, nobs=192);
estimation(order=1,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8);

80
tests/fs2000/fs2000_dseries_b.mod Normal file
View File

@ -0,0 +1,80 @@
set_time(1950Q3);
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;
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;
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_.solve_tolf = 1e-12;
fsdataset = dseries('fsdat_simul_dseries.m');
fsdataset = fsdataset(1950Q3:1950Q3+191);
data(series=fsdataset);
estimation(order=1,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8);

22
tests/gsa/ls2003.mod
View File

@ -108,7 +108,7 @@ disp(' ');
disp('MC FILTERING(rmse=1), TO MAP THE FIT FROM PRIORS');
disp('Press ENTER to continue'); pause(5);
dynare_sensitivity(nodisplay, datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible
dynare_sensitivity(nodisplay, datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible
load_stab=1, // load prior sample
istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error)
stab=0, // don't plot again stability analysis results
@ -121,7 +121,7 @@ disp('BY USING THE COMBINED CALL');
disp(' ');
disp('dynare_sensitivity(redform=1,')
disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R),')
disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,')
disp('datafile=data_ca1.m,first_obs=8,nobs=79,prefilter=1,')
disp('istart_rmse=2, rmse=1);')
disp(' ');
disp('Press ENTER to continue'); pause(5);
@ -131,7 +131,7 @@ disp('Press ENTER to continue'); pause(5);
//namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
//namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
//namlagendo=(R), // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
//datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
//istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error)
//rmse=1, // do rmse analysis
//);
@ -144,13 +144,13 @@ disp(' ');
disp('Press ENTER to continue'); pause(5);
// run this to generate posterior mode and Metropolis files if not yet done
estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,
estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2,
prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=4, mh_drop=0.6, nodisplay,
bayesian_irf, filtered_vars, smoother) y_obs R_obs pie_obs dq de;
// run this to produce posterior samples of filtered, smoothed and irf variables, if not yet done
//estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3,
//estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3,
// mh_replic=0, mode_file=ls2003_mode, mode_compute=0, load_mh_file, bayesian_irf,
// filtered_vars, smoother, mh_drop=0.6);
@ -163,7 +163,7 @@ disp('Press ENTER to continue'); pause(5);
dynare_sensitivity(nodisplay, pprior=0,Nsam=2048,neighborhood_width=0.2,
mode_file=ls2003_mode, // specifies the mode file where the mode and Hessian are stored
datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
rmse=1);
disp(' ');
@ -182,7 +182,7 @@ disp('RMSE ANALYSIS FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE');
disp(' ');
disp('Press ENTER to continue'); pause(5);
dynare_sensitivity(nodisplay, mode_file=ls2003_mode,
datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
pprior=0,
stab=0,
rmse=1,
@ -195,12 +195,12 @@ disp('THE LAST TWO CALLS COULD BE DONE TOGETHER');
disp('BY USING THE COMBINED CALL');
disp(' ');
disp('dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,')
disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,')
disp('datafile=data_ca1.m,first_obs=8,nobs=79,prefilter=1,')
disp('rmse=1, alpha2_rmse=0, alpha_rmse=0);')
disp(' ');
disp('Press ENTER to continue'); pause(5);
//dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,
//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
//datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
//rmse=1,
//alpha2_rmse=0, // no correlation analysis
//alpha_rmse=0 // no Smirnov sensitivity analysis
@ -210,10 +210,10 @@ disp(' ');
disp('RMSE ANALYSIS FOR POSTERIOR MCMC sample (ppost=1)');
disp('Needs a call to dynare_estimation to load all MH environment');
disp('Press ENTER to continue'); pause(5);
//estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2, mode_file=ls2003_mode, load_mh_file,
//estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, mode_file=ls2003_mode, load_mh_file,
// prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=0, mh_drop=0.6);
dynare_sensitivity(nodisplay, stab=0, // no need for stability analysis since the posterior sample is surely OK
datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
datafile='data_ca1.m',first_obs=8,nobs=79,prefilter=1,
rmse=1,ppost=1);