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 1211 additions and 746 deletions
function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
%function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
function [ny, nx, posterior, prior, forecast_data] = toolbox(nlags)
%function [ny, nx, posterior, prior, forecast_data] = toolbox(nlags)
% bvar_toolbox Routines shared between BVAR methods
% Computes several things for the estimations of a BVAR(nlags)
%
......@@ -42,7 +42,7 @@ function [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags)
% - bvar_prior_{tau,decay,lambda,mu,omega,flat,train}
% Copyright © 2003-2007 Christopher Sims
% Copyright © 2007-2017 Dynare Team
% Copyright © 2007-2023 Dynare Team
%
% This file is part of Dynare.
%
......
function [forcs, e] = mcforecast3(cL, H, mcValue, shocks, forcs, T, R, mv, mu)
function [forcs, e] = get_shock_paths(cL, H, mcValue, shocks, forcs, T, R, mv, mu)
% [forcs, e] = get_shock_paths(cL, H, mcValue, shocks, forcs, T, R, mv, mu)
% Computes the shock values for constrained forecasts necessary to keep
% endogenous variables at their constrained paths
%
......
function plot_icforecast(Variables,periods,options_,oo_)
function plot(Variables,periods,options_,oo_)
% plot(Variables,periods,options_,oo_)
% Build plots for the conditional forecasts.
%
% INPUTS
......@@ -11,7 +12,7 @@ function plot_icforecast(Variables,periods,options_,oo_)
% None.
%
% SPECIAL REQUIREMENTS
% This routine has to be called after imcforecast.m.
% This routine has to be called after conditional_forecasts.run
% Copyright © 2006-2023 Dynare Team
%
......@@ -44,7 +45,7 @@ else
end
if periods>forecast_periods
fprintf('\nplot_icforecast:: Number of periods for plotting exceeds forecast horizon. Setting periods to %d.\n',forecast_periods-1)
fprintf('\nconditional_forecasts.plot:: Number of periods for plotting exceeds forecast horizon. Setting periods to %d.\n',forecast_periods-1)
periods=forecast_periods;
end
......
function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
function forecasts=run(M_,options_,oo_,bayestopt_,estim_params_,constrained_paths, constrained_vars, options_cond_fcst)
% run(constrained_paths, constrained_vars, options_cond_fcst)
% Computes conditional forecasts.
%
% INPUTS
% M_ [structure] describing the model
% options_ [structure] describing the options
% oo_: [structure] storing the results
% bayestopt_ [structure] describing the priors
% estim_params_ [structure] characterizing parameters to be estimated
% - constrained_paths [double] m*p array, where m is the number of constrained endogenous variables and p is the number of constrained periods.
% - constrained_vars [integer] m*1 array, indices in M_.endo_names of the constrained variables.
% - options_cond_fcst [structure] containing the options. The fields are:
......@@ -19,16 +24,16 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
% + conf_sig [double] scalar in [0,1], probability mass covered by the confidence bands.
%
% OUTPUTS
% None.
% - forecasts [structure] results of the conditional forecast exercise
%
% SPECIAL REQUIREMENTS
% This routine has to be called after an estimation statement or an estimated_params block.
%
% REMARKS
% [1] Results are stored in oo_.conditional_forecast.
% [2] Use the function plot_icforecast to plot the results.
% [2] Use the function conditional_forecasts.plot to plot the results.
% Copyright © 2006-2020 Dynare Team
% Copyright © 2006-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -45,8 +50,6 @@ function imcforecast(constrained_paths, constrained_vars, options_cond_fcst)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global options_ oo_ M_ bayestopt_ estim_params_
if ~isfield(options_cond_fcst,'parameter_set') || isempty(options_cond_fcst.parameter_set)
if isfield(oo_,'posterior_mode')
options_cond_fcst.parameter_set = 'posterior_mode';
......@@ -77,7 +80,7 @@ end
if estimated_model
if options_.prefilter
error('imcforecast:: Conditional forecasting does not support the prefiltering option')
error('conditional_forecast: Conditional forecasting does not support the prefiltering option')
end
if ischar(options_cond_fcst.parameter_set)
switch options_cond_fcst.parameter_set
......@@ -100,19 +103,19 @@ if estimated_model
xparam = bayestopt_.p1;
graph_title='Prior Mean';
otherwise
disp('imcforecast:: If the input argument is a string, then it has to be equal to:')
disp('conditional_forecast: If the input argument is a string, then it has to be equal to:')
disp(' ''calibration'', ')
disp(' ''posterior_mode'', ')
disp(' ''posterior_mean'', ')
disp(' ''posterior_median'', ')
disp(' ''prior_mode'' or')
disp(' ''prior_mean''.')
error('imcforecast:: Wrong argument type!')
error('conditional_forecast: Wrong argument type!')
end
else
xparam = options_cond_fcst.parameter_set;
if length(xparam)~=length(M_.params)
error('imcforecast:: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
error('conditional_forecast: The dimension of the vector of parameters doesn''t match the number of estimated parameters!')
end
end
set_parameters(xparam);
......@@ -123,7 +126,6 @@ if estimated_model
missing_value = dataset_info.missing.state;
%store qz_criterium
qz_criterium_old=options_.qz_criterium;
options_=select_qz_criterium_value(options_);
options_smoothed_state_uncertainty_old = options_.smoothed_state_uncertainty;
[atT, ~, ~, ~,ys, ~, ~, ~, ~, ~, ~, ~, ~, ~,oo_,bayestopt_] = ...
......@@ -155,12 +157,11 @@ if estimated_model
trend = constant(oo_.dr.order_var,:);
InitState(:,1) = atT(:,end);
else
qz_criterium_old=options_.qz_criterium;
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
graph_title='Calibration';
if ~isfield(oo_.dr,'kstate')
if ~isfield(oo_.dr,'ghx')
error('You need to call stoch_simul before conditional_forecast')
end
end
......@@ -218,7 +219,7 @@ n2 = size(options_cond_fcst.controlled_varexo,1);
constrained_vars = oo_.dr.inv_order_var(constrained_vars); % must be in decision rule order
if n1 ~= n2
error('imcforecast:: The number of constrained variables doesn''t match the number of controlled shocks')
error('conditional_forecast: The number of constrained variables doesn''t match the number of controlled shocks')
end
% Get indices of controlled varexo.
......@@ -246,7 +247,7 @@ FORCS1_shocks = zeros(n1,cL,options_cond_fcst.replic);
for b=1:options_cond_fcst.replic %conditional forecast using cL set to constrained values
shocks = sQ*randn(ExoSize,options_cond_fcst.periods);
shocks(controlled_varexo,:) = zeros(n1, options_cond_fcst.periods);
[FORCS1(:,:,b), FORCS1_shocks(:,:,b)] = mcforecast3(cL,options_cond_fcst.periods,constrained_paths,shocks,FORCS1(:,:,b),T,R,mv, mu);
[FORCS1(:,:,b), FORCS1_shocks(:,:,b)] = conditional_forecasts.get_shock_paths(cL,options_cond_fcst.periods,constrained_paths,shocks,FORCS1(:,:,b),T,R,mv, mu);
FORCS1(:,:,b)=FORCS1(:,:,b)+trend; %add trend
end
if max(max(max(abs(bsxfun(@minus,FORCS1(constrained_vars,1+1:1+cL,:),trend(constrained_vars,1:cL)+constrained_paths)))))>1e-4
......@@ -292,7 +293,7 @@ FORCS2(:,1,:) = repmat(InitState,1,options_cond_fcst.replic); %set initial stead
for b=1:options_cond_fcst.replic %conditional forecast using cL set to 0
shocks = sQ*randn(ExoSize,options_cond_fcst.periods);
shocks(controlled_varexo,:) = zeros(n1, options_cond_fcst.periods);
FORCS2(:,:,b) = mcforecast3(0,options_cond_fcst.periods,constrained_paths,shocks,FORCS2(:,:,b),T,R,mv, mu)+trend;
FORCS2(:,:,b) = conditional_forecasts.get_shock_paths(0,options_cond_fcst.periods,constrained_paths,shocks,FORCS2(:,:,b),T,R,mv, mu)+trend;
end
mFORCS2 = mean(FORCS2,3);
......@@ -308,7 +309,3 @@ for i = 1:EndoSize
end
forecasts.graph.title = graph_title;
forecasts.graph.fname = M_.fname;
%reset qz_criterium
options_.qz_criterium=qz_criterium_old;
oo_.conditional_forecast = forecasts;
......@@ -28,7 +28,7 @@ function nls(eqname, params, data, range, optimizer, varargin)
% equation must have NaN values in the object.
% [4] It is assumed that the residual is additive.
% Copyright © 2021-2022 Dynare Team
% Copyright © 2021-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -179,12 +179,6 @@ write_residuals_routine(lhs, rhs, eqname, ipnames_, M_);
% Create a routine for evaluating the sum of squared residuals of the nonlinear model
write_ssr_routine(lhs, rhs, eqname, ipnames_, M_);
% Workaround for Octave bug https://savannah.gnu.org/bugs/?46282
% Octave will randomly fail to read the ssr_* file generated in the +folder
if isoctave && octave_ver_less_than('7')
path(path)
end
% Create a function handle returning the sum of square residuals for a given vector of parameters.
ssrfun = @(p) feval([M_.fname '.ssr_' eqname], p, DATA, M_, oo_);
......@@ -298,17 +292,17 @@ end
%
if is_gauss_newton
[params1, SSR, exitflag] = gauss_newton(resfun, params0);
[params1, SSR] = gauss_newton(resfun, params0);
elseif is_lsqnonlin
if ismember('levenberg-marquardt', varargin)
% Levenberg Marquardt does not handle boundary constraints.
[params1, SSR, ~, exitflag] = lsqnonlin(resfun, params0, [], [], optimset(varargin{:}));
[params1, SSR] = lsqnonlin(resfun, params0, [], [], optimset(varargin{:}));
else
[params1, SSR, ~, exitflag] = lsqnonlin(resfun, params0, bounds(:,1), bounds(:,2), optimset(varargin{:}));
[params1, SSR] = lsqnonlin(resfun, params0, bounds(:,1), bounds(:,2), optimset(varargin{:}));
end
else
% Estimate the parameters by minimizing the sum of squared residuals.
[params1, SSR, exitflag] = dynare_minimize_objective(ssrfun, params0, ...
[params1, SSR] = dynare_minimize_objective(ssrfun, params0, ...
minalgo, ...
options_, ...
bounds, ...
......
......@@ -27,7 +27,7 @@ function [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p,
%
% Copyright © 2005 European Commission
% Copyright © 2012-2017 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -50,32 +50,34 @@ if nargin==0
return
end
OutMatrix=[];
if nargin < 5, Group=[]; end
NumGroups = size(Group,2);
if nargin < 4 | isempty(p)
if nargin < 4 || isempty(p)
p = 4;
end
Delt = p/(2*p-2);
if NumGroups ~ 0
if NumGroups ~= 0
sizea = NumGroups; % Number of groups
GroupMat=Group;
GroupMat = GroupMat';
else
sizea = NumFact;
end
r=size(Sample,1)/(sizea+1); % Number of trajectories
if NumGroups > 0
OutMatrix = NaN(sizea*size(Output,2),1);
else
OutMatrix = NaN(sizea*size(Output,2),3);
end
% For Each Output
for k=1:size(Output,2)
OutValues=Output(:,k);
% For each r trajectory
for i=1:r
% For each step j in the trajectory
% Read the orientation matrix fact for the r-th sampling
% Read the corresponding output values
......@@ -86,13 +88,11 @@ for k=1:size(Output,2)
% For each point of the fixed trajectory compute the values of the Morris function. The function
% is partitioned in four parts, from order zero to order 4th.
% SAmeas=NaN();
for j=1:sizea % For each point in the trajectory i.e for each factor
% matrix of factor which changes
if NumGroups ~ 0
if NumGroups ~= 0
AuxFind (:,1) = A(:,j);
% AuxFind(find(A(:,j)),1)=1;
% Pippo = sum((Group - repmat(AuxFind,1,NumGroups)),1);
% Change_factor(j,i) = find(Pippo==0);
Change_factor = find(abs(AuxFind)>1e-010);
% If we deal with groups we can only estimate the new mu*
% measure since factors in the same groups can move in
......@@ -100,7 +100,6 @@ for k=1:size(Output,2)
% Morris mu cannopt be applied.
% In the new version the elementary effect is defined with
% the absolute value.
%SAmeas(find(GroupMat(Change_factor(j,i),:)),i) = abs((Single_OutValues(j) - Single_OutValues(j+1) )/Delt); %(2/3));
SAmeas(i,Change_factor') = abs((Single_OutValues(j) - Single_OutValues(j+1) )/Delt);
else
Change_factor(j,i) = find(Single_Sample(j+1,:)-Single_Sample(j,:));
......@@ -116,12 +115,17 @@ for k=1:size(Output,2)
end %for i=1:r
if NumGroups ~ 0
if NumGroups ~= 0
SAmeas = SAmeas';
end
% Compute Mu AbsMu and StDev
if any(any(isnan(SAmeas)))
AbsMu=NaN(NumFact,1);
if NumGroups == 0
Mu=NaN(NumFact,1);
StDev=NaN(NumFact,1);
end
for j=1:NumFact
SAm = SAmeas(j,:);
SAm = SAm(find(~isnan(SAm)));
......@@ -143,8 +147,8 @@ for k=1:size(Output,2)
% Define the output Matrix - if we have groups we cannot define the old
% measure mu, only mu* makes sense
if NumGroups > 0
OutMatrix = [OutMatrix; AbsMu];
OutMatrix((k-1)*sizea+1:k*sizea,:) = AbsMu;
else
OutMatrix = [OutMatrix; AbsMu, Mu, StDev];
OutMatrix((k-1)*sizea+1:k*sizea,:) = [AbsMu, Mu, StDev];
end
end % For Each Output
function [Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
%[Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
% Inputs: k (1,1) := number of factors examined or number of groups examined.
% Inputs:
% k (1,1) := number of factors examined or number of groups examined.
% In case the groups are chosen the number of factors is stores in NumFact and
% sizea becomes the number of created groups.
% NumFact (1,1) := number of factors examined in the case when groups are chosen
......@@ -34,7 +35,8 @@ function [Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
% AuxMat(sizeb,sizea) := Delta*0.5*((2*B - A) * DD0 + A) in Morris, 1991. The AuxMat is used as in Morris design
% for single factor analysis, while it constitutes an intermediate step for the group analysis.
%
% Output: Outmatrix(sizeb*r, sizea) := for the entire sample size computed In(i,j) matrices
% Outputs:
% Outmatrix(sizeb*r, sizea) := for the entire sample size computed In(i,j) matrices
% OutFact(sizea*r,1) := for the entire sample size computed Fact(i,1) vectors
%
% Note: B0 is constructed as in Morris design when groups are not considered. When groups are considered the routine
......@@ -56,7 +58,7 @@ function [Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
%
% Copyright © 2005 European Commission
% Copyright © 2012-2017 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -80,25 +82,23 @@ Delta = p/(2*p-2);
NumFact = sizea;
GroupNumber = size(GroupMat,2);
if GroupNumber ~ 0;
if GroupNumber ~= 0
sizea = size(GroupMat,2);
end
sizeb = sizea + 1;
sizec = 1;
Outmatrix = [];
OutFact = [];
% For each i generate a trajectory
for i=1:r
% Construct DD0 - OLD VERSION - it does not need communication toolbox
% RAND(N,M) is an NXM matrix with random entries, chosen from a uniform distribution on the interval (0.0,1.0).
% Note that DD0 tells if the factor have to be increased or ddecreased
% by Delta.
randmult = ones(k,1);
v = rand(k,1);
randmult (find(v < 0.5))=-1;
randmult (v < 0.5)=-1;
randmult = repmat(randmult,1,k);
DD0 = randmult .* eye(k);
......@@ -133,7 +133,7 @@ for i=1:r
% When groups are present the random permutation is done only on B. The effect is the same since
% the added part (A0*x0') is completely random.
if GroupNumber ~ 0
if GroupNumber ~= 0
B = B * (GroupMat*P0')';
end
......@@ -147,7 +147,7 @@ for i=1:r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a --> Define the random vector x0 for the factors. Note that x0 takes value in the hypercube
% [0,...,1-Delta]*[0,...,1-Delta]*[0,...,1-Delta]*[0,...,1-Delta]
MyInt = repmat([0:(1/(p-1)):(1-Delta)],NumFact,1); % Construct all possible values of the factors
MyInt = repmat(0:(1/(p-1)):(1-Delta),NumFact,1); % Construct all possible values of the factors
% OLD VERSION - it needs communication toolbox
% w = randint(NumFact,1,[1,size(MyInt,2)]);
......@@ -158,9 +158,9 @@ for i=1:r
% 3) check in which interval the random numbers fall
% 4) generate the corresponding integer
v = repmat(rand(NumFact,1),1,size(MyInt,2)+1); % 1)
IntUsed = repmat([0:1/size(MyInt,2):1],NumFact,1); % 2)
IntUsed = repmat(0:1/size(MyInt,2):1,NumFact,1); % 2)
DiffAuxVec = IntUsed - v; % 3)
w=NaN(1,size(DiffAuxVec,1));
for ii = 1:size(DiffAuxVec,1)
w(1,ii) = max(find(DiffAuxVec(ii,:)<0)); % 4)
end
......@@ -168,7 +168,7 @@ for i=1:r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% b --> Compute the matrix B*, here indicated as B0. Each row in B0 is a
% trajectory for Morris Calculations. The dimension of B0 is (Numfactors+1,Numfactors)
if GroupNumber ~ 0
if GroupNumber ~= 0
B0 = (A0*x0' + AuxMat);
else
B0 = (A0*x0' + AuxMat)*P0;
......@@ -183,6 +183,7 @@ for i=1:r
% Create the Factor vector. Each component of this vector indicate which factor or group of factor
% has been changed in each step of the trajectory.
Fact=NaN(1,sizea);
for j=1:sizea
Fact(1,j) = find(P0(j,:));
end
......@@ -190,5 +191,4 @@ for i=1:r
Outmatrix = [Outmatrix; In];
OutFact = [OutFact; Fact'];
end
\ No newline at end of file
function sout = boxplot (data,notched,symbol,vertical,maxwhisker)
% sout = boxplot (data,notched,symbol,vertical,maxwhisker)
% Creates a box plot
function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
%
% Copyright © 2010-2017 Dynare Team
% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -21,7 +19,6 @@ function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% % % % endif
if nargin < 5 || isempty(maxwhisker), maxwhisker = 1.5; end
if nargin < 4 || isempty(vertical), vertical = 1; end
if nargin < 3 || isempty(symbol), symbol = ['+','o']; end
......@@ -32,7 +29,7 @@ if length(symbol)==1, symbol(2)=symbol(1); end
if notched==1, notched=0.25; end
a=1-notched;
% ## figure out how many data sets we have
% % figure out how many data sets we have
if iscell(data)
nc = length(data);
else
......@@ -40,11 +37,11 @@ else
nc = size(data,2);
end
% ## compute statistics
% ## s will contain
% ## 1,5 min and max
% ## 2,3,4 1st, 2nd and 3rd quartile
% ## 6,7 lower and upper confidence intervals for median
% compute statistics
% s will contain
% 1,5 min and max
% 2,3,4 1st, 2nd and 3rd quartile
% 6,7 lower and upper confidence intervals for median
s = zeros(7,nc);
box = zeros(1,nc);
whisker_x = ones(2,1)*[1:nc,1:nc];
......@@ -55,44 +52,36 @@ outliers2_x = [];
outliers2_y = [];
for i=1:nc
% ## Get the next data set from the array or cell array
% Get the next data set from the array or cell array
if iscell(data)
col = data{i}(:);
else
col = data(:,i);
end
% ## Skip missing data
% Skip missing data
% % % % % % % col(isnan(col) | isna (col)) = [];
col(isnan(col)) = [];
% ## Remember the data length
% Remember the data length
nd = length(col);
box(i) = nd;
if (nd > 1)
% ## min,max and quartiles
% s(1:5,i) = statistics(col)(1:5);
% min,max and quartiles
s(1,i)=min(col);
s(5,i)=max(col);
s(2,i)=myprctilecol(col,25);
s(3,i)=myprctilecol(col,50);
s(4,i)=myprctilecol(col,75);
% ## confidence interval for the median
% confidence interval for the median
est = 1.57*(s(4,i)-s(2,i))/sqrt(nd);
s(6,i) = max([s(3,i)-est, s(2,i)]);
s(7,i) = min([s(3,i)+est, s(4,i)]);
% ## whiskers out to the last point within the desired inter-quartile range
% whiskers out to the last point within the desired inter-quartile range
IQR = maxwhisker*(s(4,i)-s(2,i));
whisker_y(:,i) = [min(col(col >= s(2,i)-IQR)); s(2,i)];
whisker_y(:,nc+i) = [max(col(col <= s(4,i)+IQR)); s(4,i)];
% ## outliers beyond 1 and 2 inter-quartile ranges
% outliers beyond 1 and 2 inter-quartile ranges
outliers = col((col < s(2,i)-IQR & col >= s(2,i)-2*IQR) | (col > s(4,i)+IQR & col <= s(4,i)+2*IQR));
outliers2 = col(col < s(2,i)-2*IQR | col > s(4,i)+2*IQR);
outliers_x = [outliers_x; i*ones(size(outliers))];
......@@ -100,41 +89,37 @@ for i=1:nc
outliers2_x = [outliers2_x; i*ones(size(outliers2))];
outliers2_y = [outliers2_y; outliers2];
elseif (nd == 1)
% ## all statistics collapse to the value of the point
% all statistics collapse to the value of the point
s(:,i) = col;
% ## single point data sets are plotted as outliers.
% single point data sets are plotted as outliers.
outliers_x = [outliers_x; i];
outliers_y = [outliers_y; col];
else
% ## no statistics if no points
% no statistics if no points
s(:,i) = NaN;
end
end
% % % % if isempty(outliers2_y)
% % % % outliers2_y=
% ## Note which boxes don't have enough stats
% Note which boxes don't have enough stats
chop = find(box <= 1);
% ## Draw a box around the quartiles, with width proportional to the number of
% ## items in the box. Draw notches if desired.
% Draw a box around the quartiles, with width proportional to the number of
% items in the box. Draw notches if desired.
box = box*0.23/max(box);
quartile_x = ones(11,1)*[1:nc] + [-a;-1;-1;1;1;a;1;1;-1;-1;-a]*box;
quartile_y = s([3,7,4,4,7,3,6,2,2,6,3],:);
% ## Draw a line through the median
% Draw a line through the median
median_x = ones(2,1)*[1:nc] + [-a;+a]*box;
% median_x=median(col);
median_y = s([3,3],:);
% ## Chop all boxes which don't have enough stats
% Chop all boxes which don't have enough stats
quartile_x(:,chop) = [];
quartile_y(:,chop) = [];
whisker_x(:,[chop,chop+nc]) = [];
whisker_y(:,[chop,chop+nc]) = [];
median_x(:,chop) = [];
median_y(:,chop) = [];
% % % %
% ## Add caps to the remaining whiskers
% Add caps to the remaining whiskers
cap_x = whisker_x;
cap_x(1,:) =cap_x(1,:)- 0.05;
cap_x(2,:) =cap_x(2,:)+ 0.05;
......@@ -144,12 +129,14 @@ cap_y = whisker_y([1,1],:);
% #whisker_x,whisker_y
% #median_x,median_y
% #cap_x,cap_y
%
% ## Do the plot
% Do the plot
mm=min(min(data));
MM=max(max(data));
if isnan(mm), mm=0; MM=0; end
if isnan(mm)
mm=0;
MM=0;
end
if vertical
plot (quartile_x, quartile_y, 'b', ...
......@@ -161,17 +148,30 @@ if vertical
set(gca,'XTick',1:nc);
set(gca, 'XLim', [0.5, nc+0.5]);
set(gca, 'YLim', [mm-(MM-mm)*0.05-eps, MM+(MM-mm)*0.05+eps]);
else
% % % % % plot (quartile_y, quartile_x, "b;;",
% % % % % whisker_y, whisker_x, "b;;",
% % % % % cap_y, cap_x, "b;;",
% % % % % median_y, median_x, "r;;",
% % % % % outliers_y, outliers_x, [symbol(1),"r;;"],
% % % % % outliers2_y, outliers2_x, [symbol(2),"r;;"]);
end
if nargout
sout=s;
end
% % % endfunction
\ No newline at end of file
function y = myprctilecol(x,p)
xx = sort(x);
[m,n] = size(x);
if m==1 | n==1
m = max(m,n);
if m == 1
y = x*ones(length(p),1);
return
end
n = 1;
q = 100*(0.5:m - 0.5)./m;
xx = [min(x); xx(:); max(x)];
else
q = 100*(0.5:m - 0.5)./m;
xx = [min(x); xx; max(x)];
end
q = [0 q 100];
y = interp1(q,xx,p);
\ No newline at end of file
function h = cumplot(x)
%function h =cumplot(x)
% Inputs:
% - x [double] data series
%
% Outputs:
% - h [handle] figure handle
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
......@@ -26,9 +31,5 @@ function h = cumplot(x)
n=length(x);
x=[-inf; sort(x); Inf];
y=[0:n n]./n;
h0 = stairs(x,y);
grid on,
if nargout
h=h0;
end
h = stairs(x,y);
grid on
\ No newline at end of file
function [yy, xdir, isig, lam]=log_trans_(y0,xdir0,isig,lam)
function [yy, xdir, isig, lam]=log_transform(y0,xdir0,isig,lam)
% [yy, xdir, isig, lam]=log_transform(y0,xdir0,isig,lam)
% Conduct automatic log transformation lam(yy/isig+lam)
% Inputs:
% - y0 [double] series to transform
% - xdir [char] string indating the type of transformation:
% - log: standard log transformation
% - minuslog: log of minus (y0)
% - logsquared: log of y0^2
% - logskew: log of y0 shifted by lam
% - isig [double] scaling factor for y0
% - lam [double] shifting for y0
%
% Outputs:
% - yy [double] transformed series
% - xdir [char] string indating the type of transformation:
% - log: standard log transformation
% - minuslog: log of minus (y0)
% - logsquared: log of y0^2
% - logskew: log of y0 shifted by lam
% - isig [double] scaling factor for y0
% - lam [double] shifting for y0
%
% Notes: takes either one or four arguments. For one argument, the log
% transformation is conducted. For four arguments, the inverse
% transformation is applied.
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
......@@ -31,14 +56,18 @@ end
if nargin==1
xdir0='';
end
f=@(lam,y)gsa_skewness(log(y+lam));
f=@(lam,y)gsa.skewness(log(y+lam));
isig=1;
if ~(max(y0)<0 | min(y0)>0)
if gsa_skewness(y0)<0,
if ~(max(y0)<0 || min(y0)>0)
if gsa.skewness(y0)<0
isig=-1;
y0=-y0;
end
if isoctave
n=hist(y0,10);
else
n=histcounts(y0,10);
end
if n(1)>20*n(end)
try
lam=fzero(f,[-min(y0)+10*eps -min(y0)+abs(median(y0))],[],y0);
......@@ -48,7 +77,7 @@ if ~(max(y0)<0 | min(y0)>0)
if abs(yl(1))<abs(yl(2))
lam=-min(y0)+eps;
else
lam = -min(y0)+abs(median(y0)); %abs(100*(1+min(y0)));
lam = -min(y0)+abs(median(y0));
end
end
yy = log(y0+lam);
......@@ -63,10 +92,8 @@ else
if max(y0)<0
isig=-1;
y0=-y0;
%yy=log(-y0);
xdir=[xdir0,'_minuslog'];
elseif min(y0)>0
%yy=log(y0);
xdir=[xdir0,'_log'];
end
try
......@@ -77,7 +104,7 @@ else
if abs(yl(1))<abs(yl(2))
lam=-min(y0)+eps;
else
lam = -min(y0)+abs(median(y0)); %abs(100*(1+min(y0)));
lam = -min(y0)+abs(median(y0));
end
end
lam = max(lam,0);
......
function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters, BayesInfo)
function map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_, bayestopt_)
% map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_, bayestopt_)
% Inputs:
% - OutputDirectoryName [string] name of the output directory
% - M_ [structure] describing the model
% - options_ [structure] describing the options
% - oo_ [structure] storing the results
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2014-2016 European Commission
% Copyright © 2014-2018 Dynare Team
% Copyright © 2014-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -22,33 +30,32 @@ function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResult
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
fname_ = Model.fname;
fname_ = M_.fname;
np = EstimatedParameters.np;
nshock = EstimatedParameters.nvx + EstimatedParameters.nvn + EstimatedParameters.ncx + EstimatedParameters.ncn;
np = estim_params_.np;
nshock = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
pnames=cell(np,1);
pnames_tex=cell(np,1);
for jj=1:np
if DynareOptions.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj, DynareOptions.TeX, Model, EstimatedParameters, DynareOptions);
pnames_tex{jj,1} = strrep(param_name_tex_temp,'$','');
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj, options_.TeX, M_, estim_params_, options_.varobs);
pnames_tex{jj,1} = param_name_tex_temp;
pnames{jj,1} = param_name_temp;
else
param_name_temp = get_the_name(nshock+jj, DynareOptions.TeX, Model, EstimatedParameters, DynareOptions);
param_name_temp = get_the_name(nshock+jj, options_.TeX, M_, estim_params_, options_.varobs);
pnames{jj,1} = param_name_temp;
end
end
pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
indx_irf = [];
indx_moment = [];
init = ~DynareOptions.opt_gsa.load_stab;
init = ~options_.opt_gsa.load_stab;
options_mcf.pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
options_mcf.pvalue_corr = DynareOptions.opt_gsa.pvalue_corr;
options_mcf.alpha2 = DynareOptions.opt_gsa.alpha2_stab;
options_mcf.pvalue_ks = options_.opt_gsa.pvalue_ks;
options_mcf.pvalue_corr = options_.opt_gsa.pvalue_corr;
options_mcf.alpha2 = options_.opt_gsa.alpha2_stab;
options_mcf.param_names = pnames;
if DynareOptions.TeX
if options_.TeX
options_mcf.param_names_tex = pnames_tex;
end
options_mcf.fname_ = fname_;
......@@ -57,24 +64,24 @@ options_mcf.OutputDirectoryName = OutputDirectoryName;
skipline()
disp('Sensitivity analysis for calibration criteria')
if DynareOptions.opt_gsa.ppost
filetoload=dir([Model.dname filesep 'metropolis' filesep fname_ '_param_irf*.mat']);
if options_.opt_gsa.ppost
filetoload=dir([M_.dname filesep 'metropolis' filesep fname_ '_param_irf*.mat']);
lpmat=[];
for j=1:length(filetoload)
load([Model.dname filesep 'metropolis' filesep fname_ '_param_irf',int2str(j),'.mat'])
load([M_.dname filesep 'metropolis' filesep fname_ '_param_irf',int2str(j),'.mat'],'stock')
lpmat = [lpmat; stock];
clear stock
end
type = 'post';
else
if DynareOptions.opt_gsa.pprior
if options_.opt_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior'];
load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
load(filetoload,'lpmat','lpmat0')
lpmat = [lpmat0 lpmat];
type = 'prior';
else
filetoload=[OutputDirectoryName '/' fname_ '_mc'];
load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
load(filetoload,'lpmat','lpmat0')
lpmat = [lpmat0 lpmat];
type = 'mc';
end
......@@ -83,31 +90,31 @@ end
npar = size(pnames,1);
nshock = np - npar;
nbr_irf_restrictions = size(DynareOptions.endogenous_prior_restrictions.irf,1);
nbr_moment_restrictions = size(DynareOptions.endogenous_prior_restrictions.moment,1);
nbr_irf_restrictions = size(options_.endogenous_prior_restrictions.irf,1);
nbr_moment_restrictions = size(options_.endogenous_prior_restrictions.moment,1);
if init
mat_irf=cell(nbr_irf_restrictions,1);
for ij=1:nbr_irf_restrictions
mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3}));
mat_irf{ij}=NaN(Nsam,length(options_.endogenous_prior_restrictions.irf{ij,3}));
end
mat_moment=cell(nbr_moment_restrictions,1);
for ij=1:nbr_moment_restrictions
mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3}));
mat_moment{ij}=NaN(Nsam,length(options_.endogenous_prior_restrictions.moment{ij,3}));
end
irestrictions = [1:Nsam];
irestrictions = 1:Nsam;
h = dyn_waitbar(0,'Please wait...');
for j=1:Nsam
Model = set_all_parameters(lpmat(j,:)',EstimatedParameters,Model);
M_ = set_all_parameters(lpmat(j,:)',estim_params_,M_);
if nbr_moment_restrictions
[Tt,Rr,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state);
[Tt,Rr,~,info,oo_.dr, M_.params] = dynare_resolve(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
else
[Tt,Rr,SteadyState,info,DynareResults.dr, Model.params] = dynare_resolve(Model,DynareOptions,DynareResults.dr, DynareResults.steady_state, DynareResults.exo_steady_state, DynareResults.exo_det_steady_state,'restrict');
[Tt,Rr,~,info,oo_.dr, M_.params] = dynare_resolve(M_,options_,oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state,'restrict');
end
if info(1)==0
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
[~, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if ~isempty(info_irf)
for ij=1:nbr_irf_restrictions
mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
......@@ -123,14 +130,16 @@ if init
else
irestrictions(j)=0;
end
if mod(j,3)==0
dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
end
end
dyn_waitbar_close(h);
irestrictions=irestrictions(find(irestrictions));
xmat=lpmat(irestrictions,:);
skipline()
endo_prior_restrictions=DynareOptions.endogenous_prior_restrictions;
endo_prior_restrictions=options_.endogenous_prior_restrictions;
save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
else
load([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
......@@ -181,7 +190,7 @@ if ~isempty(indx_irf)
maxijv=0;
for ij=1:nbr_irf_restrictions
if length(endo_prior_restrictions.irf{ij,3})>maxijv
maxij=ij;maxijv=length(endo_prior_restrictions.irf{ij,3});
maxijv=length(endo_prior_restrictions.irf{ij,3});
end
plot_indx(ij) = find(strcmp(irf_couples,all_irf_couples(ij,:)));
time_matrix{plot_indx(ij)} = [time_matrix{plot_indx(ij)} endo_prior_restrictions.irf{ij,3}];
......@@ -189,8 +198,8 @@ if ~isempty(indx_irf)
iplot_indx = ones(size(plot_indx));
indx_irf = indx_irf(irestrictions,:);
if ~DynareOptions.nograph
h1=dyn_figure(DynareOptions.nodisplay,'name',[type ' evaluation of irf restrictions']);
if ~options_.nograph
h1=dyn_figure(options_.nodisplay,'name',[type ' evaluation of irf restrictions']);
nrow=ceil(sqrt(nbr_irf_couples));
ncol=nrow;
if nrow*(nrow-1)>nbr_irf_couples
......@@ -203,7 +212,7 @@ if ~isempty(indx_irf)
indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij);
for ik=1:size(mat_irf{ij},2)
[Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_irf{ij}(:,ik),0,DynareOptions.mh_conf_sig);
posterior_moments(mat_irf{ij}(:,ik),options_.mh_conf_sig);
irf_mean{plot_indx(ij)} = [irf_mean{plot_indx(ij)}; Mean];
irf_median{plot_indx(ij)} = [irf_median{plot_indx(ij)}; Median];
irf_var{plot_indx(ij)} = [irf_var{plot_indx(ij)}; Var];
......@@ -217,10 +226,10 @@ if ~isempty(indx_irf)
aleg = [aleg,'-' ,num2str(endo_prior_restrictions.irf{ij,3}(end))];
iplot_indx(ij)=0;
end
if ~DynareOptions.nograph && length(time_matrix{plot_indx(ij)})==1
if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1
set(0,'currentfigure',h1),
subplot(nrow,ncol, plot_indx(ij)),
hc = cumplot(mat_irf{ij}(:,ik));
hc = gsa.cumplot(mat_irf{ij}(:,ik));
a=axis;
delete(hc);
x1val=max(endo_prior_restrictions.irf{ij,4}(1),a(1));
......@@ -228,47 +237,34 @@ if ~isempty(indx_irf)
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
hold all,
set(hp,'FaceColor', [0.7 0.8 1])
hc = cumplot(mat_irf{ij}(:,ik));
hc = gsa.cumplot(mat_irf{ij}(:,ik));
set(hc,'color','k','linewidth',2)
hold off,
% hold off,
title([endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'],'interpreter','none'),
%set(legend_h,'Xlim',[0 1]);
% if ij==maxij
% leg1 = num2str(endo_prior_restrictions.irf{ij,3}(:));
% [legend_h,object_h,plot_h,text_strings]=legend(leg1);
% Position=get(legend_h,'Position');Position(1:2)=[-0.055 0.95-Position(4)];
% set(legend_h,'Position',Position);
% end
end
% hc = get(h,'Children');
%for i=2:2:length(hc)
%end
end
indx1 = find(indx_irf(:,ij)==0);
indx2 = find(indx_irf(:,ij)~=0);
atitle0=[endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the ',type,' support matches IRF ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.irf{ij,4})
% aname=[type '_irf_calib_',int2str(ij)];
aname=[type '_irf_calib_',endo_prior_restrictions.irf{ij,1},'_vs_',endo_prior_restrictions.irf{ij,2},'_',aleg];
atitle=[type ' IRF Calib: Parameter(s) driving ',endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'IRF restriction';
options_mcf.nobeha_title = 'NO IRF restriction';
if options_.TeX
options_mcf.beha_title_latex = 'IRF restriction';
options_mcf.nobeha_title_latex = 'NO IRF restriction';
end
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions);
gsa.monte_carlo_filtering_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
% indplot=find(proba<pvalue_ks);
% if ~isempty(indplot)
% stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
% end
end
for ij=1:nbr_irf_couples
if length(time_matrix{ij})>1
if ~DynareOptions.nograph
if ~options_.nograph
set(0,'currentfigure',h1);
subplot(nrow,ncol, ij)
itmp = (find(plot_indx==ij));
......@@ -282,7 +278,6 @@ if ~isempty(indx_irf)
tmp(temp_index,:) = endo_prior_restrictions.irf{itmp(ir),4};
end
end
% tmp = cell2mat(endo_prior_restrictions.irf(itmp,4));
tmp(isinf(tmp(:,1)),1)=a(3);
tmp(isinf(tmp(:,2)),2)=a(4);
hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],[tmp(:,1); tmp(end:-1:1,2)],'c');
......@@ -295,7 +290,6 @@ if ~isempty(indx_irf)
hold off
axis([max(1,a(1)) a(2:4)])
box on
%set(gca,'xtick',sort(time_matrix{ij}))
itmp = min(itmp);
title([endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}],'interpreter','none'),
end
......@@ -309,23 +303,27 @@ if ~isempty(indx_irf)
aleg = 'ALL';
atitle0=[endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}, '(', leg,')'];
fprintf(['%4.1f%% of the ',type,' support matches IRF restrictions ',atitle0,'\n'],length(indx1)/length(irestrictions)*100)
% aname=[type '_irf_calib_',int2str(ij)];
aname=[type '_irf_calib_',endo_prior_restrictions.irf{itmp,1},'_vs_',endo_prior_restrictions.irf{itmp,2},'_',aleg];
atitle=[type ' IRF Calib: Parameter(s) driving ',endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'IRF restriction';
options_mcf.nobeha_title = 'NO IRF restriction';
if options_.TeX
options_mcf.beha_title_latex = 'IRF restriction';
options_mcf.nobeha_title_latex = 'NO IRF restriction';
end
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions);
gsa.monte_carlo_filtering_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
end
end
end
if ~DynareOptions.nograph
dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],DynareOptions.nodisplay,DynareOptions.graph_format);
create_TeX_loader(DynareOptions,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],[type ' evaluation of irf restrictions'],'irf_restrictions',type,DynareOptions.figures.textwidth*min(ij/ncol,1))
if ~options_.nograph
dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],[type ' evaluation of irf restrictions'],'irf_restrictions',type,options_.figures.textwidth*min(ij/ncol,1))
end
skipline()
end
......@@ -361,24 +359,24 @@ if ~isempty(indx_moment)
skipline()
%get parameter names including standard deviations
np=size(BayesInfo.name,1);
np=size(bayestopt_.name,1);
name=cell(np,1);
name_tex=cell(np,1);
for jj=1:np
if DynareOptions.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(jj,DynareOptions.TeX,Model,EstimatedParameters,DynareOptions);
name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_.varobs);
name_tex{jj,1} = param_name_tex_temp;
name{jj,1} = param_name_temp;
else
param_name_temp = get_the_name(jj,DynareOptions.TeX,Model,EstimatedParameters,DynareOptions);
param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_.varobs);
name{jj,1} = param_name_temp;
end
end
options_mcf.param_names = name;
if DynareOptions.TeX
if options_.TeX
options_mcf.param_names_tex = name_tex;
end
options_mcf.param_names = BayesInfo.name;
options_mcf.param_names = bayestopt_.name;
all_moment_couples = cellstr([char(endo_prior_restrictions.moment(:,1)) char(endo_prior_restrictions.moment(:,2))]);
moment_couples = unique(all_moment_couples);
nbr_moment_couples = size(moment_couples,1);
......@@ -396,7 +394,7 @@ if ~isempty(indx_moment)
for ij=1:nbr_moment_restrictions
endo_prior_restrictions.moment{ij,3} = sort(endo_prior_restrictions.moment{ij,3});
if length(endo_prior_restrictions.moment{ij,3})>maxijv
maxij=ij;maxijv=length(endo_prior_restrictions.moment{ij,3});
maxijv=length(endo_prior_restrictions.moment{ij,3});
end
plot_indx(ij) = find(strcmp(moment_couples,all_moment_couples(ij,:)));
time_matrix{plot_indx(ij)} = [time_matrix{plot_indx(ij)} endo_prior_restrictions.moment{ij,3}];
......@@ -404,8 +402,8 @@ if ~isempty(indx_moment)
iplot_indx = ones(size(plot_indx));
indx_moment = indx_moment(irestrictions,:);
if ~DynareOptions.nograph
h2=dyn_figure(DynareOptions.nodisplay,'name',[type ' evaluation of moment restrictions']);
if ~options_.nograph
h2=dyn_figure(options_.nodisplay,'name',[type ' evaluation of moment restrictions']);
nrow=ceil(sqrt(nbr_moment_couples));
ncol=nrow;
if nrow*(nrow-1)>nbr_moment_couples
......@@ -419,7 +417,7 @@ if ~isempty(indx_moment)
indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij);
for ik=1:size(mat_moment{ij},2)
[Mean,Median,Var,HPD,Distrib] = ...
posterior_moments(mat_moment{ij}(:,ik),0,DynareOptions.mh_conf_sig);
posterior_moments(mat_moment{ij}(:,ik),options_.mh_conf_sig);
moment_mean{plot_indx(ij)} = [moment_mean{plot_indx(ij)}; Mean];
moment_median{plot_indx(ij)} = [moment_median{plot_indx(ij)}; Median];
moment_var{plot_indx(ij)} = [moment_var{plot_indx(ij)}; Var];
......@@ -433,10 +431,10 @@ if ~isempty(indx_moment)
aleg = [aleg,'_' ,num2str(endo_prior_restrictions.moment{ij,3}(end))];
iplot_indx(ij)=0;
end
if ~DynareOptions.nograph && length(time_matrix{plot_indx(ij)})==1
if ~options_.nograph && length(time_matrix{plot_indx(ij)})==1
set(0,'currentfigure',h2);
subplot(nrow,ncol,plot_indx(ij)),
hc = cumplot(mat_moment{ij}(:,ik));
hc = gsa.cumplot(mat_moment{ij}(:,ik));
a=axis; delete(hc),
% hist(mat_moment{ij}),
x1val=max(endo_prior_restrictions.moment{ij,4}(1),a(1));
......@@ -444,43 +442,34 @@ if ~isempty(indx_moment)
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
set(hp,'FaceColor', [0.7 0.8 1])
hold all
hc = cumplot(mat_moment{ij}(:,ik));
hc = gsa.cumplot(mat_moment{ij}(:,ik));
set(hc,'color','k','linewidth',2)
hold off
title([endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
% if ij==maxij
% leg1 = num2str(endo_prior_restrictions.moment{ij,3}(:));
% [legend_h,object_h,plot_h,text_strings]=legend(leg1);
% Position=get(legend_h,'Position');Position(1:2)=[-0.055 0.95-Position(4)];
% set(legend_h,'Position',Position);
% end
end
indx1 = find(indx_moment(:,ij)==0);
indx2 = find(indx_moment(:,ij)~=0);
atitle0=[endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the ',type,' support matches MOMENT ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.moment{ij,4})
% aname=[type '_moment_calib_',int2str(ij)];
aname=[type '_moment_calib_',endo_prior_restrictions.moment{ij,1},'_vs_',endo_prior_restrictions.moment{ij,2},'_',aleg];
atitle=[type ' MOMENT Calib: Parameter(s) driving ',endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'moment restriction';
options_mcf.nobeha_title = 'NO moment restriction';
if options_.TeX
options_mcf.beha_title_latex = 'moment restriction';
options_mcf.nobeha_title_latex = 'NO moment restriction';
end
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions);
gsa.monte_carlo_filtering_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
% [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
% indplot=find(proba<pvalue_ks);
% if ~isempty(indplot)
% stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
% end
end
for ij=1:nbr_moment_couples
time_matrix{ij} = sort(time_matrix{ij});
if length(time_matrix{ij})>1
if ~DynareOptions.nograph
if ~options_.nograph
itmp = (find(plot_indx==ij));
set(0,'currentfigure',h2);
subplot(nrow,ncol, ij)
......@@ -494,7 +483,6 @@ if ~isempty(indx_moment)
tmp(temp_index,:) = endo_prior_restrictions.moment{itmp(ir),4};
end
end
% tmp = cell2mat(endo_prior_restrictions.moment(itmp,4));
tmp(isinf(tmp(:,1)),1)=a(3);
tmp(isinf(tmp(:,2)),2)=a(4);
hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],[tmp(:,1); tmp(end:-1:1,2)],'b');
......@@ -507,7 +495,6 @@ if ~isempty(indx_moment)
hold off
axis(a)
box on
% set(gca,'xtick',sort(time_matrix{ij}))
itmp = min(itmp);
title([endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}],'interpreter','none'),
end
......@@ -521,23 +508,26 @@ if ~isempty(indx_moment)
aleg = 'ALL';
atitle0=[endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}, '(', leg,')'];
fprintf(['%4.1f%% of the ',type,' support matches MOMENT restrictions ',atitle0,'\n'],length(indx1)/length(irestrictions)*100)
% aname=[type '_moment_calib_',int2str(ij)];
aname=[type '_moment_calib_',endo_prior_restrictions.moment{itmp,1},'_vs_',endo_prior_restrictions.moment{itmp,2},'_',aleg];
atitle=[type ' MOMENT Calib: Parameter(s) driving ',endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}, '(', leg,')'];
options_mcf.amcf_name = aname;
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'moment restriction';
options_mcf.nobeha_title = 'NO moment restriction';
if options_.TeX
options_mcf.beha_title_latex = 'moment restriction';
options_mcf.nobeha_title_latex = 'NO moment restriction';
end
options_mcf.title = atitle0;
if ~isempty(indx1) && ~isempty(indx2)
mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions);
gsa.monte_carlo_filtering_analysis(xmat, indx1, indx2, options_mcf, M_, options_, bayestopt_, estim_params_);
end
end
end
end
if ~DynareOptions.nograph
dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],DynareOptions.nodisplay,DynareOptions.graph_format);
create_TeX_loader(DynareOptions,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],[type ' evaluation of moment restrictions'],'moment_restrictions',type,DynareOptions.figures.textwidth*min(ij/ncol,1))
if ~options_.nograph
dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],[type ' evaluation of moment restrictions'],'moment_restrictions',type,options_.figures.textwidth*min(ij/ncol,1))
end
skipline()
......
function map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
% map_identification(OutputDirectoryName,opt_gsa,M_,oo_,options_,estim_params_,bayestopt_)
% Inputs
% - OutputDirectoryName [string] name of the output directory
% - opt_gsa [structure] GSA options structure
% - M_ [structure] Matlab's structure describing the model
% - oo_ [structure] Matlab's structure describing the results
% - options_ [structure] Matlab's structure describing the current options
% - estim_params_ [structure] characterizing parameters to be estimated
% - bayestopt_ [structure] describing the priors
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2023 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 <https://www.gnu.org/licenses/>.
fname_ = M_.fname;
dr=oo_.dr;
nliv = opt_gsa.morris_nliv;
itrans = opt_gsa.trans_ident;
np = size(estim_params_.param_vals,1);
pnames = M_.param_names(estim_params_.param_vals(:,1));
if opt_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior'];
else
filetoload=[OutputDirectoryName '/' fname_ '_mc'];
end
load(filetoload,'lpmat','lpmat0','istable','T','yys')
if ~isempty(lpmat0)
lpmatx=lpmat0(istable,:);
else
lpmatx=[];
end
Nsam = size(lpmat,1);
nshock = size(lpmat0,2);
npT = np+nshock;
fname_ = M_.fname;
if opt_gsa.load_ident_files==0
mss = yys(bayestopt_.mfys,:);
mss = gsa.teff(mss(:,istable),Nsam,istable);
yys = gsa.teff(yys(dr.order_var,istable),Nsam,istable);
if exist('T','var')
[vdec, cc, ac] = gsa.monte_carlo_moments(T, lpmatx, dr, M_, options_, estim_params_);
else
return
end
if opt_gsa.morris==2
pdraws = identification.run(M_,oo_,options_,bayestopt_,estim_params_,options_.options_ident,[lpmatx lpmat(istable,:)]);
if ~isempty(pdraws) && max(max(abs(pdraws-[lpmatx lpmat(istable,:)])))==0
disp(['Sample check OK. Largest difference: ', num2str(max(max(abs(pdraws-[lpmatx lpmat(istable,:)]))))]),
clear pdraws;
end
clear GAM gas
end
if opt_gsa.morris~=1 && M_.exo_nbr>1
ifig=0;
for j=1:M_.exo_nbr
if mod(j,6)==1
hh_fig=dyn_figure(options_.nodisplay,'name','Variance decomposition shocks');
ifig=ifig+1;
iplo=0;
end
iplo=iplo+1;
subplot(2,3,iplo)
gsa.boxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:size(options_.varobs,1))
set(gca,'xlim',[0.5 size(options_.varobs,1)+0.5])
set(gca,'ylim',[-2 102])
for ip=1:size(options_.varobs,1)
if options_.TeX
text(ip,-4,deblank(opt_gsa.varobs_tex(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-4,deblank(options_.varobs(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
ylabel(' ')
title(M_.exo_names{j},'interpreter','none')
if mod(j,6)==0 || j==M_.exo_nbr
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],ifig,'Variance decomposition shocks','vdec_exo',options_.figures.textwidth*min(iplo/3,1))
end
end
end
for j=1:size(cc,1)
cc(j,j,:)=gsa.standardize_columns(squeeze(log(cc(j,j,:))))./2;
end
[vdec, ~, ir_vdec, ic_vdec] = gsa.teff(vdec,Nsam,istable);
[cc, ~, ir_cc, ic_cc] = gsa.teff(cc,Nsam,istable);
[ac, ~, ir_ac, ic_ac] = gsa.teff(ac,Nsam,istable);
nc1= size(T,2);
endo_nbr = M_.endo_nbr;
nstatic = M_.nstatic;
nspred = M_.nspred;
iv = (1:endo_nbr)';
ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
dr.ghx = T(:, 1:(nc1-M_.exo_nbr),1);
dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, 1);
[Aa,Bb] = kalman_transition_matrix(dr,iv,ic);
A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
if ~isempty(lpmatx)
M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(1,:));
end
A(:,:,1)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
for j=2:length(istable)
dr.ghx = T(:, 1:(nc1-M_.exo_nbr),j);
dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, j);
[Aa,Bb] = kalman_transition_matrix(dr, iv, ic);
if ~isempty(lpmatx)
M_=gsa.set_shocks_param(M_,estim_params_,lpmatx(j,:));
end
A(:,:,j)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
end
clear T
clear lpmatx
[yt, j0]=gsa.teff(A,Nsam,istable);
yt = [yys yt];
if opt_gsa.morris==2
clear TAU A
else
clear A
end
save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'ac','cc','vdec','yt','mss')
else %load identification files
load([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'ac','cc','vdec','yt','mss')
end
if opt_gsa.morris==1
if ~isempty(vdec)
if opt_gsa.load_ident_files==0
SAMorris=NaN(npT,3,size(vdec,2));
for i=1:size(vdec,2)
[~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], vdec(:,i),nliv);
end
SAvdec = squeeze(SAMorris(:,1,:))';
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec')
end
hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition');
gsa.boxplot(SAvdec,[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
ydum = get(gca,'ylim');
set(gca,'ylim',[0 ydum(2)])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:npT
if options_.TeX
[~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_.varobs);
text(ip,-2,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-2,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
title('Elementary effects variance decomposition')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morris_vdec'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_morris_vdec'],1,'Screening identification: variance decomposition','morris_vdec',1)
else
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'vdec')
end
if opt_gsa.load_ident_files==0
ccac = [mss cc ac];
SAMorris=NaN(npT,3,size(ccac,2));
for i=1:size(ccac,2)
[~, SAMorris(:,:,i)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], [ccac(:,i)],nliv);
end
SAcc = squeeze(SAMorris(:,1,:))';
SAcc = SAcc./(max(SAcc,[],2)*ones(1,npT));
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAcc','cc','ir_cc','ic_cc','-append')
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'ac','ir_ac','ic_ac','-append')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAcc','cc','ir_cc','ic_cc')
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'ac','ir_ac','ic_ac')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: theoretical moments');
gsa.boxplot(SAcc,[],'.',[],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])
for ip=1:npT
if options_.TeX
[~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_.varobs);
text(ip,-0.02,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
title('Elementary effects in the moments')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morris_moments'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_morris_moments'],1,'Screening identification: theoretical moments','morris_moments',1)
if opt_gsa.load_ident_files==0
SAMorris=NaN(npT,3,j0);
for j=1:j0
[~, SAMorris(:,:,j)] = gsa.Morris_Measure_Groups(npT, [lpmat0 lpmat], yt(:,j),nliv);
end
SAM = squeeze(SAMorris(1:end,1,:));
SAnorm=NaN(npT,j0);
irex=NaN(j0);
for j=1:j0
SAnorm(:,j)=SAM(:,j)./max(SAM(:,j));
irex(j)=length(find(SAnorm(:,j)>0.01));
end
SAMmu = squeeze(SAMorris(1:end,2,:));
SAmunorm=NaN(npT,j0);
for j=1:j0
SAmunorm(:,j)=SAMmu(:,j)./max(SAM(:,j)); % normalised w.r.t. mu*
end
SAMsig = squeeze(SAMorris(1:end,3,:));
SAsignorm=NaN(npT,j0);
for j=1:j0
SAsignorm(:,j)=SAMsig(:,j)./max(SAMsig(:,j));
end
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAnorm','SAmunorm','SAsignorm','-append')
else
load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model');
gsa.boxplot(SAnorm',[],'.',[],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
if options_.TeX
[~, param_name_tex_temp]= get_the_name(ip,options_.TeX,M_,estim_params_,options_.varobs);
text(ip,-0.02,param_name_tex_temp,'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,bayestopt_.name{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
title('Elementary effects in the model')
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_morris_par'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_morris_par'],1,'Screening identification: model','morris_par',1)
elseif opt_gsa.morris==3
return
elseif opt_gsa.morris==2 % ISKREV stuff
return
else
error('gsa/map_identification: unsupported option morris=%u',opt_gsa.morris)
end
function []=create_TeX_loader(options_,figpath,ifig_number,caption,label_name,scale_factor)
if nargin<6
scale_factor=1;
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([figpath '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by map_ident_.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s}\n',scale_factor,strrep(figpath,'\','/'));
fprintf(fidTeX,'\\caption{%s.}',caption);
fprintf(fidTeX,'\\label{Fig:%s:%u}\n',label_name,ifig_number);
fprintf(fidTeX,'\\end{figure}\n\n');
fprintf(fidTeX,'%% End Of TeX file. \n');
fclose(fidTeX);
end
function yr = trank(y)
% yr is the rank transformation of y
yr=NaN(size(y));
[nr, nc] = size(y);
for j=1:nc
[~, is]=sort(y(:,j));
yr(is,j)=[1:nr]'./nr;
end
function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% function [rmse_MC, ixx] = filt_mc_(OutDir)
function [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_)
% [rmse_MC, ixx] = monte_carlo_filtering(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_
% Inputs:
% - OutputDirectoryName [string] name of the output directory
% - options_gsa_ [structure] GSA options
% - dataset_ [dseries] object storing the dataset
% - dataset_info [structure] storing informations about the sample.
% - M_ [structure] Matlab's structure describing the model
% - oo_ [structure] storing the results
% - options_ [structure] Matlab's structure describing the current options
% - bayestopt_ [structure] describing the priors
% - estim_params_ [structure] characterizing parameters to be estimated
%
% Outputs:
% - rmse_MC [double] RMSE by nvar matrix of the RMSEs
% - ixx [double] RMSE by nvar matrix of sorting
% indices (descending order of RMSEs)
%
% Notes: the R^2 definition is 1-var(ymodel-ydata)/var(ydata). It ranges
% bewteen (-inf, 1], with negative values indicating that themodel is a worse
% predictor than the sample mean of the data
% inputs (from opt_gsa structure)
% vvarvecm = options_gsa_.var_rmse;
% loadSA = options_gsa_.load_rmse;
......@@ -7,7 +27,6 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% alpha = options_gsa_.alpha_rmse;
% alpha2 = options_gsa_.alpha2_rmse;
% istart = options_gsa_.istart_rmse;
% alphaPC = 0.5;
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
......@@ -31,9 +50,6 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global bayestopt_ estim_params_ M_ options_ oo_
% options_gsa_=options_.opt_gsa;
vvarvecm = options_gsa_.var_rmse;
if options_.TeX
vvarvecm_tex = options_gsa_.var_rmse_tex;
......@@ -46,13 +62,10 @@ alpha = options_gsa_.alpha_rmse;
alpha2 = 0;
pvalue = options_gsa_.alpha2_rmse;
istart = max(2,options_gsa_.istart_rmse);
alphaPC = 0.5;
fname_ = M_.fname;
lgy_ = M_.endo_names;
dr_ = oo_.dr;
skipline(2)
skipline(1)
disp('Starting sensitivity analysis')
disp('for the fit of EACH observed series ...')
skipline()
......@@ -61,12 +74,12 @@ if ~options_.nograph
a=dir([OutDir,filesep,'*.*']);
tmp1='0';
if options_.opt_gsa.ppost
tmp=['_rmse_post'];
tmp='_rmse_post';
else
if options_.opt_gsa.pprior
tmp=['_rmse_prior'];
tmp='_rmse_prior';
else
tmp=['_rmse_mc'];
tmp='_rmse_mc';
end
if options_gsa_.lik_only
tmp1 = [tmp,'_post_SA'];
......@@ -75,17 +88,23 @@ if ~options_.nograph
end
for j=1:length(a)
if strmatch([fname_,tmp],a(j).name)
if options_.debug
disp(a(j).name)
end
delete([OutDir,filesep,a(j).name])
end
if strmatch([fname_,tmp1],a(j).name)
if options_.debug
disp(a(j).name)
end
delete([OutDir,filesep,a(j).name])
end
end
disp('done !')
end
[param_names,param_names_tex]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
nshock=estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
npar=estim_params_.np;
if ~isempty(options_.mode_file)
......@@ -94,10 +113,12 @@ end
if options_.opt_gsa.ppost
c=load([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'xparam1');
xparam1_mean=c.xparam1;
xparam1=c.xparam1;
clear c
elseif ~isempty(options_.mode_file) && exist([M_.dname filesep 'Output' filesep fname_,'_mean.mat'])==2
elseif ~isempty(options_.mode_file) && exist([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'file')==2
c=load([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'xparam1');
xparam1_mean=c.xparam1;
xparam1=c.xparam1;
clear c
end
......@@ -124,31 +145,11 @@ if loadSA
end
end
if ~loadSA
if exist('xparam1','var')
M_ = set_all_parameters(xparam1,estim_params_,M_);
ys_mode=evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
end
if exist('xparam1_mean','var')
M_ = set_all_parameters(xparam1_mean,estim_params_,M_);
ys_mean=evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
end
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
load([DirectoryName filesep M_.fname '_data.mat']);
filfilt = dir([DirectoryName filesep M_.fname '_filter_step_ahead*.mat']);
temp_smooth_file_list = dir([DirectoryName filesep M_.fname '_smooth*.mat']);
jfile=0;
for j=1:length(temp_smooth_file_list)
if isempty(strfind(temp_smooth_file_list(j).name,'smoothed')),
jfile=jfile+1;
filsmooth(jfile)=temp_smooth_file_list(j);
end
end
filupdate = dir([DirectoryName filesep M_.fname '_update*.mat']);
filparam = dir([DirectoryName filesep M_.fname '_param*.mat']);
x=[];
......@@ -156,11 +157,11 @@ if ~loadSA
sto_ys=[];
for j=1:length(filparam)
if isempty(strmatch([M_.fname '_param_irf'],filparam(j).name))
load([DirectoryName filesep filparam(j).name]);
x=[x; stock];
logpo2=[logpo2; stock_logpo];
sto_ys=[sto_ys; stock_ys];
clear stock stock_logpo stock_ys;
temp=load([DirectoryName filesep filparam(j).name]); % from prior_posterior_statistics_core
x=[x; temp.stock];
logpo2=[logpo2; temp.stock_logpo];
sto_ys=[sto_ys; temp.stock_ys];
clear temp;
end
end
nruns=size(x,1);
......@@ -168,38 +169,41 @@ if ~loadSA
if options_.opt_gsa.ppost || (options_.opt_gsa.ppost==0 && options_.opt_gsa.lik_only==0)
skipline()
disp('Computing RMSE''s...')
jxj=NaN(length(vvarvecm),1);
js=NaN(length(vvarvecm),1);
yss=NaN(length(vvarvecm),gend,size(sto_ys,1));
for i = 1:length(vvarvecm)
vj = vvarvecm{i};
jxj(i) = strmatch(vj, lgy_(dr_.order_var), 'exact');
js(i) = strmatch(vj, lgy_, 'exact');
jxj(i) = strmatch(vj, M_.endo_names(oo_.dr.order_var), 'exact');
js(i) = strmatch(vj, M_.endo_names, 'exact');
yss(i,:,:)=repmat(sto_ys(:,js(i))',[gend,1]);
end
if exist('xparam1','var')
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);% + kron(ys_mode(js),ones(1,gend)));
yobs = transpose( ahat(jxj,:));% + kron(ys_mode(js),ones(1,gend)));
[~,~,~,ahat,~,~,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);
yobs = transpose( ahat(jxj,:));
rmse_mode = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2));
r2_mode = 1-sum((yobs(istart:end,:)-y0(istart:end,:)).^2)./sum(yobs(istart:end,:).^2);
end
y0=-yss;
y0=-yss; %demean everything using the theoretical mean, i.e. steady state
nbb=0;
for j=1:length(filfilt)
load([DirectoryName filesep M_.fname '_filter_step_ahead',num2str(j),'.mat']);
nb = size(stock,4);
y0(:,:,nbb+1:nbb+nb)=y0(:,:,nbb+1:nbb+nb)+reshape(stock(1,js,1:gend,:),[length(js) gend nb]);
temp=load([DirectoryName filesep M_.fname '_filter_step_ahead',num2str(j),'.mat']);
nb = size(temp.stock,4);
y0(:,:,nbb+1:nbb+nb)=y0(:,:,nbb+1:nbb+nb)+reshape(temp.stock(1,js,1:gend,:),[length(js) gend nb]);
nbb=nbb+nb;
clear stock;
clear temp;
end
yobs=-yss;
nbb=0;
for j=1:length(filupdate)
load([DirectoryName filesep M_.fname '_update',num2str(j),'.mat']);
nb = size(stock,3);
yobs(:,:,nbb+1:nbb+nb)=yobs(:,:,nbb+1:nbb+nb)+reshape(stock(js,1:gend,:),[length(js) gend nb]);
temp=load([DirectoryName filesep M_.fname '_update',num2str(j),'.mat']);
nb = size(temp.stock,3);
yobs(:,:,nbb+1:nbb+nb)=yobs(:,:,nbb+1:nbb+nb)+reshape(temp.stock(js,1:gend,:),[length(js) gend nb]);
nbb=nbb+nb;
clear stock;
clear temp;
end
y0M=mean(y0,2);
rmse_MC=zeros(nruns,length(js));
r2_MC=zeros(nruns,length(js));
for j=1:nruns
......@@ -207,14 +211,15 @@ if ~loadSA
r2_MC(j,:) = 1-mean((yobs(:,istart:end,j)'-y0(:,istart:end,j)').^2)./mean((yobs(:,istart:end,j)').^2);
end
if exist('xparam1_mean','var')
[alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);% + kron(ys_mean(js),ones(1,gend)));
yobs = transpose( ahat(jxj,:));% + kron(ys_mean(js),ones(1,gend)));
[~,~,~,ahat,~,~,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_);
y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);
yobs = transpose( ahat(jxj,:));
rmse_pmean = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2));
r2_pmean = 1-mean((yobs(istart:end,:)-y0(istart:end,:)).^2)./mean(yobs(istart:end,:).^2);
end
clear stock_filter;
end
lnprior=NaN(nruns,1);
for j=1:nruns
lnprior(j,1) = priordens(x(j,:)',bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
end
......@@ -242,7 +247,7 @@ if ~loadSA
end
end
end
else
else % loadSA
if options_.opt_gsa.lik_only && options_.opt_gsa.ppost==0
load([OutDir,filesep,fnamtmp, '.mat'],'x','logpo2','likelihood');
else
......@@ -252,27 +257,27 @@ else
nruns=size(x,1);
nfilt=floor(pfilt*nruns);
end
% smirnov tests
% Smirnov tests
nfilt0 = nfilt*ones(length(vvarvecm), 1);
logpo2=logpo2(:);
if ~options_.opt_gsa.ppost
[dum, ipost]=sort(-logpo2);
[dum, ilik]=sort(-likelihood);
[~, ipost]=sort(-logpo2);
[~, ilik]=sort(-likelihood);
end
% visual scatter analysis!
if options_.opt_gsa.ppost
tmp_title='R2 Posterior:';
atitle='R2 Posterior:';
tmp_title='R2 Scatter plot: Posterior';
atitle='R2 Scatter plot: Posterior';
asname='r2_post';
else
if options_.opt_gsa.pprior
tmp_title='R2 Prior:';
atitle='R2 Prior:';
tmp_title='R2 Scatter plot: Prior';
atitle='R2 Scatter plot: Prior';
asname='r2_prior';
else
tmp_title='R2 MC:';
atitle='R2 MC:';
tmp_title='R2 Scatter plot: MC';
atitle='R2 Scatter plot: MC';
asname='r2_mc';
end
end
......@@ -283,7 +288,7 @@ options_scatter.OutputDirectoryName = OutDir;
options_scatter.amcf_name = asname;
options_scatter.amcf_title = atitle;
options_scatter.title = tmp_title;
scatter_analysis(r2_MC, x,options_scatter, options_);
gsa.scatter_analysis(r2_MC, x,options_scatter, options_);
% end of visual scatter analysis
if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
......@@ -297,13 +302,10 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.pvalue_ks = alpha;
options_mcf.pvalue_corr = pvalue;
options_mcf.alpha2 = alpha2;
options_mcf.param_names = param_names;
if options_.TeX
[pnames,pnames_tex]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
options_mcf.param_names = pnames;
options_mcf.param_names_tex = pnames_tex;
options_mcf.param_names_tex = param_names_tex;
else
[pnames]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
options_mcf.param_names = pnames;
options_mcf.param_names_tex = {};
end
options_mcf.fname_ = fname_;
......@@ -313,7 +315,12 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.title = atitle;
options_mcf.beha_title = 'better posterior kernel';
options_mcf.nobeha_title = 'worse posterior kernel';
mcf_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, options_);
if options_.TeX
options_mcf.beha_title_latex = 'better posterior kernel';
options_mcf.nobeha_title_latex = 'worse posterior kernel';
end
gsa.monte_carlo_filtering_analysis(x, ipost(1:nfilt), ipost(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
if options_.opt_gsa.pprior
anam = 'rmse_prior_lik';
atitle = 'RMSE prior: Log Likelihood Kernel';
......@@ -326,15 +333,20 @@ if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
options_mcf.title = atitle;
options_mcf.beha_title = 'better likelihood';
options_mcf.nobeha_title = 'worse likelihood';
mcf_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, options_);
if options_.TeX
options_mcf.beha_title_latex = 'better likelihood';
options_mcf.nobeha_title_latex = 'worse likelihood';
end
gsa.monte_carlo_filtering_analysis(x, ilik(1:nfilt), ilik(nfilt+1:end), options_mcf, M_, options_, bayestopt_, estim_params_);
else
if options_.opt_gsa.ppost
rmse_txt=rmse_pmean;
r2_txt=r2_pmean;
else
if options_.opt_gsa.pprior || ~exist('rmse_pmean')
if exist('rmse_mode')
if options_.opt_gsa.pprior || ~exist('rmse_pmean','var')
if exist('rmse_mode','var')
rmse_txt=rmse_mode;
r2_txt=r2_mode;
else
......@@ -346,18 +358,19 @@ else
r2_txt=r2_pmean;
end
end
ixx=NaN(size(rmse_MC,1),length(vvarvecm));
for i = 1:length(vvarvecm)
[dum, ixx(:,i)] = sort(rmse_MC(:,i));
[~, ixx(:,i)] = sort(rmse_MC(:,i));
end
PP = ones(npar+nshock, length(vvarvecm));
PPV = ones(length(vvarvecm), length(vvarvecm), npar+nshock);
SS = zeros(npar+nshock, length(vvarvecm));
for j = 1:npar+nshock
for i = 1:length(vvarvecm)
[H, P, KSSTAT] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
[H1, P1, KSSTAT1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
[H2, P2, KSSTAT2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
if H1 & H2==0
[~, P] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha);
[H1] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1);
[H2] = gsa.smirnov_test(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1);
if H1==0 && H2==0
SS(j,i)=1;
elseif H1==0
SS(j,i)=-1;
......@@ -369,7 +382,7 @@ else
for i = 1:length(vvarvecm)
for l = 1:length(vvarvecm)
if l~=i && PP(j,i)<alpha && PP(j,l)<alpha
[H,P,KSSTAT] = smirnov(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
[~,P] = gsa.smirnov_test(x(ixx(1:nfilt0(i),i),j),x(ixx(1:nfilt0(l),l),j), alpha);
PPV(i,l,j) = P;
elseif l==i
PPV(i,l,j) = PP(j,i);
......@@ -394,13 +407,17 @@ else
hh_fig=dyn_figure(options_.nodisplay,'name',[temp_name,' ',int2str(ifig)]);
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(lnprior(ixx(1:nfilt0(i),i)));
h=gsa.cumplot(lnprior(ixx(1:nfilt0(i),i)));
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(lnprior);
hold on, h=gsa.cumplot(lnprior);
set(h,'color','k','linewidth',1)
h=cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
h=gsa.cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
set(h,'color','red','linewidth',2)
if options_.TeX
title(vvarvecm_tex{i},'interpreter','latex')
else
title(vvarvecm{i},'interpreter','none')
end
if mod(i,9)==0 || i==length(vvarvecm)
if ~isoctave
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Log-prior for BETTER R2','Color','Blue','horizontalalignment','center');
......@@ -442,13 +459,17 @@ else
hh_fig = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(likelihood(ixx(1:nfilt0(i),i)));
h=gsa.cumplot(likelihood(ixx(1:nfilt0(i),i)));
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(likelihood);
hold on, h=gsa.cumplot(likelihood);
set(h,'color','k','linewidth',1)
h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
h=gsa.cumplot(likelihood(ixx(nfilt0(i)+1:end,i)));
set(h,'color','red','linewidth',2)
if options_.TeX
title(vvarvecm_tex{i},'interpreter','latex')
else
title(vvarvecm{i},'interpreter','none')
end
if options_.opt_gsa.ppost==0
set(gca,'xlim',[min( likelihood(ixx(1:nfilt0(i),i)) ) max( likelihood(ixx(1:nfilt0(i),i)) )])
end
......@@ -493,13 +514,17 @@ else
hh_fig = dyn_figure(options_.nodisplay,'Name',[temp_name,' ',int2str(ifig)]);
end
subplot(3,3,i-9*(ifig-1))
h=cumplot(logpo2(ixx(1:nfilt0(i),i)));
h=gsa.cumplot(logpo2(ixx(1:nfilt0(i),i)));
set(h,'color','blue','linewidth',2)
hold on, h=cumplot(logpo2);
hold on, h=gsa.cumplot(logpo2);
set(h,'color','k','linewidth',1)
h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
h=gsa.cumplot(logpo2(ixx(nfilt0(i)+1:end,i)));
set(h,'color','red','linewidth',2)
if options_.TeX
title(vvarvecm_tex{i},'interpreter','latex')
else
title(vvarvecm{i},'interpreter','none')
end
if options_.opt_gsa.ppost==0
set(gca,'xlim',[min( logpo2(ixx(1:nfilt0(i),i)) ) max( logpo2(ixx(1:nfilt0(i),i)) )])
end
......@@ -529,15 +554,6 @@ else
end
end
end
if options_.TeX
[pnames,pnames_tex]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
param_names = pnames;
param_names_tex = pnames_tex;
else
[pnames]=get_LaTeX_parameter_names(M_,options_,estim_params_,bayestopt_);
param_names = pnames;
param_names_tex = {};
end
skipline()
title_string='RMSE over the MC sample:';
data_mat=[min(rmse_MC)' max(rmse_MC)'];
......@@ -549,7 +565,7 @@ else
end
invar = find( std(rmse_MC)./mean(rmse_MC)<=0.0001 );
if ~isempty(invar)
skipline(2)
skipline(1)
disp('RMSE is not varying significantly over the MC sample for the following variables:')
disp(vvarvecm{invar})
disp('These variables are excluded from SA')
......@@ -561,8 +577,7 @@ else
rmse_MC = rmse_MC(:,ivar);
skipline()
disp(['Sample filtered the ',num2str(pfilt*100),'% best RMSE''s for each observed series ...' ])
skipline(2)
disp('RMSE ranges after filtering:')
skipline(1)
title_string='RMSE ranges after filtering:';
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
headers = {'Variable'; 'min'; 'max'; 'min'; 'max'; 'posterior mode'};
......@@ -589,7 +604,7 @@ else
else
values_length = max(ceil(max(max(log10(abs(data_mat(isfinite(data_mat))))))),1)+val_precis+1;
end
if any(data_mat) < 0 %add one character for minus sign
if any(data_mat < 0) %add one character for minus sign
values_length = values_length+1;
end
headers_length = cellofchararraymaxlength(headers(2:end));
......@@ -598,7 +613,6 @@ else
else
val_width = max(headers_length, values_length)+2;
end
value_format = sprintf('%%%d.%df',val_width,val_precis);
header_string_format = sprintf('%%%ds',val_width);
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
optional_header=sprintf([label_format_leftbound,header_string_format,header_string_format,header_string_format,header_string_format],'','',['best ',num2str(pfilt*100),'% filtered'],'','remaining 90%');
......@@ -610,7 +624,7 @@ else
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
optional_header={[' & \multicolumn{2}{c}{best ',num2str(pfilt*100),' filtered} & \multicolumn{2}{c}{remaining 90\%}\\']};
else
optional_header={[' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\']};
optional_header={' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\'};
end
dyn_latex_table(M_, options_, title_string, 'RMSE_ranges_after_filtering', headers_tex, vvarvecm_tex, data_mat, 0, val_width, val_precis, optional_header);
end
......@@ -657,7 +671,7 @@ else
else
values_length = max(ceil(max(max(log10(abs(data_mat(isfinite(data_mat))))))),1)+val_precis+1;
end
if any(data_mat) < 0 %add one character for minus sign
if any(data_mat < 0) %add one character for minus sign
values_length = values_length+1;
end
headers_length = cellofchararraymaxlength(headers(2:end));
......@@ -666,7 +680,6 @@ else
else
val_width = max(headers_length, values_length)+2;
end
value_format = sprintf('%%%d.%df',val_width,val_precis);
header_string_format = sprintf('%%%ds',val_width);
if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior
......@@ -679,7 +692,7 @@ else
if ~options_.opt_gsa.ppost && options_.opt_gsa.pprior
optional_header = {[' & \multicolumn{2}{c}{best ',num2str(pfilt*100),' filtered} & \multicolumn{2}{c}{remaining 90\%}\\']};
else
optional_header = {[' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\']};
optional_header = {' & \multicolumn{2}{c}{best filtered} & \multicolumn{2}{c}{remaining}\\'};
end
dyn_latex_table(M_, options_, title_string, 'R2_ranges_after_filtering', headers_tex, vvarvecm_tex, data_mat, 0, val_width, val_precis, optional_header);
end
......@@ -690,16 +703,15 @@ else
SP(ns,j)=ones(size(ns));
SS(:,j)=SS(:,j).*SP(:,j);
end
for j=1:npar+nshock %estim_params_.np,
nsp=NaN(npar+nshock,1);
for j=1:npar+nshock
nsp(j)=length(find(SP(j,:)));
end
snam0=param_names(find(nsp==0));
snam1=param_names(find(nsp==1));
snam2=param_names(find(nsp>1));
snam=param_names(find(nsp>0));
snam0=param_names(nsp==0);
snam1=param_names(nsp==1);
snam2=param_names(nsp>1);
nsnam=(find(nsp>1));
skipline(2)
skipline(1)
disp('These parameters do not affect significantly the fit of ANY observed series:')
disp(char(snam0))
skipline()
......@@ -708,7 +720,6 @@ else
skipline()
disp('These parameters affect MORE THAN ONE observed series: trade off exists!')
disp(char(snam2))
pnam=bayestopt_.name;
% plot trade-offs
if ~options_.nograph
a00=jet(length(vvarvecm));
......@@ -740,8 +751,12 @@ else
options_mcf.amcf_title = [atitle ' ' vvarvecm{iy}];
options_mcf.beha_title = ['better fit of ' vvarvecm{iy}];
options_mcf.nobeha_title = ['worse fit of ' vvarvecm{iy}];
if options_.TeX
options_mcf.beha_title_latex = ['better fit of ' vvarvecm_tex{iy}];
options_mcf.nobeha_title_latex = ['worse fit of ' vvarvecm_tex{iy}];
end
options_mcf.title = ['the fit of ' vvarvecm{iy}];
mcf_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, options_);
gsa.monte_carlo_filtering_analysis(x, ixx(1:nfilt0(iy),iy), ixx(nfilt0(iy)+1:end,iy), options_mcf, M_, options_, bayestopt_, estim_params_);
end
for iy = 1:length(vvarvecm)
ipar = find(any(squeeze(PPV(iy,:,:))<alpha));
......@@ -749,53 +764,61 @@ else
hh_fig = dyn_figure(options_.nodisplay,'name',[temp_name,' observed variable ', vvarvecm{iy}]);
for j=1+5*(ix-1):min(length(ipar),5*ix)
subplot(2,3,j-5*(ix-1))
h0=cumplot(x(:,ipar(j)));
h0=gsa.cumplot(x(:,ipar(j)));
set(h0,'color',[0 0 0])
hold on,
iobs=find(squeeze(PPV(iy,:,ipar(j)))<alpha);
for i = 1:length(vvarvecm)
if any(iobs==i) || i==iy
h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),ipar(j)));
if ~isoctave
hcmenu = uicontextmenu;
uimenu(hcmenu,'Label',vvarvecm{i});
set(h0,'uicontextmenu',hcmenu)
end
else
h0=cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),ipar(j))*NaN);
end
set(h0,'color',a00(i,:),'linewidth',2)
end
ydum=get(gca,'ylim');
if exist('xparam1')
if exist('xparam1','var')
xdum=xparam1(ipar(j));
h1=plot([xdum xdum],ydum);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
end
xlabel('')
title([pnam{ipar(j)}],'interpreter','none')
if options_.TeX
title([param_names_tex{ipar(j)}],'interpreter','latex')
else
title([param_names{ipar(j)}],'interpreter','none')
end
end
if isoctave
legend(vertcat('base',vvarvecm),'location','eastoutside');
else
h0=legend(vertcat('base',vvarvecm));
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none');
if options_.TeX
h0=legend(vertcat('base',vvarvecm_tex),'interpreter','latex');
else
h0=legend(vertcat('base',vvarvecm),'interpreter','none');
end
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3]);
end
if options_.opt_gsa.ppost
dyn_saveas(hh_fig,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable $',vvarvecm_tex{iy} '$'],['rmse_post_' vvarvecm{iy}],1)
create_TeX_loader(options_,[ OutDir filesep fname_ '_rmse_post_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable ',vvarvecm_tex{iy} ],['rmse_post_' vvarvecm{iy}],1)
end
else
if options_.opt_gsa.pprior
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],ix,[temp_name,' observed variable $',vvarvecm_tex{iy} '$'],['rmse_prior_' vvarvecm{iy}],1)
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_prior_' vvarvecm{iy} '_' int2str(ix) ],ix,[temp_name,' observed variable ',vvarvecm_tex{iy}],['rmse_prior_' vvarvecm{iy}],1)
end
else
dyn_saveas(hh_fig,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],options_.nodisplay,options_.graph_format);
if options_.TeX
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable $',vvarvecm_tex{iy} '$'],['rmse_mc_' vvarvecm{iy}],1)
create_TeX_loader(options_,[OutDir filesep fname_ '_rmse_mc_' vvarvecm{iy} '_' int2str(ix)],ix,[temp_name,' observed variable ',vvarvecm_tex{iy}],['rmse_mc_' vvarvecm{iy}],1)
end
end
end
......@@ -806,15 +829,15 @@ else
hh_fig = dyn_figure(options_.nodisplay,'name',[temp_name,' estimated params and shocks ',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)));
h0=gsa.cumplot(x(:,nsnam(j)));
set(h0,'color',[0 0 0])
hold on,
npx=find(SP(nsnam(j),:)==0);
for i = 1:length(vvarvecm)
if any(npx==i)
h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
else
h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
h0=gsa.cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
if ~isoctave
hcmenu = uicontextmenu;
uimenu(hcmenu,'Label', vvarvecm{i});
......@@ -824,20 +847,28 @@ else
set(h0,'color',a00(i,:),'linewidth',2)
end
ydum=get(gca,'ylim');
if exist('xparam1')
if exist('xparam1','var')
xdum=xparam1(nsnam(j));
h1=plot([xdum xdum],ydum);
set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
end
xlabel('')
title([pnam{nsnam(j)}],'interpreter','none')
if options_.TeX
title([param_names_tex{nsnam(j)}],'interpreter','latex')
else
title([param_names{nsnam(j)}],'interpreter','none')
end
end
%subplot(3,2,6)
if isoctave
legend(vertcat('base',vvarvecm),'location','eastoutside');
else
h0=legend(vertcat('base',vvarvecm));
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none');
if options_.TeX
h0=legend(vertcat('base',vvarvecm_tex),'interpreter','latex');
else
h0=legend(vertcat('base',vvarvecm),'interpreter','none');
end
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3]);
end
if options_.opt_gsa.ppost
dyn_saveas(hh_fig,[ OutDir filesep fname_ '_rmse_post_params_' int2str(ix)],options_.nodisplay,options_.graph_format);
......@@ -885,11 +916,11 @@ pnames=cell(np,1);
pnames_tex=cell(np,1);
for ii=1:length(bayestopt_.name)
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(ii,options_.TeX,M_,estim_params_,options_);
pnames_tex{ii,1} = strrep(param_name_tex_temp,'$','');
[param_name_temp, param_name_tex_temp]= get_the_name(ii,options_.TeX,M_,estim_params_,options_.varobs);
pnames_tex{ii,1} = param_name_tex_temp;
pnames{ii,1} = param_name_temp;
else
param_name_temp = get_the_name(ii,options_.TeX,M_,estim_params_,options_);
param_name_temp = get_the_name(ii,options_.TeX,M_,estim_params_,options_.varobs);
pnames{ii,1} = param_name_temp;
end
end
function indmcf = mcf_analysis(lpmat, ibeha, inobeha, options_mcf, DynareOptions)
function indmcf = monte_carlo_filtering_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
% indmcf = monte_carlo_filtering_analysis(lpmat, ibeha, inobeha, options_mcf, M_, options_, bayestopt_, estim_params_)
% Inputs:
% - lpmat [double] Monte Carlo matrix
% - ibeha [integer] index of behavioural runs
% - inobeha [integer] index of non-behavioural runs
% - options_gsa_ [structure] GSA options_
% - M_ [structure] describing the model
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
% - estim_params_ [structure] characterizing parameters to be estimated
%
% Outputs:
% - indmcf [double] results of matrix
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
%
% Copyright © 2014 European Commission
% Copyright © 2016-2018 Dynare Team
% Copyright © 2016-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -28,7 +41,7 @@ pvalue_corr = options_mcf.pvalue_corr;
alpha2 = options_mcf.alpha2;
param_names = options_mcf.param_names;
if DynareOptions.TeX
if options_.TeX
if ~isfield(options_mcf,'param_names_tex')
param_names_tex = options_mcf.param_names;
else
......@@ -41,6 +54,10 @@ amcf_name = options_mcf.amcf_name;
amcf_title = options_mcf.amcf_title;
beha_title = options_mcf.beha_title;
nobeha_title = options_mcf.nobeha_title;
if options_.TeX
beha_title_latex = options_mcf.beha_title_latex;
nobeha_title_latex = options_mcf.nobeha_title_latex;
end
title = options_mcf.title;
fname_ = options_mcf.fname_;
xparam1=[];
......@@ -49,18 +66,18 @@ if isfield(options_mcf,'xparam1')
end
OutputDirectoryName = options_mcf.OutputDirectoryName;
[proba, dproba] = stab_map_1(lpmat, ibeha, inobeha, [],0);
[proba, dproba] = gsa.stability_mapping_univariate(lpmat, ibeha, inobeha, [],fname_, options_, bayestopt_.name, estim_params_,0);
indmcf=find(proba<pvalue_ks);
[tmp,jtmp] = sort(proba(indmcf),2,'ascend');
[~,jtmp] = sort(proba(indmcf),1,'ascend');
indmcf = indmcf(jtmp);
if ~isempty(indmcf)
skipline()
headers = {'Parameter','d-stat','p-value'};
labels = param_names(indmcf);
data_mat=[dproba(indmcf)' proba(indmcf)'];
data_mat=[dproba(indmcf) proba(indmcf)];
options_temp.noprint=0;
dyntable(options_temp,['Smirnov statistics in driving ', title],headers,labels,data_mat,size(labels,2)+2,16,3);
if DynareOptions.TeX
if options_.TeX
labels_TeX=param_names_tex(indmcf);
M_temp.dname=OutputDirectoryName ;
M_temp.fname=fname_;
......@@ -68,19 +85,31 @@ if ~isempty(indmcf)
end
end
if length(ibeha)>10 && length(inobeha)>10
indcorr1 = stab_map_2(lpmat(ibeha,:),alpha2, pvalue_corr, beha_title);
indcorr2 = stab_map_2(lpmat(inobeha,:),alpha2, pvalue_corr, nobeha_title);
if options_.TeX
indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title, beha_title_latex);
indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title, nobeha_title_latex);
else
indcorr1 = gsa.stability_mapping_bivariate(lpmat(ibeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, beha_title);
indcorr2 = gsa.stability_mapping_bivariate(lpmat(inobeha,:),alpha2, pvalue_corr, M_, options_, bayestopt_, estim_params_, nobeha_title);
end
indcorr = union(indcorr1(:), indcorr2(:));
indcorr = indcorr(~ismember(indcorr(:),indmcf));
indmcf = [indmcf(:); indcorr(:)];
end
if ~isempty(indmcf) && ~DynareOptions.nograph
if ~isempty(indmcf) && ~options_.nograph
skipline()
xx=[];
if ~ isempty(xparam1), xx=xparam1(indmcf); end
scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, DynareOptions, ...
if ~ isempty(xparam1)
xx=xparam1(indmcf);
end
if options_.TeX
gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
beha_title, nobeha_title, beha_title_latex, nobeha_title_latex)
else
gsa.scatter_mcf(lpmat(ibeha,indmcf),lpmat(inobeha,indmcf), param_names_tex(indmcf), ...
'.', [fname_,'_',amcf_name], OutputDirectoryName, amcf_title,xx, options_, ...
beha_title, nobeha_title)
end
end
function [vdec, cc, ac] = mc_moments(mm, ss, dr)
function [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_, estim_params_)
% [vdec, cc, ac] = monte_carlo_moments(mm, ss, dr, M_, options_,estim_params_)
% Conduct Monte Carlo simulation of second moments for GSA
% Inputs:
% - dr [structure] decision rules
% - M_ [structure] model structure
% - options_ [structure] Matlab's structure describing the current options
% - estim_params_ [structure] characterizing parameters to be estimated
%
% Outputs:
% - vdec [double] variance decomposition matrix
% - cc [double] vector of unique elements of cross correlation matrix
% - ac [cell] autocorrelation matrix
% Copyright © 2012-2018 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -17,31 +30,31 @@ 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 <https://www.gnu.org/licenses/>.
global options_ M_ estim_params_ oo_
[nr1, nc1, nsam] = size(mm);
[~, nc1, nsam] = size(mm);
nobs=length(options_.varobs);
disp('Computing theoretical moments ...')
disp('monte_carlo_moments: 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
oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
oo_.dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
dr.ghx = mm(:, 1:(nc1-M_.exo_nbr),j);
dr.ghu = mm(:, (nc1-M_.exo_nbr+1):end, j);
if ~isempty(ss)
set_shocks_param(ss(j,:));
M_=gsa.set_shocks_param(M_,estim_params_,ss(j,:));
end
[vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs);
[vdec(:,:,j), corr, autocorr] = gsa.th_moments(dr,options_,M_);
cc(:,:,j)=triu(corr);
dum=[];
dum=NaN(nobs,nobs*options_.ar);
for i=1:options_.ar
dum=[dum, autocorr{i}];
dum(:,(i-1)*nobs+1:i*nobs)=autocorr{i};
end
ac(:,:,j)=dum;
if mod(j,3)==0
dyn_waitbar(j/nsam,h)
end
end
dyn_waitbar_close(h)
skipline()
disp('... done !')
function pdraw = prior_draw_gsa(init,rdraw)
function pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw)
% pdraw = prior_draw(M_,bayestopt_,options_,estim_params_,init,rdraw)
% Draws from the prior distributions for use with Sensitivity Toolbox for DYNARE
%
% INPUTS
% o init [integer] scalar equal to 1 (first call) or 0.
% o rdraw
% - M_ [structure] describing the model
% - bayestopt_ [structure] describing the priors
% - options_ [structure] describing the options
% - estim_params_ [structure] characterizing parameters to be estimated
% - init [integer] scalar equal to 1 (first call) or 0.
% - rdraw
%
% OUTPUTS
% o pdraw [double] draw from the joint prior density.
......@@ -36,7 +41,6 @@ function pdraw = prior_draw_gsa(init,rdraw)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global bayestopt_ options_ estim_params_ M_
persistent npar pshape p6 p7 p3 p4 lbcum ubcum
if init
......@@ -49,7 +53,7 @@ if init
pdraw = zeros(npar,1);
lbcum = zeros(npar,1);
ubcum = ones(npar,1);
[~,~,~,lb,ub,~] = set_prior(estim_params_,M_,options_); %Prepare bounds
[~,~,~,lb,ub] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
% Set prior bounds
bounds = prior_bounds(bayestopt_, options_.prior_trunc);
......@@ -64,29 +68,29 @@ if init
% set bounds for cumulative probabilities
for i = 1:npar
switch pshape(i)
case 5% Uniform prior.
p4(i) = min(p4(i),bounds.ub(i));
p3(i) = max(p3(i),bounds.lb(i));
case 3% Gaussian prior.
lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 2% Gamma prior.
lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
case 1% Beta distribution (TODO: generalized beta distribution)
lbcum(i) = betainc((bounds.lb(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
ubcum(i) = betainc((bounds.ub(i)-p3(i))./(p4(i)-p3(i)),p6(i),p7(i));
case 2% Gamma prior.
lbcum(i) = gamcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = gamcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
case 3% Gaussian prior.
lbcum(i) = 0.5 * erfc(-(bounds.lb(i)-p6(i))/p7(i) ./ sqrt(2));
ubcum(i) = 0.5 * erfc(-(bounds.ub(i)-p6(i))/p7(i) ./ sqrt(2));
case 4% INV-GAMMA1 distribution
% TO BE CHECKED
lbcum(i) = gamcdf(1/(bounds.ub(i)-p3(i))^2,p7(i)/2,2/p6(i));
ubcum(i) = gamcdf(1/(bounds.lb(i)-p3(i))^2,p7(i)/2,2/p6(i));
case 5% Uniform prior.
p4(i) = min(p4(i),bounds.ub(i));
p3(i) = max(p3(i),bounds.lb(i));
case 6% INV-GAMMA2 distribution
% TO BE CHECKED
lbcum(i) = gamcdf(1/(bounds.ub(i)-p3(i)),p7(i)/2,2/p6(i));
ubcum(i) = gamcdf(1/(bounds.lb(i)-p3(i)),p7(i)/2,2/p6(i));
case 8
lbcum(i) = weibcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = weibcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
lbcum(i) = wblcdf(bounds.lb(i)-p3(i),p6(i),p7(i));
ubcum(i) = wblcdf(bounds.ub(i)-p3(i),p6(i),p7(i));
otherwise
% Nothing to do here.
end
......@@ -94,7 +98,7 @@ if init
return
end
pdraw=NaN(size(rdraw,1),npar);
for i = 1:npar
rdraw(:,i) = rdraw(:,i).*(ubcum(i)-lbcum(i))+lbcum(i);
switch pshape(i)
......
function xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% This procedure transforms x vectors into cumulative values
% pshape: 0 is point mass, both para and p2 are ignored
% 1 is BETA(mean,stdd)
......@@ -11,7 +11,7 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% 8 is WEIBULL(s, k)
% Adapted by M. Ratto from MJ priordens.m
% Copyright © 2012-2015 Dynare Team
% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
......@@ -28,6 +28,7 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
xcum=NaN(size(para));
for i=1:length(pshape)
switch pshape(i)
case 1 % (generalized) BETA Prior
......
function redform_map(dirname,options_gsa_)
%function redform_map(dirname)
% inputs (from opt_gsa structure
% anamendo = options_gsa_.namendo;
% anamlagendo = options_gsa_.namlagendo;
% anamexo = options_gsa_.namexo;
% iload = options_gsa_.load_redform;
% pprior = options_gsa_.pprior;
% ilog = options_gsa_.logtrans_redform;
% threshold = options_gsa_.threshold_redform;
% ksstat = options_gsa_.ksstat_redform;
% alpha2 = options_gsa_.alpha2_redform;
function reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
% reduced_form_mapping(dirname,options_gsa_,M_,estim_params_,options_,bayestopt_,oo_)
% Inputs:
% - dirname [string] name of the output directory
% - options_gsa_ [structure] GSA options_
% - M_ [structure] describing the model
% - estim_params_ [structure] characterizing parameters to be estimated
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
% - oo_ [structure] storing the results
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
......@@ -33,23 +31,16 @@ function redform_map(dirname,options_gsa_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ oo_ estim_params_ options_ bayestopt_
% options_gsa_ = options_.opt_gsa;
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
anamexo = options_gsa_.namexo;
anamendo_tex = options_gsa_.namendo_tex;
anamlagendo_tex = options_gsa_.namlagendo_tex;
anamexo_tex = options_gsa_.namexo_tex;
iload = options_gsa_.load_redform;
pprior = options_gsa_.pprior;
ilog = options_gsa_.logtrans_redform;
threshold = options_gsa_.threshold_redform;
% ksstat = options_gsa_.ksstat_redform;
alpha2 = options_gsa_.alpha2_redform;
alpha2=0;
pvalue_ks = options_gsa_.ksstat_redform;
pvalue_corr = options_gsa_.alpha2_redform;
np = estim_params_.np;
nshock = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
......@@ -57,11 +48,11 @@ pnames=cell(np,1);
pnames_tex=cell(np,1);
for jj=1:np
if options_.TeX
[param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
pnames_tex{jj,1} = strrep(param_name_tex_temp,'$','');
[param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_.varobs);
pnames_tex{jj,1} = param_name_tex_temp;
pnames{jj,1} = param_name_temp;
else
param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_.varobs);
pnames{jj,1} = param_name_temp;
end
end
......@@ -93,14 +84,14 @@ end
options_mcf.fname_ = M_.fname;
options_mcf.OutputDirectoryName = adir;
if ~exist('T')
stab_map_(dirname,options_gsa_);
if ~exist('T','var')
gsa.stability_mapping(dirname,options_gsa_,M_,oo_,options_,bayestopt_,estim_params_);
if pprior
load([dirname,filesep,M_.fname,'_prior'],'T');
else
load([dirname,filesep,M_.fname,'_mc'],'T');
end
if ~exist('T')
if ~exist('T','var')
disp('The model is too large!')
disp('Reduced form mapping stopped!')
return
......@@ -109,8 +100,6 @@ end
if isempty(dir(adir))
mkdir(adir)
end
adir0=pwd;
%cd(adir)
nspred=size(T,2)-M_.exo_nbr;
x0=lpmat(istable,:);
......@@ -121,7 +110,7 @@ else
xx0=lpmat0(istable,:);
nshocks=size(xx0,2);
end
[kn, np]=size(x0);
[~, np]=size(x0);
offset = length(bayestopt_.pshape)-np;
if options_gsa_.prior_range
pshape=5*(ones(np,1));
......@@ -144,23 +133,22 @@ options_map.pshape = pshape;
options_map.pd = pd;
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
lpmat=[];
lpmat0=[];
js=0;
for j = 1:length(anamendo)
namendo = anamendo{j};
namendo_tex = anamendo_tex{j};
iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact');
ifig = 0;
iplo = 0;
for jx = 1:length(anamexo)
namexo = anamexo{jx};
namexo_tex = anamexo_tex{jx};
iexo=strmatch(namexo, M_.exo_names, 'exact');
skipline()
disp(['[', namendo,' vs ',namexo,']'])
if ~isempty(iexo)
%y0=squeeze(T(iendo,iexo+nspred,istable));
y0=squeeze(T(iendo,iexo+nspred,:));
if (max(y0)-min(y0))>1.e-10
if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph
......@@ -194,20 +182,23 @@ for j = 1:length(anamendo)
end
if ~options_.nograph
hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs ', namexo]);
hc = cumplot(y0);
hc = gsa.cumplot(y0);
a=axis; delete(hc);
% hist(mat_moment{ij}),
x1val=max(threshold(1),a(1));
x2val=min(threshold(2),a(2));
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
set(hp,'FaceColor', [0.7 0.8 1])
hold all,
hc = cumplot(y0);
hc = gsa.cumplot(y0);
set(hc,'color','k','linewidth',2)
hold off,
if options_.TeX
title([namendo_tex,' vs ', namexo_tex ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','latex')
else
title([namendo,' vs ', namexo ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','none')
end
dyn_saveas(hf,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],['Reduced Form Mapping (Monte Carlo Filtering): ',strrep(namendo,'_','\_'),' vs ', strrep(namexo,'_','\_')],[type '_' namendo,'_vs_', namexo])
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],['Reduced Form Mapping (Monte Carlo Filtering): ',namendo_tex,' vs ', namexo_tex],[type '_' namendo,'_vs_', namexo])
end
si(:,js) = NaN(np,1);
delete([xdir, '/*threshold*.*'])
......@@ -219,19 +210,23 @@ for j = 1:length(anamendo)
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'inside threshold';
options_mcf.nobeha_title = 'outside threshold';
if options_.TeX
options_mcf.beha_title_latex = 'inside threshold';
options_mcf.nobeha_title_latex = 'outside threshold';
end
options_mcf.title = atitle0;
options_mcf.OutputDirectoryName = xdir;
if ~isempty(iy) && ~isempty(iyc)
fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_);
icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
lpmat=x0(iy,:);
if nshocks
lpmat0=xx0(iy,:);
end
istable=[1:length(iy)];
istable=1:length(iy);
save([xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo '_threshold' ],'lpmat','lpmat0','istable','y0','x0','xx0','iy','iyc')
lpmat=[]; lpmat0=[]; istable=[];
lpmat0=[];
if length(iy)<=10 || length(iyc)<=10
icheck = []; % do the generic plot in any case
end
......@@ -255,12 +250,10 @@ for j = 1:length(anamendo)
end
atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namexo];
options_mcf.title = atitle0;
indmcf = redform_mcf(y0, x0, options_mcf, options_);
redform_mcf(y0, x0, options_mcf, options_, M_.fname, bayestopt_.name, estim_params_);
end
end
else
[yy, xdir] = log_trans_(y0,xdir0);
atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namexo];
aname=[type '_' namendo '_vs_' namexo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namexo];
......@@ -276,24 +269,31 @@ for j = 1:length(anamendo)
figure(hh_fig)
subplot(3,3,iplo),
if ilog
[saso, iso] = sort(-silog(:,js));
[~, iso] = sort(-silog(:,js));
bar([silog(iso(1:min(np,10)),js)])
logflag='log';
else
[saso, iso] = sort(-si(:,js));
[~, iso] = sort(-si(:,js));
bar(si(iso(1:min(np,10)),js))
logflag='';
end
%set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
if options_.TeX
text(ip,-0.02,deblank(pnames_tex(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if options_.TeX
title([logflag,' ',namendo_tex,' vs ',namexo_tex],'interpreter','none')
else
title([logflag,' ',namendo,' vs ',namexo],'interpreter','none')
end
if iplo==9
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namexo,'_','\_')],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],1)
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',namendo_tex,' vs ',namexo_tex],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],1)
end
end
......@@ -310,12 +310,12 @@ for j = 1:length(anamendo)
iplo=0;
for je=1:length(anamlagendo)
namlagendo = anamlagendo{je};
namlagendo_tex = anamlagendo_tex{je};
ilagendo=strmatch(namlagendo, M_.endo_names(oo_.dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
skipline()
disp(['[', namendo,' vs lagged ',namlagendo,']'])
if ~isempty(ilagendo)
%y0=squeeze(T(iendo,ilagendo,istable));
y0=squeeze(T(iendo,ilagendo,:));
if (max(y0)-min(y0))>1.e-10
if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph
......@@ -331,9 +331,9 @@ for j = 1:length(anamendo)
if isempty(dir(xdir0))
mkdir(xdir0)
end
atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs ', namlagendo];
aname=[type '_' namendo '_vs_' namlagendo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo];
atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs lagged', namlagendo];
aname=[type '_' namendo '_vs_lag_' namlagendo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs lagged',namlagendo];
options_map.amap_name = aname;
options_map.amap_title = atitle;
options_map.figtitle = atitle0;
......@@ -349,20 +349,23 @@ for j = 1:length(anamendo)
end
if ~options_.nograph
hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs lagged ', namlagendo]);
hc = cumplot(y0);
hc = gsa.cumplot(y0);
a=axis; delete(hc);
% hist(mat_moment{ij}),
x1val=max(threshold(1),a(1));
x2val=min(threshold(2),a(2));
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
set(hp,'FaceColor', [0.7 0.8 1])
hold all,
hc = cumplot(y0);
hc = gsa.cumplot(y0);
set(hc,'color','k','linewidth',2)
hold off,
hold off
if options_.TeX
title([namendo_tex,' vs lagged ', namlagendo_tex ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','latex')
else
title([namendo,' vs lagged ', namlagendo ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','none')
end
dyn_saveas(hf,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],['Reduced Form Mapping (Monte Carlo Filtering): ',strrep(namendo,'_','\_'),' vs lagged ', strrep(namlagendo,'_','\_')],[type '_' namendo,'_vs_', namlagendo],1)
create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],['Reduced Form Mapping (Monte Carlo Filtering): ',namendo_tex,' vs lagged ', namlagendo_tex],[type '_' namendo,'_vs_', namlagendo],1)
end
delete([xdir, '/*threshold*.*'])
......@@ -374,24 +377,28 @@ for j = 1:length(anamendo)
options_mcf.amcf_title = atitle;
options_mcf.beha_title = 'inside threshold';
options_mcf.nobeha_title = 'outside threshold';
if options_.TeX
options_mcf.beha_title_latex = 'inside threshold';
options_mcf.nobeha_title_latex = 'outside threshold';
end
options_mcf.title = atitle0;
options_mcf.OutputDirectoryName = xdir;
if ~isempty(iy) && ~isempty(iyc)
fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_);
icheck = gsa.monte_carlo_filtering_analysis(x0, iy, iyc, options_mcf, M_, options_, bayestopt_, estim_params_);
lpmat=x0(iy,:);
if nshocks
lpmat0=xx0(iy,:);
end
istable=[1:length(iy)];
istable=1:length(iy);
save([xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo '_threshold' ],'lpmat','lpmat0','istable','y0','x0','xx0','iy','iyc')
lpmat=[]; lpmat0=[]; istable=[];
if length(iy)<=10 || length(iyc)<=10,
lpmat0=[];
if length(iy)<=10 || length(iyc)<=10
icheck = []; % do the generic plot in any case
end
else
icheck = [];
end
......@@ -412,11 +419,10 @@ for j = 1:length(anamendo)
end
atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namlagendo];
options_mcf.title = atitle0;
indmcf = redform_mcf(y0, x0, options_mcf, options_);
redform_mcf(y0, x0, options_mcf, options_, M_.fname, bayestopt_.name, estim_params_);
end
end
else
[yy, xdir] = log_trans_(y0,xdir0);
atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namlagendo];
aname=[type '_' namendo '_vs_' namlagendo];
atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo];
......@@ -432,24 +438,27 @@ for j = 1:length(anamendo)
figure(hh_fig),
subplot(3,3,iplo),
if ilog
[saso, iso] = sort(-silog(:,js));
[~, iso] = sort(-silog(:,js));
bar([silog(iso(1:min(np,10)),js)])
logflag='log';
else
[saso, iso] = sort(-si(:,js));
[~, iso] = sort(-si(:,js));
bar(si(iso(1:min(np,10)),js))
logflag='';
end
%set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
if options_.TeX
text(ip,-0.02,deblank(pnames_tex(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
title([logflag,' ',namendo,' vs ',namlagendo,'(-1)'],'interpreter','none')
if iplo==9
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namlagendo,'_','\_'),'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],1)
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',namendo_tex,' vs ',namlagendo_tex,'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],1)
end
end
......@@ -460,48 +469,37 @@ for j = 1:length(anamendo)
end
if iplo<9 && iplo>0 && ifig && ~options_.nograph
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namlagendo,'_','\_'),'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],options_.figures.textwidth*min(iplo/3,1));
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',namendo_tex,' vs ',namlagendo_tex,'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],options_.figures.textwidth*min(iplo/3,1));
end
end
if isempty(threshold) && ~options_.nograph
hh_fig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA');
if ilog==0
hh_fig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(si)
% boxplot(si','whis',10,'symbol','r.')
myboxplot(si',[],'.',[],10)
gsa.boxplot(si',[],'.',[],10)
else
gsa.boxplot(silog',[],'.',[],10)
end
xlabel(' ')
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np)
set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
if ilog==0
title('Reduced form GSA')
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_gsa'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa'],'Reduced Form GSA','redform_gsa')
else
hh_fig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(silog)
% boxplot(silog','whis',10,'symbol','r.')
myboxplot(silog',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
xlabel(' ')
set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
title('Reduced form GSA - Log-transformed elements')
dyn_saveas(hh_fig,[dirname,filesep,M_.fname,'_redform_gsa_log'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa_log'],'Reduced form GSA - Log-transformed elements','redform_gsa_log')
end
end
function si = redform_private(x0, y0, options_map, options_)
np=size(x0,2);
x00=x0;
ilog = options_map.log_trans;
......@@ -515,7 +513,7 @@ if options_map.prior_range
x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3));
end
else
x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
x0=gsa.priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
end
if ilog
......@@ -531,14 +529,12 @@ if iload==0
nest=max(50,nrun/2);
nest=min(250,nest);
nfit=min(1000,nrun);
% dotheplots = (nfit<=nest);
% gsa_ = gsa_sdp(y0(1:nest), x0(1:nest,:), 2, [],[-1 -1 -1 -1 -1 0],[],0,[fname,'_est'], pnames);
[ys,is] = sort(y0);
[~,is] = sort(y0);
istep = ceil(nrun/nest);
if istep>1
iest = is(floor(istep/2):istep:end);
nest = length(iest);
irest = is(setdiff([1:nrun],[floor(istep/2):istep:nrun]));
irest = is(setdiff(1:nrun,floor(istep/2):istep:nrun));
istep = ceil(length(irest)/(nfit-nest));
ifit = union(iest, irest(1:istep:end));
else
......@@ -550,29 +546,34 @@ if iload==0
ifit = union(ifit, irest(end));
end
nfit=length(ifit);
% ifit = union(iest, irest(randperm(nrun-nest,nfit-nest)));
% ifit = iest;
% nfit=nest;
ipred = setdiff([1:nrun],ifit);
ipred = setdiff(1:nrun,ifit);
if ilog
[y1, tmp, isig, lam] = log_trans_(y0(iest));
[~, ~, isig, lam] = gsa.log_transform(y0(iest));
y1 = log(y0*isig+lam);
end
if ~options_.nograph
hh_fig=dyn_figure(options_.nodisplay,'name',options_map.figtitle);
subplot(221)
if ilog
if isoctave
hist(y1,30)
else
histogram(y1,30)
end
else
if isoctave
hist(y0,30)
else
histogram(y0,30)
end
end
title(options_map.title,'interpreter','none')
subplot(222)
if ilog
hc = cumplot(y1);
hc = gsa.cumplot(y1);
else
hc = cumplot(y0);
hc = gsa.cumplot(y0);
end
set(hc,'color','k','linewidth',2)
title([options_map.title ' CDF'],'interpreter','none')
......@@ -582,15 +583,7 @@ if iload==0
if ilog
[gsa22, gsa1, gsax] = ss_anova_log(y1(iest), x0(iest,:), isig, lam, gsa0);
end
% if (gsa1.out.bic-gsa0.out.bic) < 10,
% y00=y0;
% gsa00=gsa0;
% gsa0=gsa1;
% y0=y1;
% ilog=1;
% end
if nfit>nest
% gsa_ = gsa_sdp(y0(1:nfit), x0(1:nfit,:), -2, gsa_.nvr*nest^3/nfit^3,[-1 -1 -1 -1 -1 0],[],0,fname, pnames);
nvr = gsa0.nvr*nest^3/nfit^3;
nvr(gsa0.stat<2) = gsa0.nvr(gsa0.stat<2)*nest^5/nfit^5;
gsa_ = ss_anova(y0(ifit), x0(ifit,:), 1, 0, 2, nvr);
......@@ -601,33 +594,6 @@ if iload==0
nvrx = gsax.nvr*nest^3/nfit^3;
nvrx(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5;
[gsa22, gsa1, gsax] = ss_anova_log(y1(ifit), x0(ifit,:), isig, lam, gsa0, [nvr1' nvrx']);
% gsa1 = ss_anova(y1(ifit), x0(ifit,:), 1, 0, 2, nvr);
% gsa2=gsa1;
% gsa2.y = gsa0.y;
% gsa2.fit = (exp(gsa1.fit)-lam)*isig;
% gsa2.f0 = mean(gsa2.fit);
% gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
% gsa2.out.bic = gsa2.out.bic-nest*log(gsa1.out.SSE)+nest*log(gsa2.out.SSE);
% gsa2.r2 = 1-cov(gsa2.fit-gsa2.y)/cov(gsa2.y);
% for j=1:np,
% gsa2.fs(:,j) = exp(gsa1.fs(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
% gsa2.f(:,j) = exp(gsa1.f(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
% gsa2.si(j) = var(gsa2.f(:,j))/var(gsa2.y);
% end
% nvr = gsax.nvr*nest^3/nfit^3;
% nvr(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5;
% gsax = ss_anova([gsa2.y-gsa2.fit], x0(ifit,:), 1, 0, 2, nvr);
% gsa22=gsa2;
% gsa22.fit = gsa2.fit+gsax.fit;
% gsa22.f0 = mean(gsa22.fit);
% gsa22.out.SSE = sum((gsa22.fit-gsa22.y).^2);
% gsa22.out.bic = nest*log(gsa22.out.SSE/nest) + (gsax.out.df+gsa2.out.df-1)*log(nest);
% gsa22.r2 = 1-sum((gsa22.fit-gsa22.y).^2)/sum((gsa22.y-mean(gsa22.y)).^2);
% for j=1:np,
% gsa22.fs(:,j) = gsa2.fs(:,j)+gsax.fs(:,j);
% gsa22.f(:,j) = gsa2.f(:,j)+gsax.f(:,j);
% gsa22.si(j) = var(gsa22.f(:,j))/var(gsa22.y);
% end
gsa_ = gsa22;
end
else
......@@ -638,31 +604,23 @@ if iload==0
end
end
save([fname,'_map.mat'],'gsa_')
[sidum, iii]=sort(-gsa_.si);
[~, iii]=sort(-gsa_.si);
gsa_.x0=x00(ifit,:);
if ~options_.nograph
hmap=gsa_sdp_plot(gsa_,[fname '_map'],pnames,iii(1:min(12,np)));
set(hmap,'name',options_map.amap_title);
end
gsa_.x0=x0(ifit,:);
% copyfile([fname,'_est.mat'],[fname,'.mat'])
if ~options_.nograph
figure(hh_fig);
subplot(223),
plot(y0(ifit),[gsa_.fit y0(ifit)],'.'),
r2 = gsa_.r2;
% if ilog,
% plot(y00(ifit),[log_trans_(gsa_.fit,'',isig,lam) y00(ifit)],'.'),
% r2 = 1 - cov(log_trans_(gsa_.fit,'',isig,lam)-y00(ifit))/cov(y00(ifit));
% else
% plot(y0(ifit),[gsa_.fit y0(ifit)],'.'),
% r2 = gsa_.r2;
% end
title(['Learning sample fit - R2=' num2str(r2,2)],'interpreter','none')
if nfit<nrun
if ilog
yf = ss_anova_fcast(x0(ipred,:), gsa1);
yf = log_trans_(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
yf = gsa.log_transform(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
else
yf = ss_anova_fcast(x0(ipred,:), gsa_);
end
......@@ -680,8 +638,6 @@ if iload==0
end
end
else
% gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames);
% gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames);
load([fname,'_map.mat'],'gsa_')
if ~options_.nograph
yf = ss_anova_fcast(x0, gsa_);
......@@ -690,10 +646,8 @@ else
title([namy,' vs ', namx,' pred'],'interpreter','none')
dyn_saveas(hh_fig,[fname '_pred'],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[fname '_pred'],options_map.title,[namy,' vs ', namx,' pred'])
end
end
% si = gsa_.multivariate.si;
si = gsa_.si;
return
......@@ -703,7 +657,7 @@ function gsa2 = log2level_map(gsa1, isig, lam)
nest=length(gsa1.y);
np = size(gsa1.x0,2);
gsa2=gsa1;
gsa2.y = log_trans_(gsa1.y,'',isig,lam);
gsa2.y = gsa.log_transform(gsa1.y,'',isig,lam);
gsa2.fit = (exp(gsa1.fit)-lam)*isig;
gsa2.f0 = mean(gsa2.fit);
gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
......@@ -730,6 +684,8 @@ else
end
gsa2 = log2level_map(gsa1, isig, lam);
if nargin >=5 && ~isempty(gsa0)
nvr2=NaN(np,1);
nvr0=NaN(np,1);
for j=1:np
nvr2(j) = var(diff(gsa2.fs(:,j),2));
nvr0(j) = var(diff(gsa0.fs(:,j),2));
......@@ -760,26 +716,24 @@ end
return
function indmcf = redform_mcf(y0, x0, options_mcf, options_)
function indmcf = redform_mcf(y0, x0, options_mcf, options_, fname, parnames, estim_params_)
hh_fig=dyn_figure(options_.nodisplay,'name',options_mcf.amcf_title);
[post_mean, post_median, post_var, hpd_interval, post_deciles, ...
density] = posterior_moments(y0,1,0.9);
[~, ~, ~, ~, post_deciles] = posterior_moments(y0,0.9);
post_deciles = [-inf; post_deciles; inf];
for jt=1:10
indy{jt}=find( (y0>post_deciles(jt)) & (y0<=post_deciles(jt+1)));
leg{jt}=[int2str(jt) '-dec'];
end
[proba, dproba] = stab_map_1(x0, indy{1}, indy{end}, [],0);
[proba] = gsa.stability_mapping_univariate(x0, indy{1}, indy{end}, [], fname, options_, parnames, estim_params_,0);
indmcf=find(proba<options_mcf.pvalue_ks);
if isempty(indmcf)
[tmp,jtmp] = sort(proba,2,'ascend');
[~,jtmp] = sort(proba,1,'ascend');
indmcf = jtmp(1);
% indmcf = jtmp(1:min(2,length(proba)));
end
[tmp,jtmp] = sort(proba(indmcf),2,'ascend');
[~,jtmp] = sort(proba(indmcf),1,'ascend');
indmcf = indmcf(jtmp);
nbr_par = length(indmcf);
nrow=ceil(sqrt(nbr_par+1));
......@@ -793,12 +747,16 @@ for jx=1:nbr_par
subplot(nrow,ncol,jx)
hold off
for jt=1:10
h=cumplot(x0(indy{jt},indmcf(jx)));
h=gsa.cumplot(x0(indy{jt},indmcf(jx)));
set(h,'color', cmap(jt,:), 'linewidth', 2)
hold all
end
if options_.TeX
title(options_mcf.param_names_tex(indmcf(jx),:),'interpreter','latex')
else
title(options_mcf.param_names(indmcf(jx),:),'interpreter','none')
end
end
hleg = legend(leg);
aa=get(hleg,'Position');
aa(1)=1-aa(3)-0.02;
......@@ -824,7 +782,7 @@ if nargin<5
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([figpath '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by redform_map.m (Dynare).\n');
fprintf(fidTeX,'%% TeX eps-loader file generated by reduced_form_mapping.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
......
function redform_screen(dirname, options_gsa_)
%function redform_map(dirname, options_gsa_)
% inputs (from opt_gsa structure
% anamendo = options_gsa_.namendo;
% anamlagendo = options_gsa_.namlagendo;
% anamexo = options_gsa_.namexo;
% iload = options_gsa_.load_redform;
function reduced_form_screening(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
% reduced_form_screening(dirname, options_gsa_, estim_params_, M_, dr, options_, bayestopt_)
% Conduct reduced form screening
% Inputs:
% - dirname [string] name of the output directory
% - options_gsa_ [structure] GSA options_
% - estim_params [structure] describing the estimated parameters
% - M_ [structure] describing the model
% - dr [structure] decision rules
% - options_ [structure] describing the options
% - bayestopt_ [structure] describing the priors
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
......@@ -28,17 +32,21 @@ function redform_screen(dirname, options_gsa_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ oo_ estim_params_ options_ bayestopt_
% options_gsa_ = options_.opt_gsa;
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
anamexo = options_gsa_.namexo;
iload = options_gsa_.load_redform;
anamendo_tex = options_gsa_.namendo_tex;
anamlagendo_tex = options_gsa_.namlagendo_tex;
anamexo_tex = options_gsa_.namexo_tex;
nliv = options_gsa_.morris_nliv;
pnames = M_.param_names(estim_params_.param_vals(:,1));
if options_.TeX
for par_iter=1:size(estim_params_.param_vals(:,1),1)
[~,tex_names{par_iter,1}]=get_the_name(estim_params_.param_vals(par_iter,1),options_.TeX, M_, estim_params_, options_.varobs);
end
end
if nargin==0
dirname='';
end
......@@ -54,35 +62,45 @@ nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
js=0;
for j=1:size(anamendo,1)
namendo = deblank(anamendo(j,:));
iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact');
namendo = anamendo{j,:};
namendo_tex = anamendo_tex{j,:};
iendo = strmatch(namendo, M_.endo_names(dr.order_var), 'exact');
iplo=0;
ifig=0;
for jx=1:size(anamexo,1)
namexo = deblank(anamexo(jx,:));
namexo = anamexo{jx};
namexo_tex = anamexo_tex{jx};
iexo = strmatch(namexo, M_.exo_names, 'exact');
if ~isempty(iexo)
y0=teff(T(iendo,iexo+nspred,:), kn, istable);
y0=gsa.teff(T(iendo,iexo+nspred,:), kn, istable);
if ~isempty(y0)
if mod(iplo,9)==0
ifig = ifig+1;
hh_fig = dyn_figure(options_.nodisplay, 'name', [namendo,[' vs. shocks '], int2str(ifig)]);
hh_fig = dyn_figure(options_.nodisplay, 'name', [namendo,' vs. shocks ', int2str(ifig)]);
iplo = 0;
end
iplo = iplo+1;
js = js+1;
subplot(3, 3, iplo)
[SAmeas, SAMorris] = Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0, nliv);
[~, SAMorris] = gsa.Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0, nliv);
SAM = squeeze(SAMorris(nshock+1:end,1));
SA(:,js) = SAM./(max(SAM)+eps);
[saso, iso] = sort(-SA(:,js));
[~, iso] = sort(-SA(:,js));
bar(SA(iso(1:min(np,10)),js))
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
if options_.TeX
text(ip,-0.02,tex_names(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,pnames(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if options_.TeX
title([namendo_tex,' vs. ',namexo_tex],'interpreter','latex')
else
title([namendo,' vs. ',namexo],'interpreter','none')
end
if iplo==9
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)],ifig,[namendo,' vs. shocks ',int2str(ifig)],[namendo,'_vs_shock'],1)
......@@ -99,11 +117,12 @@ for j=1:size(anamendo,1)
iplo=0;
ifig=0;
for je=1:size(anamlagendo,1)
namlagendo=deblank(anamlagendo(je,:));
ilagendo=strmatch(namlagendo, M_.endo_names(oo_.dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
namlagendo=anamlagendo{je};
namlagendo_tex=anamlagendo_tex{je};
ilagendo=strmatch(namlagendo, M_.endo_names(dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact');
if ~isempty(ilagendo)
y0=teff(T(iendo,ilagendo,:),kn,istable);
y0=gsa.teff(T(iendo,ilagendo,:),kn,istable);
if ~isempty(y0)
if mod(iplo,9)==0
ifig=ifig+1;
......@@ -113,19 +132,26 @@ for j=1:size(anamendo,1)
iplo=iplo+1;
js=js+1;
subplot(3,3,iplo),
[SAmeas, SAMorris] = Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0,nliv);
[~, SAMorris] = gsa.Morris_Measure_Groups(np+nshock, [lpmat0 lpmat], y0,nliv);
SAM = squeeze(SAMorris(nshock+1:end,1));
SA(:,js)=SAM./(max(SAM)+eps);
[saso, iso] = sort(-SA(:,js));
[~, iso] = sort(-SA(:,js));
bar(SA(iso(1:min(np,10)),js))
%set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
set(gca,'xticklabel',' ','fontsize',10)
set(gca,'xlim',[0.5 10.5])
for ip=1:min(np,10)
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
if options_.TeX
text(ip,-0.02,tex_names(iso(ip)),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,pnames{iso(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
if options_.TeX
title([namendo_tex,' vs. ',namlagendo_tex,'(-1)'],'interpreter','latex')
else
title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
end
if iplo==9
dyn_saveas(hh_fig,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],options_.nodisplay,options_.graph_format);
create_TeX_loader(options_,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)],ifig,[namendo,' vs. lagged endogenous ',int2str(ifig)],[namendo,'_vs_lags'],1)
......@@ -140,15 +166,17 @@ for j=1:size(anamendo,1)
end
hh_fig=dyn_figure(options_.nodisplay,'Name','Reduced form screening');
%bar(SA)
% boxplot(SA','whis',10,'symbol','r.')
myboxplot(SA',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
gsa.boxplot(SA',[],'.',[],10)
set(gca,'xticklabel',' ','fontsize',10,'xtick',1:np)
set(gca,'xlim',[0.5 np+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
for ip=1:np
text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
if options_.TeX
text(ip,-0.02,tex_names(ip),'rotation',90,'HorizontalAlignment','right','interpreter','latex')
else
text(ip,-0.02,pnames{ip},'rotation',90,'HorizontalAlignment','right','interpreter','none')
end
end
xlabel(' ')
ylabel('Elementary Effects')
......@@ -163,7 +191,7 @@ if nargin<6
end
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
fidTeX = fopen([figpath '.tex'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by redform_screen.m (Dynare).\n');
fprintf(fidTeX,'%% TeX eps-loader file generated by reduced_form_screening.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
fprintf(fidTeX,'\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
......
function x0=dynare_sensitivity(options_gsa)
function x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
% x0=run(M_,oo_,options_,bayestopt_,estim_params_,options_gsa)
% Frontend to the Sensitivity Analysis Toolbox for DYNARE
% Inputs:
% - M_ [structure] Matlab's structure describing the model
% - oo_ [structure] Matlab's structure describing the results
% - options_ [structure] Matlab's structure describing the current options
% - bayestopt_ [structure] describing the priors
% - estim_params_ [structure] characterizing parameters to be estimated
% - options_gsa [structure] Matlab's structure describing the GSA options
%
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
% M. Ratto (2008), Analysing DSGE Models with Global Sensitivity Analysis,
% Computational Economics (2008), 31, pp. 115–139
% Copyright © 2008-2018 Dynare Team
% Copyright © 2008-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -21,14 +30,11 @@ function x0=dynare_sensitivity(options_gsa)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
global M_ options_ oo_ bayestopt_ estim_params_
if options_.dsge_var
error('Identification does not support DSGE-VARs at the current stage')
end
fname_ = M_.fname;
lgy_ = M_.endo_names;
x0=[];
% check user defined options
......@@ -43,7 +49,6 @@ end
if isfield(options_gsa,'morris') && options_gsa.morris==1
if isfield(options_gsa,'identification') && options_gsa.identification==0
% options_gsa.redform=1;
end
if isfield(options_gsa,'ppost') && options_gsa.ppost
error('sensitivity:: Morris is incompatible with posterior sampling')
......@@ -89,9 +94,6 @@ if options_.order~=1
options_.order = 1;
end
original_prior_trunc = options_.prior_trunc;
original_qz_criterium = options_.qz_criterium;
if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
if isempty(options_gsa.datafile) && options_gsa.rmse
disp('The data file and all relevant estimation options ')
......@@ -99,6 +101,9 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
disp('must be specified for RMSE analysis!');
error('Sensitivity anaysis error!')
end
if isfield(options_gsa,'nobs')
options_.nobs=options_gsa.nobs;
end
if ~isempty(options_.nobs) && length(options_.nobs)~=1
error('dynare_sensitivity does not support recursive estimation. Please specify nobs as a scalar, not a vector.')
end
......@@ -106,9 +111,6 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
if isfield(options_gsa,'first_obs')
options_.first_obs=options_gsa.first_obs;
end
if isfield(options_gsa,'nobs')
options_.nobs=options_gsa.nobs;
end
if isfield(options_gsa,'presample')
options_.presample=options_gsa.presample;
end
......@@ -130,7 +132,7 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_) || options_gsa.rmse
options_.mode_compute = 0;
options_.filtered_vars = 1;
options_.plot_priors = 0;
[dataset_,dataset_info,xparam1,hh, M_, options_, oo_, estim_params_, bayestopt_] = ...
[dataset_,dataset_info,~,~, 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
......@@ -146,7 +148,12 @@ if M_.exo_nbr==0
error('dynare_sensitivity does not support having no varexo in the model. As a workaround you could define a dummy exogenous variable.')
end
[make,my,day,punk,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
[~,~,~,~,oo_.dr,M_.params] = dynare_resolve(M_,options_,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
if isfield(oo_.dr,'eigval') && any(abs(oo_.dr.eigval-1)<abs(1-options_.qz_criterium)) && options_.qz_criterium<1
fprintf('\ngsa: The model features a unit root, but qz_criterium<1. Check whether that is intended.')
fprintf('\ngsa: If not, use the diffuse_filter option.\n')
end
options_gsa = set_default_option(options_gsa,'identification',0);
if options_gsa.identification
......@@ -166,7 +173,6 @@ if options_gsa.identification
options_.options_ident.load_ident_files = options_gsa.load_ident_files;
options_.options_ident.useautocorr = options_gsa.useautocorr;
options_.options_ident.ar = options_gsa.ar;
options_ident=options_.options_ident;
else
options_ident=[];
options_ident = set_default_option(options_ident,'load_ident_files',options_gsa.load_ident_files);
......@@ -193,7 +199,6 @@ options_gsa = set_default_option(options_gsa,'load_stab',0);
options_gsa = set_default_option(options_gsa,'alpha2_stab',0);
options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
options_gsa = set_default_option(options_gsa,'pvalue_corr',1.e-5);
%options_gsa = set_default_option(options_gsa,'load_mh',0);
% REDFORM mapping
options_gsa = set_default_option(options_gsa,'redform',0);
options_gsa = set_default_option(options_gsa,'load_redform',0);
......@@ -202,8 +207,29 @@ options_gsa = set_default_option(options_gsa,'threshold_redform',[]);
options_gsa = set_default_option(options_gsa,'ksstat_redform',0.001);
options_gsa = set_default_option(options_gsa,'alpha2_redform',1.e-5);
options_gsa = set_default_option(options_gsa,'namendo',{});
options_gsa = set_default_option(options_gsa,'namlagendo',[]);
options_gsa = set_default_option(options_gsa,'namlagendo',{});
options_gsa = set_default_option(options_gsa,'namexo',{});
options_gsa = set_default_option(options_gsa,'namendo_tex',{});
options_gsa = set_default_option(options_gsa,'namlagendo_tex',{});
options_gsa = set_default_option(options_gsa,'namexo_tex',{});
if strmatch(':',options_gsa.namendo,'exact')
options_gsa.namendo = M_.endo_names(1:M_.orig_endo_nbr);
end
if strmatch(':',options_gsa.namexo,'exact')
options_gsa.namexo = M_.exo_names;
end
if strmatch(':',options_gsa.namlagendo,'exact')
options_gsa.namlagendo = M_.endo_names(1:M_.orig_endo_nbr);
end
if options_.TeX
[~,Locb]=ismember(options_gsa.namendo,M_.endo_names);
options_gsa.namendo_tex=cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(Locb), 'UniformOutput', false);
[~,Locb]=ismember(options_gsa.namlagendo,M_.endo_names);
options_gsa.namlagendo_tex=cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(Locb), 'UniformOutput', false);
[~,Locb]=ismember(options_gsa.namexo,M_.exo_names);
options_gsa.namexo_tex=cellfun(@(x) horzcat('$', x, '$'), M_.exo_names_tex(Locb), 'UniformOutput', false);
end
% RMSE mapping
options_gsa = set_default_option(options_gsa,'load_rmse',0);
options_gsa = set_default_option(options_gsa,'lik_only',0);
......@@ -212,8 +238,9 @@ options_gsa = set_default_option(options_gsa,'var_rmse', options_.varobs);
options_gsa.var_rmse_tex={};
for ii=1:length(options_gsa.var_rmse)
temp_name = M_.endo_names_tex{strmatch(options_gsa.var_rmse{ii}, M_.endo_names, 'exact')};
options_gsa.var_rmse_tex = vertcat(options_gsa.var_rmse_tex, temp_name);
options_gsa.var_rmse_tex = vertcat(options_gsa.var_rmse_tex, ['$' temp_name '$']);
end
options_gsa.varobs_tex = cellfun(@(x) horzcat('$', x, '$'), M_.endo_names_tex(options_.varobs_id), 'UniformOutput', false);
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);
......@@ -245,31 +272,22 @@ if options_gsa.morris==1
options_gsa.pprior=1;
end
options_gsa.ppost=0;
%options_gsa.stab=1;
options_gsa.glue=0;
options_gsa.rmse=0;
options_gsa.load_rmse=0;
options_gsa.alpha2_stab=1;
options_gsa.pvalue_ks=0;
options_gsa.pvalue_corr=0;
% if options_gsa.morris==3,
% options_gsa = set_default_option(options_gsa,'Nsam',256);
% OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
% else
OutputDirectoryName = CheckPath('gsa/screen',M_.dname);
% end
else
OutputDirectoryName = CheckPath('gsa',M_.dname);
end
% options_.opt_gsa = options_gsa;
if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform) && options_gsa.pprior
filetoload=[OutputDirectoryName '/' fname_ '_prior.mat'];
if ~exist(filetoload,'file')
disp([filetoload,' not found!'])
disp(['You asked to load a non existent analysis'])
%options_gsa.load_stab=0;
disp('You asked to load a non existent analysis')
return
else
if isempty(strmatch('bkpprior',who('-file', filetoload),'exact'))
......@@ -293,7 +311,7 @@ if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform)
end
if options_gsa.stab && ~options_gsa.ppost
x0 = stab_map_(OutputDirectoryName,options_gsa);
x0 = gsa.stability_mapping(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
if isempty(x0)
skipline()
disp('Sensitivity computations stopped: no parameter set provided a unique solution')
......@@ -301,16 +319,13 @@ if options_gsa.stab && ~options_gsa.ppost
end
end
% reduced form
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
options_.opt_gsa = options_gsa;
if ~isempty(options_gsa.moment_calibration) || ~isempty(options_gsa.irf_calibration)
map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_,bayestopt_);
gsa.map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_,bayestopt_);
end
if options_gsa.identification
map_ident_(OutputDirectoryName,options_gsa);
gsa.map_identification(OutputDirectoryName,options_gsa,M_,oo_,options_,estim_params_,bayestopt_);
end
if options_gsa.redform && ~isempty(options_gsa.namendo)
......@@ -318,7 +333,7 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
filnam = dir([M_.dname filesep 'metropolis' filesep '*param_irf*.mat']);
lpmat=[];
for j=1:length(filnam)
load ([M_.dname filesep 'metropolis' filesep M_.fname '_param_irf' int2str(j) '.mat'])
load ([M_.dname filesep 'metropolis' filesep M_.fname '_param_irf' int2str(j) '.mat'],'stock')
lpmat=[lpmat; stock];
end
clear stock
......@@ -336,35 +351,24 @@ if options_gsa.redform && ~isempty(options_gsa.namendo)
save([OutputDirectoryName filesep M_.fname '_mc.mat'],'lpmat','lpmat0','istable','iunstable','iwrong','iindeterm')
options_gsa.load_stab=1;
x0 = stab_map_(OutputDirectoryName,options_gsa);
end
if strmatch(':',options_gsa.namendo,'exact')
options_gsa.namendo = M_.endo_names(1:M_.orig_endo_nbr);
end
if strmatch(':',options_gsa.namexo,'exact')
options_gsa.namexo = M_.exo_names;
end
if strmatch(':',options_gsa.namlagendo,'exact')
options_gsa.namlagendo = M_.endo_names(1:M_.orig_endo_nbr);
x0 = gsa.stability_mapping(OutputDirectoryName,options_gsa,M_,oo_,options_,bayestopt_,estim_params_);
end
% options_.opt_gsa = options_gsa;
if options_gsa.morris==1
redform_screen(OutputDirectoryName,options_gsa);
gsa.reduced_form_screening(OutputDirectoryName,options_gsa, estim_params_, M_, oo_.dr, options_, bayestopt_);
else
% check existence of the SS_ANOVA toolbox
if isempty(options_gsa.threshold_redform) && ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2)
fprintf('\nThe "SS-ANOVA-R: MATLAB Toolbox for the estimation of Smoothing Spline ANOVA models with Recursive algorithms" is missing.\n')
fprintf('To obtain it, go to:\n\n')
fprintf('https://ec.europa.eu/jrc/en/macro-econometric-statistical-software/ss-anova-r-downloads \n\n')
fprintf('https://joint-research-centre.ec.europa.eu/system/files/2025-01/ss_anova_recurs.zip \n\n')
fprintf('and follow the instructions there.\n')
fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
error('SS-ANOVA-R Toolbox missing!')
end
redform_map(OutputDirectoryName,options_gsa);
gsa.reduced_form_mapping(OutputDirectoryName,options_gsa,M_,estim_params_,options_,bayestopt_,oo_);
end
end
% RMSE mapping
% function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
options_.opt_gsa = options_gsa;
if options_gsa.rmse
if ~options_gsa.ppost
......@@ -391,7 +395,6 @@ if options_gsa.rmse
options_.forecast=0;
options_.filtered_vars=0;
end
% dynare_MC([],OutputDirectoryName,data,rawdata,data_info);
if options_gsa.pprior
TmpDirectoryName = ([M_.dname filesep 'gsa' filesep 'prior']);
else
......@@ -408,37 +411,18 @@ if options_gsa.rmse
delete([TmpDirectoryName filesep filparam(j).name]);
end
end
end
oo_=prior_posterior_statistics('gsa',dataset_, dataset_info,M_,oo_,options_,estim_params_,bayestopt_,'gsa::mcmc');
if options_.bayesian_irf
oo_=PosteriorIRF('gsa',options_,estim_params_,oo_,M_,bayestopt_,dataset_,dataset_info,'gsa::mcmc');
end
options_gsa.load_rmse=0;
% else
% if options_gsa.load_rmse==0,
% disp('You already saved a MC filter/smoother analysis ')
% disp('Do you want to overwrite ?')
% pause;
% if options_gsa.pprior
% delete([OutputDirectoryName,'/',fname_,'_prior_*.mat'])
% else
% delete([OutputDirectoryName,'/',fname_,'_mc_*.mat'])
% end
% dynare_MC([],OutputDirectoryName);
% options_gsa.load_rmse=0;
% end
end
end
clear a;
% filt_mc_(OutputDirectoryName,data_info);
filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info);
gsa.monte_carlo_filtering(OutputDirectoryName,options_gsa,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_);
end
options_.opt_gsa = options_gsa;
options_.prior_trunc=original_prior_trunc;
options_.qz_criterium=original_qz_criterium ;
if options_gsa.glue
dr_ = oo_.dr;
......@@ -453,11 +437,10 @@ if options_gsa.glue
end
end
if ~exist('x','var')
disp(['No RMSE analysis is available for current options'])
disp(['No GLUE file prepared'])
disp('No RMSE analysis is available for current options')
disp('No GLUE file prepared')
return,
end
nruns=size(x,1);
gend = options_.nobs;
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
......@@ -465,28 +448,20 @@ if options_gsa.glue
rawdata = log(rawdata);
end
if options_.prefilter == 1
%data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
data = transpose(rawdata-ones(gend,1)*mean(rawdata,1));
else
data = transpose(rawdata);
end
Obs.data = data;
Obs.time = [1:gend];
Obs.time = 1:gend;
Obs.num = gend;
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');
jxj = strmatch(vj,M_.endo_names(dr_.order_var),'exact');
if ~options_gsa.ppost
% y0=zeros(gend+1,nruns);
% nb = size(stock_filter,3);
% y0 = squeeze(stock_filter(:,jxj,:)) + ...
% kron(stock_ys(js,:),ones(size(stock_filter,1),1));
% Out(j).data = y0';
% Out(j).time = [1:size(y0,1)];
Out(j).data = jxj;
Out(j).time = [pwd,'/',OutputDirectoryName];
else
......@@ -500,17 +475,7 @@ if options_gsa.glue
Lik(j).isam = 1;
Lik(j).data = rmse_MC(:,j)';
if ~options_gsa.ppost
% y0 = squeeze(stock_smooth(:,jxj,:)) + ...
% kron(stock_ys(js,:),ones(size(stock_smooth,1),1));
% Out1(j).name = vj;
% Out1(j).ini = 'yes';
% Out1(j).time = [1:size(y0,1)];
% Out1(j).data = y0';
Out1=Out;
else
Out1=Out;
end
ismoo(j)=jxj;
end
......@@ -520,10 +485,6 @@ if options_gsa.glue
jsmoo=jsmoo+1;
vj = M_.endo_names{dr_.order_var(j)};
if ~options_gsa.ppost
% y0 = squeeze(stock_smooth(:,j,:)) + ...
% kron(stock_ys(j,:),ones(size(stock_smooth,1),1));
% Out1(jsmoo).time = [1:size(y0,1)];
% Out1(jsmoo).data = y0';
Out1(jsmoo).data = j;
Out1(jsmoo).time = [pwd,'/',OutputDirectoryName];
else
......@@ -546,36 +507,24 @@ if options_gsa.glue
end
Sam.name = bayestopt_.name;
Sam.dim = [size(x) 0];
Sam.data = [x];
Sam.data = x;
Rem.id = 'Original';
Rem.ind= [1:size(x,1)];
Rem.ind= 1:size(x,1);
Info.dynare=M_.fname;
Info.order_var=dr_.order_var;
Out=Out1;
if options_gsa.ppost
% Info.dynare=M_.fname;
% Info.order_var=dr_.order_var;
% Out=Out1;
Info.TypeofSample='post';
save([OutputDirectoryName,'/',fname_,'_glue_post.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
%save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info')
else
if options_gsa.pprior
Info.TypeofSample='prior';
save([OutputDirectoryName,'/',fname_,'_glue_prior.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
% save([OutputDirectoryName,'/',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
% Out=Out1;
% save([OutputDirectoryName,'/',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
else
Info.TypeofSample='mc';
save([OutputDirectoryName,'/',fname_,'_glue_mc.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
% save([OutputDirectoryName,'/',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
% Out=Out1;
% save([OutputDirectoryName,'/',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
end
end
end