Use a structure for the dataset. Bug fixes.

time-shift
Stéphane Adjemian (Scylla) 2011-09-17 15:38:49 +02:00
parent e3bf749826
commit e9764d0538
6 changed files with 185 additions and 189 deletions

View File

@ -1,39 +0,0 @@
function [data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,nvarobs)
% Copyright (C) 2008-2009 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/>.
[variable_index,observation_index] = find(~isnan(data));
data_index = cell(1,gend);
missing_observations_counter = NaN(gend,1);
for obs=1:gend
idx = find(observation_index==obs);
tmp = variable_index(idx);
missing_observations_counter(obs,1) = nvarobs-length(tmp);
data_index(obs) = { tmp(:) };
end
missing_observations_counter = cumsum(missing_observations_counter);
number_of_observations = length(variable_index);
if ~missing_observations_counter
no_more_missing_observations = 0;
else
tmp = find(missing_observations_counter>=(gend*nvarobs-number_of_observations));
no_more_missing_observations = tmp(1);
end

View File

@ -55,6 +55,8 @@ addpath([dynareroot '/particle/'])
addpath([dynareroot '/gsa/'])
addpath([dynareroot '/utilities/doc/'])
addpath([dynareroot '/utilities/tests/'])
addpath([dynareroot '/utilities/dataset/'])
addpath([dynareroot '/utilities/general/'])
% For functions that exist only under some Octave versions
% or some MATLAB versions, and for which we provide some replacement functions

View File

@ -1,10 +1,10 @@
function dynare_estimation(var_list,varargin)
function dynare_estimation(var_list,dname)
% function dynare_estimation(var_list)
% runs the estimation of the model
%
%
% INPUTS
% var_list: selected endogenous variables vector
%
%
% OUTPUTS
% none
%
@ -47,19 +47,20 @@ end
nnobs = length(nobs);
horizon = options_.forecast;
if nargin<2 || ~exist(dname) || isempty(dname)
dname = M_.fname;
end
M_.dname = dname;
if nnobs > 1
if nargin > 1
dname = vargin{1};
else
dname = M_.fname;
end
for i=1:nnobs
options_.nobs = nobs(i);
dynare_estimation_1(var_list,[dname '_' int2str(nobs(i))]);
oo_recursive_{nobs(i)} = oo_;
end
else
dynare_estimation_1(var_list,varargin{:});
dynare_estimation_1(var_list,dname);
end
if nnobs > 1 && horizon > 0
@ -69,12 +70,12 @@ if nnobs > 1 && horizon > 0
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
% Take the log of the variables if needed
if options_.loglinear && ~options_.logdata % and if the data are not in logs, then...
rawdata = log(rawdata);
rawdata = log(rawdata);
end
endo_names = M_.endo_names;
n_varobs = size(options_.varobs,1);
if isempty(var_list)
var_list = endo_names;
nvar = size(endo_names,1);
@ -89,7 +90,7 @@ if nnobs > 1 && horizon > 0
else
error(['Estimation: ' var_list(i,:) ' isn''t an endogenous' ...
'variable'])
end
end
end
end
@ -100,8 +101,8 @@ if nnobs > 1 && horizon > 0
end
if ~isempty(iobs)
IdObs(j,1) = iobs;
end
end
end
end
k = 3+min(nobs(end)-nobs(1)+horizon, ...
size(rawdata,1)-nobs(1));
@ -125,7 +126,7 @@ if nnobs > 1 && horizon > 0
offsetx = 3;
else
offsetx = 0;
end
end
vname = deblank(var_list(i,:));
maxlag = M_.maximum_lag;
for j=1:nnobs

View File

