Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Loading items
Show changes
Showing
with 300 additions and 140 deletions
......@@ -60,7 +60,18 @@ while weight<1
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
if ~options_.ep.stochastic.order
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
else
switch options_.ep.stochastic.algo
case 0
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
case 1
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order);
end
end
end
info.convergence = ~flag;% Equal to one if the perfect foresight solver converged for the current value of weight.
if verbose
......@@ -133,7 +144,18 @@ if weight<1
flag = 1;
end
if flag
[flag,tmp] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
if ~options_.ep.stochastic.order
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul0,exo_simul0,pfm);
else
switch options_.ep.stochastic.algo
case 0
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
case 1
[flag,tmp] = ...
solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order);
end
end
end
info.convergence = ~flag;
if info.convergence
......
......@@ -42,27 +42,12 @@ else
end
ipred = find(lead_lag_incidence(M.maximum_lag,:))';
order_var = dr.order_var;
LQ = true;
Gy = dr.ghx(nstatic+(1:nspred),:);
Gu = dr.ghu(nstatic+(1:nspred),:);
gy(dr.order_var,:) = dr.ghx;
gu(dr.order_var,:) = dr.ghu;
if options.ramsey_policy && options.order == 1 && ~options.linear
options.order = 2;
options.qz_criterium = 1+1e-6;
[dr,info] = stochastic_solvers(oo.dr,0,M,options,oo);
Gyy = dr.ghxx(nstatic+(1:nspred),:);
Guu = dr.ghuu(nstatic+(1:nspred),:);
Gyu = dr.ghxu(nstatic+(1:nspred),:);
Gss = dr.ghs2(nstatic+(1:nspred),:);
gyy(dr.order_var,:) = dr.ghxx;
guu(dr.order_var,:) = dr.ghuu;
gyu(dr.order_var,:) = dr.ghxu;
gss(dr.order_var,:) = dr.ghs2;
LQ = false;
end
ys = oo.dr.ys;
......@@ -83,24 +68,15 @@ mexErrCheck('A_times_B_kronecker_C', err);
Wbar =U/(1-beta);
Wy = Uy*gy/(eye(nspred)-beta*Gy);
Wu = Uy*gu+beta*Wy*Gu;
if LQ
Wyy = Uyygygy/(eye(nspred*nspred)-beta*kron(Gy,Gy));
else
Wyy = (Uy*gyy+Uyygygy+beta*Wy*Gyy)/(eye(nspred*nspred)-beta*kron(Gy,Gy));
end
[Wyygugu, err] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
mexErrCheck('A_times_B_kronecker_C', err);
[Wyygygu,err] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
mexErrCheck('A_times_B_kronecker_C', err);
if LQ
Wuu = Uyygugu+beta*Wyygugu;
Wyu = Uyygygu+beta*Wyygygu;
Wss = beta*Wuu*M.Sigma_e(:)/(1-beta);
else
Wuu = Uy*guu+Uyygugu+beta*(Wy*Guu+Wyygugu);
Wyu = Uy*gyu+Uyygygu+beta*(Wy*Gyu+Wyygygu);
Wss = (Uy*gss+beta*(Wuu*M.Sigma_e(:)+Wy*Gss))/(1-beta);
end
% initialize yhat1 at the steady state
yhat1 = oo.steady_state;
if options.ramsey_policy
......
......@@ -133,7 +133,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
[chck, r, junk]= bytecode('dynamic','evaluate', z, zx, M.params, ys, 1);
mexErrCheck('bytecode', chck);
elseif options.block
[r, data] = feval([M.fname '_dynamic'], z', zx, M.params, ys, M.maximum_lag+1, data);
[r, oo.dr] = feval([M.fname '_dynamic'], z', zx, M.params, ys, M.maximum_lag+1, oo.dr);
else
iyv = M.lead_lag_incidence';
iyr0 = find(iyv(:));
......
......@@ -59,7 +59,12 @@ function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options)
info(2) = NaN;
end
updated_params_flag = max(abs(params1-params)) > 1e-12 || ~isequal(isnan(params1),isnan(params)); %checks whether numbers or NaN changed
if M.param_nbr > 0
updated_params_flag = max(abs(params1-params)) > 1e-12 ...
|| ~isequal(isnan(params1),isnan(params)); %checks whether numbers or NaN changed
else
updated_params_flag = 0
end
h_set_auxiliary_variables = str2func([M.fname '_set_auxiliary_variables']);
if isnan(updated_params_flag) || (updated_params_flag && any(isnan(params(~isnan(params))-params1(~isnan(params))))) %checks if new NaNs were added
......@@ -99,6 +104,10 @@ function [ys,params,info] = evaluate_steady_state_file(ys_init,exo_ss,M,options)
info(2) = residuals'*residuals;
return
end
if any(isnan(residuals))
info(1) = 22;
return
end
elseif ~isempty(options.steadystate_partial)
ssvar = options.steadystate_partial.ssvar;
nov = length(ssvar);
......
function plan = flip_plan(plan, exogenous, endogenous, expectation_type, date, value)
% Adds to the forecast scenario a conditional forecast shock (the path of an endogenous variable is constrained and the values compatible values of the related exogenous variable will be compued)
%
% INPUTS
% o plan [structure] A structure describing the different shocks and the implied variables, the date of the shocks and the path of the shock (forecast scenario).
% The plan structure is created by the functions init_plan, basic_plan and flip_plan
% o exogenous [string] A string containing the name of the endogenous variable with a constrained path.
% o endogenous [string] A string containing the name of the exogenous. This exogenous variable will be en endogenous variable when the conditional forecast will be perform.
% o expectation_type [string] A string indicating the type of expectation: 'surprise' for an unexpected shock, and 'perfect_foresight' for a perfectly anticpated shock.
% o date [dates] The period of the shock
% o value [array of double] A vector of double containing the constrined path on the endogenous variable
%
%
% OUTPUTS
% plan [structure] Returns a structure containing the updated forecast scenario.
%
%
% Copyright (C) 2013-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 ~ischar(expectation_type) || size(expectation_type,1) ~= 1
error(['in flip_plan the fourth argument should be a string containing the simulation type (''perfect_foresight'' or ''surprise'')']);
end
exogenous = strtrim(exogenous);
ix = find(strcmp(exogenous, plan.endo_names));
if isempty(ix)
error(['in flip_plan the second argument ' exogenous ' is not an endogenous variable']);
end;
endogenous = strtrim(endogenous);
iy = find(strcmp(endogenous, plan.exo_names));
if isempty(iy)
error(['in flip_plan the third argument ' endogenous ' is not an exogenous variable']);
end;
sdate = length(date);
if sdate > 1
if date(1) < plan.date(1) || date(end) > plan.date(end)
error(['in flip_plan the fifth argument (date=' date ') must lay inside the plan.date ' plan.date]);
end
else
if date < plan.date(1) || date > plan.date(end)
error(['in flip_plan the fifth argument (date=' date ') must lay iside the plan.date ' plan.date]);
end
end
if ~isempty(plan.shock_vars_)
common_var = find(iy == plan.shock_vars_);
if ~isempty(common_var)
common_date = intersect(date, plan.shock_date_{common_var});
if ~isempty(common_date)
if common_date.length > 1
the_dates = [cell2mat(strings(common_date(1))) ':' cell2mat(strings(common_date(end)))];
else
the_dates = cell2mat(strings(common_date));
end
error(['Impossible case: ' plan.exo_names{plan.shock_vars_(common_var)} ' is used both as a shock and as an endogenous variable to control the path of ' plan.endo_names{ix} ' at the dates ' the_dates]);
end
end
end
if isempty(plan.constrained_vars_)
plan.constrained_vars_ = ix;
plan.options_cond_fcst_.controlled_varexo = iy;
if strcmp(expectation_type, 'perfect_foresight')
plan.constrained_perfect_foresight_ = 1;
else
plan.constrained_perfect_foresight_ = 0;
end
else
plan.constrained_vars_ = [plan.constrained_vars_ ; ix];
plan.options_cond_fcst_.controlled_varexo = [plan.options_cond_fcst_.controlled_varexo ; iy];
if strcmp(expectation_type, 'perfect_foresight')
plan.constrained_perfect_foresight_ = [plan.constrained_perfect_foresight_ ; 1];
else
plan.constrained_perfect_foresight_ = [plan.constrained_perfect_foresight_ ; 0];
end
end
plan.constrained_date_{length(plan.constrained_date_) + 1} = date;
plan.constrained_paths_{length(plan.constrained_paths_) + 1} = value;
......@@ -33,7 +33,7 @@ info = 1;
MetropolisFolder = CheckPath('metropolis',M_.dname);
ModelName = M_.fname;
BaseName = [MetropolisFolder, ModelName];
BaseName = [MetropolisFolder filesep ModelName];
load_last_mh_history_file(MetropolisFolder, ModelName);
......@@ -44,7 +44,7 @@ predicted_mhname = [ BaseName '_mh' int2str(mh_number) '_blck' int2str(bk_number
all_mh_files = dir([BaseName '_mh*_blck*' ]);
[junk,idx] = sort([all_mh_files.datenum]);
mhnamme = all_mh_files(idx(end)).name;
mhname = all_mh_files(idx(end)).name;
if ~strcmpi(mhname,predicted_mhname)
info = 0;
......
......@@ -54,13 +54,32 @@ options_.order = 1;
[xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
if plt_flag
plot_priors(bayestopt_,M_,options_);
plot_priors(bayestopt_,M_,estim_params_,options_)
end
PriorNames = { 'Beta' , 'Gamma' , 'Gaussian' , 'Inverted Gamma' , 'Uniform' , 'Inverted Gamma -- 2' };
if size(M_.param_names,1)==size(M_.param_names_tex,1)% All the parameters have a TeX name.
fidTeX = fopen('priors_data.tex','w+');
fprintf(fidTeX,'%% TeX-table generated by get_prior_info (Dynare).\n');
fprintf(fidTeX,'%% Prior Information\n');
fprintf(fidTeX,['%% ' datestr(now,0)]);
fprintf(fidTeX,' \n');
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|ccccccc} \n');
fprintf(fidTeX,'\\caption{Prior information (parameters)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Prior}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Lower Bound & Upper Bound & LB Untrunc. 80\\%% HPDI & UB Untrunc. 80\\%% HPDI \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Lower Bound & Upper Bound & LB Untrunc. 80\\%% HPDI & UB Untrunc. 80\\%% HPDI \\\\ \n');
fprintf(fidTeX,'\\hline \\endhead \n');
fprintf(fidTeX,'\\hline \\multicolumn{8}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
% Column 1: a string for the name of the prior distribution.
% Column 2: the prior mean.
% Column 3: the prior standard deviation.
......@@ -81,24 +100,43 @@ if size(M_.param_names,1)==size(M_.param_names_tex,1)% All the parameters have a
case { 1 , 5 }
LowerBound = bayestopt_.p3(i);
UpperBound = bayestopt_.p4(i);
if ~isinf(bayestopt_.lb(i))
LowerBound=max(LowerBound,bayestopt_.lb(i));
end
if ~isinf(bayestopt_.ub(i))
UpperBound=min(UpperBound,bayestopt_.ub(i));
end
case { 2 , 4 , 6 }
LowerBound = bayestopt_.p3(i);
if ~isinf(bayestopt_.lb(i))
LowerBound=max(LowerBound,bayestopt_.lb(i));
end
if ~isinf(bayestopt_.ub(i))
UpperBound=bayestopt_.ub(i);
else
UpperBound = '$\infty$';
end
case 3
if isinf(bayestopt_.p3(i))
if isinf(bayestopt_.p3(i)) && isinf(bayestopt_.lb(i))
LowerBound = '$-\infty$';
else
LowerBound = bayestopt_.p3(i);
if ~isinf(bayestopt_.lb(i))
LowerBound=max(LowerBound,bayestopt_.lb(i));
end
if isinf(bayestopt_.p4(i))
end
if isinf(bayestopt_.p4(i)) && isinf(bayestopt_.ub(i))
UpperBound = '$\infty$';
else
UpperBound = bayestopt_.p4(i);
if ~isinf(bayestopt_.ub(i))
UpperBound=min(UpperBound,bayestopt_.ub(i));
end
end
otherwise
error('get_prior_info:: Dynare bug!')
end
format_string = build_format_string(bayestopt_,i);
format_string = build_format_string(PriorStandardDeviation,LowerBound,UpperBound);
fprintf(fidTeX,format_string, ...
TexName, ...
PriorShape, ...
......@@ -109,6 +147,9 @@ if size(M_.param_names,1)==size(M_.param_names_tex,1)% All the parameters have a
PriorIntervals(i,1), ...
PriorIntervals(i,2) );
end
fprintf(fidTeX,'\\end{longtable}\n ');
fprintf(fidTeX,'\\end{center}\n');
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
......@@ -191,19 +232,19 @@ end
options_.order = order;
function format_string = build_format_string(bayestopt,i)
function format_string = build_format_string(PriorStandardDeviation,LowerBound,UpperBound)
format_string = ['%s & %s & %6.4f &'];
if isinf(bayestopt.p2(i))
if ~isnumeric(PriorStandardDeviation)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if isinf(bayestopt.p3(i))
if ~isnumeric(LowerBound)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
end
if isinf(bayestopt.p4(i))
if ~isnumeric(UpperBound)
format_string = [ format_string , ' %s &'];
else
format_string = [ format_string , ' %6.4f &'];
......
......@@ -48,16 +48,13 @@ n = length(i_var);
[vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
if size(u,2) > 0
i_stat_0 = find(any(abs(A*u) < options_.Schur_vec_tol,2));
i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2));
i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2)); %only set those variances of objective function for which variance is finite
ghx = ghx(i_stat,:);
ghu = ghu(i_stat,:);
else
i_stat_0 = 1:size(ghx,2);
i_stat = (1:n)';
end
vx1 = Inf*ones(n,n);
vx1(i_stat,i_stat) = ghx(:,i_stat_0)*vx(i_stat_0,i_stat_0)*ghx(:,i_stat_0)'+ghu*Sigma_e*ghu';
vx1(i_stat,i_stat) = ghx*vx*ghx'+ghu*Sigma_e*ghu';
......@@ -372,6 +372,7 @@ options_.lik_init = 1;
options_.load_mh_file = 0;
options_.logdata = 0;
options_.loglinear = 0;
options_.logged_steady_state = 0;
options_.mh_conf_sig = 0.90;
options_.prior_interval = 0.90;
options_.mh_drop = 0.5;
......@@ -453,6 +454,7 @@ simplex.tolerance.f = 1e-4;
simplex.maxiter = 5000;
simplex.maxfcallfactor = 500;
simplex.maxfcall = [];
simplex.delta_factor=0.05;
options_.simplex = simplex;
% CMAES optimization routine.
......
......@@ -122,7 +122,7 @@ while j<=MaxNumberOfTuningSimulations
end% ... otherwise I don't move.
prtfrc = j/MaxNumberOfTuningSimulations;
if mod(j, 10)==0
dyn_waitbar(prtfrc,hh,sprintf('Acceptance ratios: %f [%f]',isux/j,jsux/jj));
dyn_waitbar(prtfrc,hh,sprintf('Acceptance ratio [during last 500]: %f [%f]',isux/j,jsux/jj));
end
if j/500 == round(j/500)
test1 = jsux/jj;
......@@ -215,7 +215,7 @@ if strcmpi(info,'LastCall')
end% ... otherwise I don't move.
prtfrc = j/MaxNumberOfTuningSimulations;
if mod(j, 10)==0
dyn_waitbar(prtfrc,hh,sprintf('Acceptance rates: %f [%f]',isux/j,jsux/jj));
dyn_waitbar(prtfrc,hh,sprintf('Acceptance ratio [during last 1000]: %f [%f]',isux/j,jsux/jj));
end
if j/1000 == round(j/1000)
test1 = jsux/jj;
......
......@@ -524,17 +524,17 @@ else
% plot trade-offs
if ~options_.nograph
a00=jet(size(vvarvecm,1));
for ix=1:ceil(length(nsnam)/5),
if options_.opt_gsa.ppost
temp_name='RMSE Posterior Tradeoffs: Log Posterior';
temp_name='RMSE Posterior Tradeoffs:';
else
if options_.opt_gsa.pprior
temp_name='RMSE Prior Tradeoffs: Log Posterior';
temp_name='RMSE Prior Tradeoffs:';
else
temp_name='RMSE MC Tradeoffs: Log Posterior';
temp_name='RMSE MC Tradeoffs;';
end
end
hh = dyn_figure(options_,'name',[temp_name,' ',int2str(ix)]);
for ix=1:ceil(length(nsnam)/5),
hh = dyn_figure(options_,'name',[temp_name,' observed variables ',int2str(ix)]);
for j=1+5*(ix-1):min(size(snam2,1),5*ix),
subplot(2,3,j-5*(ix-1))
%h0=cumplot(x(:,nsnam(j)+nshock));
......@@ -617,7 +617,7 @@ else
fnam = ['rmse_mc_',deblank(vvarvecm(i,:))];
end
end
stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,pvalue,fnam, OutDir);
stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,pvalue,fnam, OutDir,[],[temp_name ' observed variable ' deblank(vvarvecm(i,:))]);
% [pc,latent,explained] = pcacov(c0);
% %figure, bar([explained cumsum(explained)])
......
......@@ -209,7 +209,7 @@ if opt_gsa.morris==1,
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAvdec','vdec','ir_vdec','ic_vdec')
end
hh = dyn_figure(options_);
hh = dyn_figure(options_,'name','Screening identification: variance decomposition');
% boxplot(SAvdec,'whis',10,'symbol','r.')
myboxplot(SAvdec,[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
......@@ -221,7 +221,7 @@ if opt_gsa.morris==1,
text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
xlabel(' ')
title('All variance decomposition')
title('Elementary effects variance decomposition')
dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_vdec'],options_);
else
save([OutputDirectoryName,'/',fname_,'_morris_IDE'],'vdec')
......@@ -314,7 +314,7 @@ if opt_gsa.morris==1,
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'ac','ir_ac','ic_ac')
end
hh=dyn_figure(options_);
hh=dyn_figure(options_,'name','Screening identification: theoretical moments');
% boxplot(SAcc,'whis',10,'symbol','r.')
myboxplot(SAcc,[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
......@@ -326,7 +326,7 @@ if opt_gsa.morris==1,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
xlabel(' ')
title('EET All moments')
title('Elementary effects in the moments')
dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_moments'],options_);
% close(gcf),
......@@ -709,7 +709,7 @@ if opt_gsa.morris==1,
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm','SAmunorm','SAsignorm')
end
hh=dyn_figure(options_); %bar(SAnorm(:,irel))
hh=dyn_figure(options_,'name','Screening identification: model'); %bar(SAnorm(:,irel))
% boxplot(SAnorm','whis',10,'symbol','r.')
myboxplot(SAnorm',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
......@@ -725,35 +725,35 @@ if opt_gsa.morris==1,
title('Elementary effects in the model')
dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morris_par'],options_);
hh=dyn_figure(options_); %bar(SAmunorm(:,irel))
% boxplot(SAmunorm','whis',10,'symbol','r.')
myboxplot(SAmunorm',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'xlim',[0.5 npT+0.5])
set(gca,'ylim',[-1 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
xlabel(' ')
for ip=1:npT,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
xlabel(' ')
title('\mu in the model')
dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_);
hh=dyn_figure(options_); %bar(SAsignorm(:,irel))
% boxplot(SAsignorm','whis',10,'symbol','r.')
myboxplot(SAsignorm',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
set(gca,'xlim',[0.5 npT+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
xlabel(' ')
for ip=1:npT,
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
xlabel(' ')
title('\sigma in the model')
dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_);
% hh=dyn_figure(options_); %bar(SAmunorm(:,irel))
% % boxplot(SAmunorm','whis',10,'symbol','r.')
% myboxplot(SAmunorm',[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
% set(gca,'xlim',[0.5 npT+0.5])
% set(gca,'ylim',[-1 1])
% set(gca,'position',[0.13 0.2 0.775 0.7])
% xlabel(' ')
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title('\mu in the model')
% dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrismu_par'],options_);
%
% hh=dyn_figure(options_); %bar(SAsignorm(:,irel))
% % boxplot(SAsignorm','whis',10,'symbol','r.')
% myboxplot(SAsignorm',[],'.',[],10)
% set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
% set(gca,'xlim',[0.5 npT+0.5])
% set(gca,'ylim',[0 1])
% set(gca,'position',[0.13 0.2 0.775 0.7])
% xlabel(' ')
% for ip=1:npT,
% text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
% end
% xlabel(' ')
% title('\sigma in the model')
% dyn_saveas(hh,[OutputDirectoryName,'/',fname_,'_morrissig_par'],options_);
% figure, bar(SAnorm(:,irel)')
% set(gca,'xtick',[1:j0])
......
......@@ -17,19 +17,23 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
% 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_ estim_params_ oo_
[nr1, nc1, nsam] = size(mm);
nobs=size(options_.varobs,1);
disp('Computing theoretical moments ...')
h = dyn_waitbar(0,'Theoretical moments ...');
vdec = zeros(nobs,M_.exo_nbr,nsam);
cc = zeros(nobs,nobs,nsam);
ac = zeros(nobs,nobs*options_.ar,nsam);
for j=1:nsam,
dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
oo_.dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
if ~isempty(ss),
set_shocks_param(ss(j,:));
end
[vdec(:,:,j), corr, autocorr, z, zz] = th_moments(dr,options_.varobs);
[vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs);
cc(:,:,j)=triu(corr);
dum=[];
for i=1:options_.ar
......
......@@ -132,7 +132,7 @@ for j=1:size(anamendo,1)
iplo=iplo+1;
js=js+1;
xdir0 = [adir,filesep,namendo,'_vs_', namexo];
if ilog==0,
if ilog==0 || ~isempty(threshold),
if isempty(threshold)
if isempty(dir(xdir0))
mkdir(xdir0)
......@@ -233,7 +233,7 @@ for j=1:size(anamendo,1)
iplo=iplo+1;
js=js+1;
xdir0 = [adir,filesep,namendo,'_vs_', namlagendo];
if ilog==0,
if ilog==0 || ~isempty(threshold),
if isempty(threshold)
if isempty(dir(xdir0))
mkdir(xdir0)
......
......@@ -257,6 +257,9 @@ if fload==0,
%try stoch_simul([]);
try
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
if info(1)==0,
info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_);
end
infox(j,1)=info(1);
if infox(j,1)==0 && ~exist('T'),
dr_=oo_.dr;
......@@ -468,23 +471,23 @@ if pprior
aunstablename=[aname, '_unst']; aunstabletitle='Prior StabMap: Parameter driving explosiveness of solution';
awronguniname=[aname, '_wrong']; awrongunititle='Prior StabMap: Parameter driving inability to find solution';
% bivariate
auname='prior_unacceptable'; autitle='Prior Unacceptable';
aunstname='prior_unstable'; aunsttitle='Prior Unstable';
aindname='prior_indeterm'; aindtitle='Prior Indeterminacy';
awrongname='prior_wrong'; awrongtitle='Prior No Solution Found';
asname='prior_stable'; astitle='Prior Stable';
auname='prior_unacceptable'; autitle='Prior StabMap: non-existence of unique stable solution (Unacceptable)';
aunstname='prior_unstable'; aunsttitle='Prior StabMap: explosiveness of solution';
aindname='prior_indeterm'; aindtitle='Prior StabMap: Indeterminacy';
awrongname='prior_wrong'; awrongtitle='Prior StabMap: inability to find solution';
asname='prior_stable'; astitle='Prior StabMap: unique Stable Saddle-Path';
else
% univariate
aname='mc_stab'; atitle='Posterior StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)';
aindetname=[aname, '_indet']; aindettitle='Posterior StabMap: Parameter driving indeterminacy';
aunstablename=[aname, '_unst']; aunstabletitle='Posterior StabMap: Parameter driving explosiveness of solution';
awronguniname=[aname, '_wrong']; awrongunititle='Posterior StabMap: Parameter driving inability to find solution';
aname='mc_stab'; atitle='MC (around posterior mode) StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)';
aindetname=[aname, '_indet']; aindettitle='MC (around posterior mode) StabMap: Parameter driving indeterminacy';
aunstablename=[aname, '_unst']; aunstabletitle='MC (around posterior mode) StabMap: Parameter driving explosiveness of solution';
awronguniname=[aname, '_wrong']; awrongunititle='MC (around posterior mode) StabMap: Parameter driving inability to find solution';
% bivariate
auname='mc_unacceptable'; autitle='Posterior Unacceptable';
aunstname='mc_unstable'; aunsttitle='Posterior Unstable';
aindname='mc_indeterm'; aindtitle='Posterior Indeterminacy';
awrongname='mc_wrong'; awrongtitle='Posterior No Solution Found';
asname='mc_stable'; astitle='Posterior Stable';
auname='mc_unacceptable'; autitle='MC (around posterior mode) StabMap: non-existence of unique stable solution (Unacceptable)';
aunstname='mc_unstable'; aunsttitle='MC (around posterior mode) StabMap: explosiveness of solution';
aindname='mc_indeterm'; aindtitle='MC (around posterior mode) StabMap: Indeterminacy';
awrongname='mc_wrong'; awrongtitle='MC (around posterior mode) StabMap: inability to find solution';
asname='mc_stable'; astitle='MC (around posterior mode) StabMap: Unique Stable Saddle-Path';
end
delete([OutputDirectoryName,filesep,fname_,'_',aname,'_*.*']);
%delete([OutputDirectoryName,filesep,fname_,'_',aname,'_SA_*.*']);
......@@ -533,6 +536,9 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
if any(infox==30),
disp([' For ',num2str(length(find(infox==30))/Nsam*100,'%1.3f'),'\% Ergodic variance can''t be computed.'])
end
if any(infox==49),
disp([' For ',num2str(length(find(infox==49))/Nsam*100,'%1.3f'),'\% The model violates one (many) endogenous prior restriction(s).'])
end
end
skipline()
......@@ -597,7 +603,9 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
c0=corrcoef(lpmat(istable,:));
c00=tril(c0,-1);
if length(istable)>10,
stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
end
if length(iunstable)>10,
stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1,autitle);
end
......
......@@ -86,10 +86,9 @@ for j=1:npar,
fprintf(1,'%20s: corrcoef = %7.3f\n',tmp_name,c0(i2(jx),j));
if ~options_.nograph,
if mod(j2,12)==1,
ifig=ifig+1;
hh=dyn_figure(options_,'name',['Correlations in the ',figtitle,' sample ', num2str(ifig)]);
hh=dyn_figure(options_,'name',[figtitle,' sample bivariate projection ', num2str(ifig)]);
end
subplot(3,4,j2-(ifig-1)*12)
% bar(c0(i2,j)),
......
......@@ -35,11 +35,12 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
end
end
[gamma_y,ivar] = th_autocovariances(dr,ivar,M_, options_);
m = dr.ys(ivar);
[gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
m = dr.ys(ivar(stationary_vars));
i1 = find(abs(diag(gamma_y{1})) > 1e-12);
% i1 = find(abs(diag(gamma_y{1})) > 1e-12);
i1 = [1:length(ivar)];
s2 = diag(gamma_y{1});
sd = sqrt(s2);
......
......@@ -225,7 +225,11 @@ if info(1)==0,
deltaM = deltaM.*abs(params');
deltaM(params==0)=deltaM_prior(params==0);
quant = siJ./repmat(sqrt(diag(cmm)),1,nparam);
if size(quant,1)==1,
siJnorm = abs(quant).*normaliz1;
else
siJnorm = vnorm(quant).*normaliz1;
end
% siJnorm = vnorm(siJ(inok,:)).*normaliz;
quant=[];
% inok = find((abs(TAU)<1.e-8));
......@@ -238,7 +242,11 @@ if info(1)==0,
siH=siH(iy,:);
if ~isempty(iy),
quant = siH./repmat(sqrt(diag_chh(iy)),1,nparam);
if size(quant,1)==1,
siHnorm = abs(quant).*normaliz1;
else
siHnorm = vnorm(quant).*normaliz1;
end
else
siHnorm = [];
end
......@@ -254,7 +262,11 @@ if info(1)==0,
siLRE=siLRE(iy,:);
if ~isempty(iy),
quant = siLRE./repmat(sqrt(diag_clre(iy)),1,np);
if size(quant,1)==1,
siLREnorm = abs(quant).*normaliz1(offset+1:end);
else
siLREnorm = vnorm(quant).*normaliz1(offset+1:end);
end
else
siLREnorm=[];
end
......
......@@ -45,7 +45,11 @@ function [condJ, ind0, indnoJ, ixnoJ, McoJ, PcoJ, jweak, jweak_pair] = identific
npar = size(JJ,2);
indnoJ = zeros(1,npar);
if size(JJ,1)>1,
ind1 = find(vnorm(JJ)>=eps); % take non-zero columns
else
ind1 = find(abs(JJ)>=eps); % take non-zero columns
end
JJ1 = JJ(:,ind1);
[eu,ee2,ee1] = svd( JJ1, 0 );
condJ= cond(JJ1);
......
function imcforecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight)
function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
% Computes conditional forecasts.
%
% INPUTS
......@@ -26,7 +26,7 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst, con
% [1] Results are stored in a structure which is saved in a mat file called conditional_forecasts.mat.
% [2] Use the function plot_icforecast to plot the results.
% Copyright (C) 2006-2013 Dynare Team
% Copyright (C) 2006-2014 Dynare Team
%
% This file is part of Dynare.
%
......@@ -45,15 +45,6 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst, con
global options_ oo_ M_ bayestopt_
if isfield(options_cond_fcst, 'simulation_type')
if strcmp(options_cond_fcst.simulation_type, 'deterministic')
disp('deterministic condtional forecast');
det_cond_forecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight);
return;
end
end
if ~isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set)
options_cond_fcst.parameter_set = 'posterior_mode';
end
......@@ -153,9 +144,14 @@ end
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
[T,R,ys,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
sQ = sqrt(M_.Sigma_e);
if ~isdiagonal(M_.Sigma_e)
warning(sprintf('The innovations are correlated (the covariance matrix has non zero off diagonal elements), the results of the conditional forecasts will\ndepend on the ordering of the innovations (as declared after varexo) because a Cholesky decomposition is used to factorize the covariance matrix.\n\n=> It is preferable to declare the correlations in the model block (explicitly imposing the identification restrictions), unless you are satisfied\nwith the implicit identification restrictions implied by the Cholesky decomposition.'))
end
sQ = chol(M_.Sigma_e,'lower');
NumberOfStates = length(InitState);
FORCS1 = zeros(NumberOfStates,options_cond_fcst.periods+1,options_cond_fcst.replic);
......