@ -1,11 +1,11 @@
function dynare_estimation_1(var_list_,dname)
% function dynare_estimation_1(var_list_,dname)
% runs the estimation of the model
%
%
% INPUTS
% var_list_: selected endogenous variables vector
% dname: alternative directory name
%
%
% OUTPUTS
% none
%
@ -31,21 +31,20 @@ function dynare_estimation_1(var_list_,dname)
global M_ options_ oo_ estim_params_ bayestopt_
if nargin > 1
[data,rawdata,xparam1] = dynare_estimation_init(var_list_,dname);
else
[data,rawdata,xparam1] = dynare_estimation_init(var_list_);
end
[dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_, fake] = dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
%% Set various options.
options_ = set_default_option(options_,'mh_nblck',2);
data = dataset_.data;
rawdata = dataset_.rawdata;
% Set various options.
options_ = set_default_option(options_,'mh_nblck',2);
options_ = set_default_option(options_,'nodiagnostic',0);
%% Set number of observations
% Set number of observations
gend = options_.nobs;
%% Set the number of observed variables.
% Set the number of observed variables.
n_varobs = size(options_.varobs,1);
%% Get the number of parameters to be estimated.
% 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).
ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters).
@ -65,7 +64,7 @@ if ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
load(options_.mode_file);
end
%% Compute the steady state:
%% Compute the steady state:
if ~isempty(estim_params_)
set_parameters(xparam1);
end
@ -73,7 +72,7 @@ if options_.steadystate_flag% if the *_steadystate.m file is provided.
[ys,tchek] = feval([M_.fname '_steadystate'],...
[zeros(M_.exo_nbr,1);...
oo_.exo_det_steady_state]);
if size(ys,1) < M_.endo_nbr
if size(ys,1) < M_.endo_nbr
if length(M_.aux_vars) > 0
ys = add_auxiliary_variables_to_steadystate(ys,M_.aux_vars,...
M_.fname,...
@ -87,7 +86,7 @@ if options_.steadystate_flag% if the *_steadystate.m file is provided.
end
oo_.steady_state = ys;
else% if the steady state file is not provided.
[dd,info] = resol(oo_.steady_state,0);
[dd,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.steady_state = dd.ys; clear('dd');
end
if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
@ -105,7 +104,7 @@ if options_.dsge_var
if options_.noconstant
evalin('base',...
['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
'var_sample_moments(options_.first_obs,' ...
'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!
@ -116,8 +115,11 @@ if options_.dsge_var
end
end
[data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
missing_value = ~(number_of_observations == gend*n_varobs);
% [data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
missing_value = dataset_.missing.state; %~(number_of_observations == gend*n_varobs);
data_index = [];
number_of_observations = gend*n_varobs;
no_more_missing_observations = [];
initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
@ -254,15 +256,15 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
else% The covariance matrix is initialized with the prior
% covariance (a diagonal matrix) %%Except for infinite variances ;-)
varinit = 'prior';
if strcmpi(varinit,'prior')
if strcmpi(varinit,'prior')
stdev = bayestopt_.p2;
indx = find(isinf(stdev));
stdev(indx) = ones(length(indx),1)*sqrt(10);
vars = stdev.^2;
CovJump = diag(vars);
elseif strcmpi(varinit,'eye')
vars = ones(length(bayestopt_.p2),1)*0.1;
CovJump = diag(vars);
vars = ones(length(bayestopt_.p2),1)*0.1;
CovJump = diag(vars);
else
disp('gmhmaxlik :: Error!')
return
@ -312,7 +314,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
[xparam1,PostVar,Scale,PostMean] = ...
gmhmaxlik('DsgeVarLikelihood',xparam1,[lb ub],...
options_.Opt6Numb,Scale,flag,PostMean,PostVar,gend);
fval = DsgeVarLikelihood(xparam1,gend);
fval = DsgeVarLikelihood(xparam1,gend);
end
options_.mh_jscale = Scale;
mouvement = max(max(abs(PostVar-OldPostVar)));
@ -325,7 +327,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
save([M_.fname '_mode.mat'],'xparam1','hh');
save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
bayestopt_.jscale = ones(length(xparam1),1)*Scale;%??!
end
end
case 7
optim_options = optimset('display','iter','MaxFunEvals',1000000,'MaxIter',6000,'TolFun',1e-8,'TolX',1e-6);
if isfield(options_,'optim_opt')
@ -348,7 +350,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
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,fh,[],myoptions,gend,data);
case 102
%simulating annealing
@ -357,7 +359,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
LB = xparam1 - 1;
UB = xparam1 + 1;
neps=10;
% Set input parameters.
% Set input parameters.
maxy=0;
epsilon=1.0e-9;
rt_=.10;
@ -367,11 +369,11 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
maxevl=100000000;
idisp =1;
npar=length(xparam1);
disp(['size of param',num2str(length(xparam1))])
disp(['size of param',num2str(length(xparam1))])
c=.1*ones(npar,1);
%* Set input values of the input/output parameters.*
vm=1*ones(npar,1);
disp(['number of parameters= ' num2str(npar) 'max= ' num2str(maxy) 't= ' num2str(t)]);
disp(['rt_= ' num2str(rt_) 'eps= ' num2str(epsilon) 'ns= ' num2str(ns)]);
@ -384,8 +386,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
disp(['lower bound(lb)', 'initial conditions', 'upper bound(ub)' ]);
disp([LB xparam1 UB]);
disp(['c vector ' num2str(c')]);
% keyboard
% keyboard
if ~options_.dsge_var
[xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(fh,xparam1,maxy,rt_,epsilon,ns,nt ...
,neps,maxevl,LB,UB,c,idisp ,t,vm,gend,data,data_index,number_of_observations,no_more_missing_observations);
@ -480,9 +482,9 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
for i = 1:nx
tstath(i) = abs(xparam1(i))/stdh(i);
end
header_width = row_header_width(M_,estim_params_,bayestopt_);
tit1 = sprintf('%-*s %7s %8s %7s %6s %4s %6s\n',header_width-2,' ','prior mean', ...
'mode','s.d.','t-stat','prior','pstdev');
if np
@ -497,7 +499,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.p2(ip)));
eval(['oo_.posterior_mode.parameters.' name ' = xparam1(ip);']);
eval(['oo_.posterior_std.parameters.' name ' = stdh(ip);']);
eval(['oo_.posterior_std.parameters.' name ' = stdh(ip);']);
ip = ip+1;
end
end
@ -511,10 +513,10 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
header_width,name,bayestopt_.p1(ip),xparam1(ip), ...
stdh(ip),tstath(ip),pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.p2(ip)));
bayestopt_.p2(ip)));
M_.Sigma_e(k,k) = xparam1(ip)*xparam1(ip);
eval(['oo_.posterior_mode.shocks_std.' name ' = xparam1(ip);']);
eval(['oo_.posterior_std.shocks_std.' name ' = stdh(ip);']);
eval(['oo_.posterior_std.shocks_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
@ -530,7 +532,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
pnames(bayestopt_.pshape(ip)+1,:), ...
bayestopt_.p2(ip)));
eval(['oo_.posterior_mode.measurement_errors_std.' name ' = xparam1(ip);']);
eval(['oo_.posterior_std.measurement_errors_std.' name ' = stdh(ip);']);
eval(['oo_.posterior_std.measurement_errors_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
@ -549,7 +551,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
M_.Sigma_e(k1,k2) = xparam1(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
eval(['oo_.posterior_mode.shocks_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.posterior_std.shocks_corr.' NAME ' = stdh(ip);']);
eval(['oo_.posterior_std.shocks_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
@ -566,7 +568,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
header_width,bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip), ...
pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.p2(ip)));
eval(['oo_.posterior_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.posterior_std.measurement_errors_corr.' NAME ' = stdh(ip);']);
eval(['oo_.posterior_std.measurement_errors_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
@ -582,7 +584,7 @@ if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
md_Laplace = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess ...
- DsgeVarLikelihood(xparam1,gend);
end
oo_.MarginalDensity.LaplaceApproximation = md_Laplace;
oo_.MarginalDensity.LaplaceApproximation = md_Laplace;
disp(' ')
disp(sprintf('Log data density [Laplace approximation] is %f.',md_Laplace))
disp(' ')
@ -605,10 +607,10 @@ elseif ~any(bayestopt_.pshape > 0) && options_.mh_posterior_mode_estimation
disp(sprintf('%-*s %8.4f %7.4f %7.4f', ...
header_width,name,xparam1(ip),stdh(ip),tstath(ip)));
eval(['oo_.mle_mode.parameters.' name ' = xparam1(ip);']);
eval(['oo_.mle_std.parameters.' name ' = stdh(ip);']);
eval(['oo_.mle_std.parameters.' name ' = stdh(ip);']);
ip = ip+1;
end
end
end
if nvx
ip = 1;
disp('standard deviation of shocks')
@ -619,7 +621,7 @@ elseif ~any(bayestopt_.pshape > 0) && options_.mh_posterior_mode_estimation
disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)));
M_.Sigma_e(k,k) = xparam1(ip)*xparam1(ip);
eval(['oo_.mle_mode.shocks_std.' name ' = xparam1(ip);']);
eval(['oo_.mle_std.shocks_std.' name ' = stdh(ip);']);
eval(['oo_.mle_std.shocks_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
@ -631,7 +633,7 @@ elseif ~any(bayestopt_.pshape > 0) && options_.mh_posterior_mode_estimation
name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)))
eval(['oo_.mle_mode.measurement_errors_std.' name ' = xparam1(ip);']);
eval(['oo_.mle_std.measurement_errors_std.' name ' = stdh(ip);']);
eval(['oo_.mle_std.measurement_errors_std.' name ' = stdh(ip);']);
ip = ip+1;
end
end
@ -648,7 +650,7 @@ elseif ~any(bayestopt_.pshape > 0) && options_.mh_posterior_mode_estimation
M_.Sigma_e(k1,k2) = xparam1(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
eval(['oo_.mle_mode.shocks_corr.' NAME ' = xparam1(ip);']);
eval(['oo_.mle_std.shocks_corr.' NAME ' = stdh(ip);']);
eval(['oo_.mle_std.shocks_corr.' NAME ' = stdh(ip);']);
ip = ip+1;
end
end
@ -697,10 +699,10 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
bayestopt_.p2(ip),...
xparam1(ip),...
stdh(ip));
ip = ip + 1;
ip = ip + 1;
end
fprintf(fidTeX,'\\hline\\hline \n');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (parameters)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:1}\n');
fprintf(fidTeX,'\\end{table}\n');
@ -732,11 +734,11 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
bayestopt_.p1(ip),...
bayestopt_.p2(ip),...
xparam1(ip), ...
stdh(ip));
stdh(ip));
ip = ip+1;
end
fprintf(fidTeX,'\\hline\\hline \n');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (standard deviation of structural shocks)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:2}\n');
fprintf(fidTeX,'\\end{table}\n');
@ -763,15 +765,15 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
idx = strmatch(options_.varobs(estim_params_.var_endo(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,:)), ...
deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
bayestopt_.p1(ip), ...
bayestopt_.p2(ip),...
bayestopt_.p2(ip),...
xparam1(ip),...
stdh(ip));
stdh(ip));
ip = ip+1;
end
fprintf(fidTeX,'\\hline\\hline \n');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (standard deviation of measurement errors)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:3}\n');
fprintf(fidTeX,'\\end{table}\n');
@ -806,7 +808,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
ip = ip+1;
end
fprintf(fidTeX,'\\hline\\hline \n');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of structural shocks)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:4}\n');
fprintf(fidTeX,'\\end{table}\n');
@ -841,7 +843,7 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
ip = ip+1;
end
fprintf(fidTeX,'\\hline\\hline \n');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\end{tabular}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of measurement errors)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:5}\n');
fprintf(fidTeX,'\\end{table}\n');
@ -917,7 +919,7 @@ end
if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape ...
> 0) && options_.load_mh_file)) ...
|| ~options_.smoother ) && M_.endo_nbr^2*gend < 1e7 && options_.partial_information == 0 % to be fixed
|| ~options_.smoother ) && M_.endo_nbr^2*gend < 1e7 && 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,gend,data,data_index,missing_value);
oo_.Smoother.SteadyState = ys;
@ -951,7 +953,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
end
if nbplt == 1
hh = figure('Name','Smoothed shocks');
NAMES = [];
@ -1035,7 +1037,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
else
TeXNAMES = char(TeXNAMES,['$ ' deblank(texname) ' $']);
end
end
end
title(name,'Interpreter','none')
eval(['oo_.SmoothedShocks.' deblank(name) ' = innov(k,:)'';']);
end
@ -1049,14 +1051,14 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:nstar
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_SmoothedShocks%s}\n',M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Smoothed shocks.}');
fprintf(fidTeX,'\\label{Fig:SmoothedShocks:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,'\n');
end
end
end
hh = figure('Name','Smoothed shocks');
set(0,'CurrentFigure',hh)
@ -1068,7 +1070,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
subplot(lr,lc,i);
else
subplot(nr,nc,i);
end
end
plot([1 gend],[0 0],'-r','linewidth',0.5)
hold on
plot(1:gend,innov(k,:),'-k','linewidth',1)
@ -1105,7 +1107,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(NAMES,1);
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_SmoothedShocks%s}\n',M_.fname,int2str(nbplt));
fprintf(fidTeX,'\\caption{Smoothed shocks.}');
@ -1114,7 +1116,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
end
end
%%
%% Smooth observational errors...
@ -1192,7 +1194,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:number_of_plots_to_draw
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_SmoothedObservationErrors%s}\n',M_.fname,int2str(1));
fprintf(fidTeX,'\\caption{Smoothed observation errors.}');
@ -1246,14 +1248,14 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:nstar
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_SmoothedObservationErrors%s}\n',M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Smoothed observation errors.}');
fprintf(fidTeX,'\\label{Fig:SmoothedObservationErrors:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,'\n');
end
end
end
hh = figure('Name','Smoothed observation errors');
set(0,'CurrentFigure',hh)
@ -1265,7 +1267,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
subplot(lr,lc,i);
else
subplot(nr,nc,i);
end
end
plot([1 gend],[0 0],'-r','linewidth',0.5)
hold on
plot(1:gend,measurement_error(index(k),:),'-k','linewidth',1)
@ -1301,7 +1303,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(NAMES,1);
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_SmoothedObservedErrors%s}\n',M_.fname,int2str(nbplt));
fprintf(fidTeX,'\\caption{Smoothed observed errors.}');
@ -1310,9 +1312,9 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
end
end
end
end
%%
%% Historical and smoothed variabes
%%
@ -1322,7 +1324,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
end
if nbplt == 1
hh = figure('Name','Historical and smoothed variables');
NAMES = [];
@ -1365,7 +1367,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:n_varobs
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_HistoricalAndSmoothedVariables%s}\n',M_.fname,int2str(1));
fprintf(fidTeX,'\\caption{Historical and smoothed variables.}');
@ -1374,7 +1376,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
end
else
for plt = 1:nbplt-1
hh = figure('Name','Historical and smoothed variables');
@ -1407,7 +1409,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
else
TeXNAMES = char(TeXNAMES,['$ ' deblank(texname) ' $']);
end
end
end
title(name,'Interpreter','none')
end
eval(['print -depsc2 ' M_.fname '_HistoricalAndSmoothedVariables' int2str(plt) '.eps']);
@ -1420,14 +1422,14 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:nstar
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_HistoricalAndSmoothedVariables%s}\n',M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Historical and smoothed variables.}');
fprintf(fidTeX,'\\label{Fig:HistoricalAndSmoothedVariables:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,'\n');
end
end
end
hh = figure('Name','Historical and smoothed variables');
set(0,'CurrentFigure',hh)
@ -1439,7 +1441,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
subplot(lr,lc,i);
else
subplot(nr,nc,i);
end
end
plot(1:gend,yf(k,:),'-r','linewidth',1)
hold on
plot(1:gend,rawdata(:,k),'-k','linewidth',1)
@ -1476,7 +1478,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(NAMES,1);
fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
end
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_HistoricalAndSmoothedVariables%s}\n',M_.fname,int2str(nbplt));
fprintf(fidTeX,'\\caption{Historical and smoothed variables.}');
@ -1489,7 +1491,7 @@ if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.psha
end
end
if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file
if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file
forecast(var_list_,'smoother');
end

View File

@ -164,7 +164,7 @@ else% If estim_params_ is empty...
estim_params_.np = 0;
end
%% Is there a linear trend in the measurement equation?
% Is there a linear trend in the measurement equation?
if ~isfield(options_,'trend_coeffs') % No!
bayestopt_.with_trend = 0;
else% Yes!
@ -181,22 +181,22 @@ else% Yes!
end
end
%% Set the "size" of penalty.
% Set the "size" of penalty.
bayestopt_.penalty = 1e8;
%% Get informations about the variables of the model.
% Get informations about the variables of the model.
dr = set_state_space(oo_.dr,M_);
oo_.dr = dr;
nstatic = dr.nstatic; % Number of static variables.
npred = dr.npred; % Number of predetermined variables.
nspred = dr.nspred; % Number of predetermined variables in the state equation.
%% Test if observed variables are declared.
% Test if observed variables are declared.
if isempty(options_.varobs)
error('VAROBS is missing')
end
%% Setting resticted state space (observed + predetermined variables)
% Setting resticted state space (observed + predetermined variables)
var_obs_index = [];
k1 = [];
for i=1:n_varobs
@ -204,6 +204,21 @@ for i=1:n_varobs
k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')];
end
% Define union of observed and state variables
k2 = union(var_obs_index',[dr.nstatic+1:dr.nstatic+dr.npred]', 'rows');
% Set restrict_state to postion of observed + state variables in expanded state vector.
oo_.dr.restrict_var_list = k2;
% set mf0 to positions of state variables in restricted state vector for likelihood computation.
[junk,bayestopt_.mf0] = ismember([dr.nstatic+1:dr.nstatic+dr.npred]',k2);
% Set mf1 to positions of observed variables in restricted state vector for likelihood computation.
[junk,bayestopt_.mf1] = ismember(var_obs_index,k2);
% Set mf2 to positions of observed variables in expanded state vector for filtering and smoothing.
bayestopt_.mf2 = var_obs_index;
bayestopt_.mfys = k1;
[junk,ic] = intersect(k2,nstatic+(1:npred)');
oo_.dr.restrict_columns = [ic; length(k2)+(1:nspred-npred)'];
k3 = [];
k3p = [];
if options_.selected_variables_only
@ -262,7 +277,7 @@ else
end;
%% Initialization with unit-root variables.
% Initialization with unit-root variables.
if ~isempty(options_.unit_root_vars)
n_ur = size(options_.unit_root_vars,1);
i_ur = zeros(n_ur,1);
@ -292,7 +307,7 @@ if ~isempty(options_.unit_root_vars)
options_.lik_init = 3;
end % if ~isempty(options_.unit_root_vars)
%% Test if the data file is declared.
% Test if the data file is declared.
if isempty(options_.datafile)
if gsa_flag
data = [];
@ -304,48 +319,61 @@ if isempty(options_.datafile)
end
end
%% If jscale isn't specified for an estimated parameter, use global option options_.jscale, set to 0.2, by default.
% 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.
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
% Set the number of observations (nobs) and build a subsample between first_obs and nobs.
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
% Take the log of the variables if needed
if options_.loglinear % If the model is log-linearized...
if ~options_.logdata % and if the data are not in logs, then...
rawdata = log(rawdata);
end
% Load and transform data.
transformation = [];
if options_.loglinear && ~options_.logdata
transformation = @log;
end
% Test if the observations are real numbers.
if ~isreal(rawdata)
error('There are complex values in the data! Probably a wrong transformation')
end
% Test for missing observations.
options_.missing_data = any(any(isnan(rawdata)));
% Prefilter the data if needed.
if options_.prefilter == 1
if options_.missing_data
bayestopt_.mean_varobs = zeros(n_varobs,1);
for variable=1:n_varobs
rdx = find(~isnan(rawdata(:,variable)));
m = mean(rawdata(rdx,variable));
rawdata(rdx,variable) = rawdata(rdx,variable)-m;
bayestopt_.mean_varobs(variable) = m;
end
else
bayestopt_.mean_varobs = mean(rawdata,1)';
rawdata = rawdata-repmat(bayestopt_.mean_varobs',gend,1);
end
end
% Transpose the dataset array.
data = transpose(rawdata);
xls.sheet = options_.xls_sheet;
xls.range = options_.xls_range;
if nargout>3
%% Compute the steady state:
if ~isfield(options_,'nobs')
options_.nobs = [];
end
dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_obs,options_.nobs,transformation,options_.prefilter,xls);
% $$$ rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
% $$$ % Set the number of observations (nobs) and build a subsample between first_obs and nobs.
% $$$ options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
% $$$ gend = options_.nobs;
% $$$ rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
% $$$ % Take the log of the variables if needed
% $$$ if options_.loglinear % If the model is log-linearized...
% $$$ if ~options_.logdata % and if the data are not in logs, then...
% $$$ rawdata = log(rawdata);
% $$$ end
% $$$ end
% $$$ % Test if the observations are real numbers.
% $$$ if ~isreal(rawdata)
% $$$ error('There are complex values in the data! Probably a wrong transformation')
% $$$ end
% $$$ % Test for missing observations.
% $$$ options_.missing_data = any(any(isnan(rawdata)));
% $$$ % Prefilter the data if needed.
% $$$ if options_.prefilter == 1
% $$$ if options_.missing_data
% $$$ bayestopt_.mean_varobs = zeros(n_varobs,1);
% $$$ for variable=1:n_varobs
% $$$ rdx = find(~isnan(rawdata(:,variable)));
% $$$ m = mean(rawdata(rdx,variable));
% $$$ rawdata(rdx,variable) = rawdata(rdx,variable)-m;
% $$$ bayestopt_.mean_varobs(variable) = m;
% $$$ end
% $$$ else
% $$$ bayestopt_.mean_varobs = mean(rawdata,1)';
% $$$ rawdata = rawdata-repmat(bayestopt_.mean_varobs',gend,1);
% $$$ end
% $$$ end
% $$$ % Transpose the dataset array.
% $$$ data = transpose(rawdata);
if nargout>7
% Compute the steady state:
if options_.steadystate_flag% if the *_steadystate.m file is provided.
[ys,tchek] = feval([M_.fname '_steadystate'],...
[zeros(M_.exo_nbr,1);...
@ -364,7 +392,7 @@ if nargout>3
end
oo_.steady_state = ys;
else% if the steady state file is not provided.
[dd,info] = resol(oo_.steady_state,0);
[dd,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.steady_state = dd.ys; clear('dd');
end
if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
@ -373,15 +401,16 @@ if nargout>3
options_.noconstant = 0;
end
[data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
missing_value = ~(number_of_observations == gend*n_varobs);
% initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
data_info.gend = gend;
data_info.data = data;
data_info.data_index = data_index;
data_info.number_of_observations = number_of_observations;
data_info.no_more_missing_observations = no_more_missing_observations;
data_info.missing_value = missing_value;
fake = [];
% $$$ [data_index,number_of_observations,no_more_missing_observations] = describe_missing_data(data,gend,n_varobs);
% $$$ missing_value = ~(number_of_observations == gend*n_varobs);
% $$$
% $$$ % initial_estimation_checks(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations);
% $$$
% $$$ data_info.gend = gend;
% $$$ data_info.data = data;
% $$$ data_info.data_index = data_index;
% $$$ data_info.number_of_observations = number_of_observations;
% $$$ data_info.no_more_missing_observations = no_more_missing_observations;
% $$$ data_info.missing_value = missing_value;
end

View File

@ -40,7 +40,8 @@ function [A,B,ys,info] = dynare_resolve(mode)
global oo_ M_ options_
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.dr = dr;
if info(1) > 0
A = [];