diff --git a/license.txt b/license.txt
index 71083f2ab..aa420ab47 100644
--- a/license.txt
+++ b/license.txt
@@ -171,16 +171,13 @@ Comment: Written by Jessica Cariboni and Francesca Campolongo
Files: matlab/gsa/cumplot.m
matlab/gsa/filt_mc_.m
- matlab/gsa/gsa_plotmatrix.m
matlab/gsa/gsa_skewness.m
matlab/gsa/log_trans_.m
matlab/gsa/map_calibration.m
matlab/gsa/map_ident_.m
matlab/gsa/mcf_analysis.m
matlab/gsa/myboxplot.m
- matlab/gsa/myprctilecol.m
matlab/gsa/prior_draw_gsa.m
- matlab/gsa/read_data.m
matlab/gsa/redform_map.m
matlab/gsa/redform_screen.m
matlab/gsa/scatter_mcf.m
@@ -191,9 +188,8 @@ Files: matlab/gsa/cumplot.m
matlab/gsa/stand_.m
matlab/gsa/tcrit.m
matlab/gsa/teff.m
- matlab/gsa/trank.m
Copyright: 2011-2018 European Commission
- 2011-2018 Dynare Team
+ 2011-2023 Dynare Team
License: GPL-3+
Files: matlab/gsa/pick.m
diff --git a/matlab/gsa/cumplot.m b/matlab/gsa/cumplot.m
index 9de34e86c..c58567a4b 100644
--- a/matlab/gsa/cumplot.m
+++ b/matlab/gsa/cumplot.m
@@ -1,5 +1,10 @@
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
diff --git a/matlab/gsa/gsa_plotmatrix.m b/matlab/gsa/gsa_plotmatrix.m
deleted file mode 100644
index 0b5caa6bf..000000000
--- a/matlab/gsa/gsa_plotmatrix.m
+++ /dev/null
@@ -1,99 +0,0 @@
-function gsa_plotmatrix(type,varargin)
-% function gsa_plotmatrix(type,varargin)
-% extended version of the standard MATLAB plotmatrix
-%
-% Written by Marco Ratto
-% Joint Research Centre, The European Commission,
-% marco.ratto@ec.europa.eu
-
-% Copyright © 2011-2012 European Commission
-% Copyright © 2011-2017 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 .
-
-global bayestopt_ options_ M_
-
-RootDirectoryName = CheckPath('gsa',M_.dname);
-
-if options_.opt_gsa.pprior
- load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable','iunstable','iindeterm','iwrong')
-else
- load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable','iunstable','iindeterm','iwrong')
- eval(['load ' options_.mode_file ' xparam1;']');
-end
-
-iexplosive = iunstable(~ismember(iunstable,[iindeterm;iwrong]));
-
-switch type
- case 'all'
- x=[lpmat0 lpmat];
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
- case 'stable'
- x=[lpmat0(istable,:) lpmat(istable,:)];
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
- case 'nosolution'
- x=[lpmat0(iunstable,:) lpmat(iunstable,:)];
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
- case 'unstable'
- x=[lpmat0(iexplosive,:) lpmat(iexplosive,:)];
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
- case 'indeterm'
- x=[lpmat0(iindeterm,:) lpmat(iindeterm,:)];
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
- case 'wrong'
- x=[lpmat0(iwrong,:) lpmat(iwrong,:)];
- NumberOfDraws=size(x,1);
- B=NumberOfDraws;
-
-end
-
-if isempty(x)
- disp('Empty parameter set!')
- return
-end
-
-for j=1:length(varargin)
- jcol(j)=strmatch(varargin{j},bayestopt_.name,'exact');
-end
-
-[H,AX,BigA,P,PAx]=plotmatrix(x(:,jcol));
-
-for j=1:length(varargin)
- % axes(AX(1,j)), title(varargin{j})
- % axes(AX(j,1)), ylabel(varargin{j})
- % set(AX(1,j),'title',varargin{j}),
- set(get(AX(j,1),'ylabel'),'string',varargin{j})
- set(get(AX(end,j),'xlabel'),'string',varargin{j})
-end
-
-if options_.opt_gsa.pprior==0
- xparam1=xparam1(jcol);
- for j=1:length(varargin)
- for i=1:j-1
- axes(AX(j,i))
- hold on, plot(xparam1(i),xparam1(j),'*r')
- end
- for i=j+1:length(varargin)
- axes(AX(j,i))
- hold on, plot(xparam1(i),xparam1(j),'*r')
- end
- end
-end
diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m
index ec84adef2..2b7194fa2 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/gsa/map_ident_.m
@@ -365,3 +365,13 @@ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
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
diff --git a/matlab/gsa/myboxplot.m b/matlab/gsa/myboxplot.m
index 7864107e2..4d6cf60d1 100644
--- a/matlab/gsa/myboxplot.m
+++ b/matlab/gsa/myboxplot.m
@@ -1,12 +1,8 @@
function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
-% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
+% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
+% Creates a box plot
-% Written by Marco Ratto
-% Joint Research Centre, The European Commission,
-% marco.ratto@ec.europa.eu
-
-% Copyright © 2012 European Commission
-% Copyright © 2012-2017 Dynare Team
+% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -23,18 +19,17 @@ 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 .
-% % % % 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
-if nargin < 2 | isempty(notched), notched = 0; end
+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
+if nargin < 2 || isempty(notched), notched = 0; end
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
@@ -42,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];
@@ -57,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))];
@@ -102,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;
@@ -146,11 +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 vertical
plot (quartile_x, quartile_y, 'b', ...
@@ -162,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
diff --git a/matlab/gsa/myprctilecol.m b/matlab/gsa/myprctilecol.m
deleted file mode 100644
index 1baafb1e0..000000000
--- a/matlab/gsa/myprctilecol.m
+++ /dev/null
@@ -1,43 +0,0 @@
-function y = myprctilecol(x,p)
-
-% Written by Marco Ratto
-% Joint Research Centre, The European Commission,
-% marco.ratto@ec.europa.eu
-
-% Copyright © 2012 European Commission
-% Copyright © 2012-2017 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 .
-
-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
diff --git a/matlab/gsa/priorcdf.m b/matlab/gsa/priorcdf.m
index 1be3e1f14..ac0ffc1a6 100644
--- a/matlab/gsa/priorcdf.m
+++ b/matlab/gsa/priorcdf.m
@@ -1,5 +1,5 @@
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 .
+xcum=NaN(size(para));
for i=1:length(pshape)
switch pshape(i)
case 1 % (generalized) BETA Prior
diff --git a/matlab/gsa/read_data.m b/matlab/gsa/read_data.m
deleted file mode 100644
index 5af0f5204..000000000
--- a/matlab/gsa/read_data.m
+++ /dev/null
@@ -1,44 +0,0 @@
-function [gend, data] = read_data()
-% Written by Marco Ratto
-% Joint Research Centre, The European Commission,
-% marco.ratto@ec.europa.eu
-
-% Copyright © 2012-2015 European Commission
-% Copyright © 2012-2017 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 .
-
-global options_
-
-rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
-
-options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
-gend = options_.nobs;
-
-rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
-if options_.loglinear == 1 & ~options_.logdata
- rawdata = log(rawdata);
-end
-if options_.prefilter == 1
- data = transpose(rawdata-ones(gend,1)* mean(rawdata,1));
-else
- data = transpose(rawdata);
-end
-
-if ~isreal(rawdata)
- error(['There are complex values in the data. Probably a wrong' ...
- ' transformation'])
-end
diff --git a/matlab/gsa/stand_.m b/matlab/gsa/stand_.m
index 176af7f5c..c95ccf143 100644
--- a/matlab/gsa/stand_.m
+++ b/matlab/gsa/stand_.m
@@ -1,20 +1,22 @@
function [y, meany, stdy] = stand_(x)
-% STAND_ Standardise a matrix by columns
+% [y, meany, stdy] = stand_(x)
+% Standardise a matrix by columns
%
% [x,my,sy]=stand(y)
%
-% y: Time series (column matrix)
+% Inputs:
+% - x: Time series (column matrix)
%
-% x: standardised equivalent of y
-% my: Vector of mean values for each column of y
-% sy: Vector of standard deviations for each column of y
+% - y: standardised equivalent of x
+% - meany: Vector of mean values for each column of x
+% - stdy: Vector of standard deviations for each column of x
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
% marco.ratto@ec.europa.eu
% Copyright © 2012 European Commission
-% Copyright © 2012-2017 Dynare Team%
+% Copyright © 2012-2023 Dynare Team
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
@@ -34,9 +36,11 @@ if nargin==0
return
end
+meany=NaN(size(x,2),1);
+stdy=NaN(size(x,2),1);
+y=NaN(size(x));
for j=1:size(x,2)
- meany(j)=mean(x(find(~isnan(x(:,j))),j));
- stdy(j)=std(x(find(~isnan(x(:,j))),j));
+ meany(j)=mean(x(~isnan(x(:,j)),j));
+ stdy(j)=std(x(~isnan(x(:,j)),j));
y(:,j)=(x(:,j)-meany(j))./stdy(j);
-end
-% end of m-file
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/matlab/gsa/tcrit.m b/matlab/gsa/tcrit.m
index c6b17f4bc..120a6107b 100644
--- a/matlab/gsa/tcrit.m
+++ b/matlab/gsa/tcrit.m
@@ -150,4 +150,4 @@ if n<=100
t_crit=t_crit(n,ncol);
else
t_crit=t_crit(end,ncol);
-end
\ No newline at end of file
+end
diff --git a/matlab/gsa/teff.m b/matlab/gsa/teff.m
index 50c4917ee..b815237bd 100644
--- a/matlab/gsa/teff.m
+++ b/matlab/gsa/teff.m
@@ -37,15 +37,12 @@ if ndim==3
[ir, ic]=(find( (tmax-tmin)>1.e-8));
j0 = length(ir);
yt=zeros(Nsam, j0);
-
for j=1:j0
y0=squeeze(T(ir(j),ic(j),:));
- %y1=ones(size(lpmat,1),1)*NaN;
y1=ones(Nsam,1)*NaN;
y1(istable,1)=y0;
yt(:,j)=y1;
end
-
else
tmax=max(T,[],2);
tmin=min(T,[],2);
@@ -53,7 +50,4 @@ else
j0 = length(ir);
yt=NaN(Nsam, j0);
yt(istable,:)=T(ir,:)';
-
-
end
-%clear y0 y1;
diff --git a/matlab/gsa/trank.m b/matlab/gsa/trank.m
deleted file mode 100644
index 41a200a94..000000000
--- a/matlab/gsa/trank.m
+++ /dev/null
@@ -1,34 +0,0 @@
-function yr = trank(y)
-% yr = trank(y);
-% yr is the rank transformation of y
-%
-% Written by Marco Ratto
-% Joint Research Centre, The European Commission,
-% marco.ratto@ec.europa.eu
-%
-% Reference:
-% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
-
-% Copyright © 2012 European Commission
-% Copyright © 2012-2017 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 .
-
-[nr, nc] = size(y);
-for j=1:nc
- [dum, is]=sort(y(:,j));
- yr(is,j)=[1:nr]'./nr;
-end
diff --git a/matlab/myboxplot.m b/matlab/myboxplot.m
deleted file mode 100644
index 73eb78987..000000000
--- a/matlab/myboxplot.m
+++ /dev/null
@@ -1,177 +0,0 @@
-
-function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
-
-% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
-
-%
-% Copyright © 2010-2017 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 .
-
-% % % % 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
-if nargin < 2 || isempty(notched), notched = 0; end
-
-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
-if iscell(data)
- nc = length(data);
-else
- % if isvector(data), data = data(:); end
- 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
-s = zeros(7,nc);
-box = zeros(1,nc);
-whisker_x = ones(2,1)*[1:nc,1:nc];
-whisker_y = zeros(2,2*nc);
-outliers_x = [];
-outliers_y = [];
-outliers2_x = [];
-outliers2_y = [];
-
-for i=1:nc
- % ## 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
- % % % % % % % col(isnan(col) | isna (col)) = [];
- col(isnan(col)) = [];
-
- % ## 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);
- 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
- 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
- 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 = 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))];
- outliers_y = [outliers_y; outliers];
- 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
- s(:,i) = col;
- % ## single point data sets are plotted as outliers.
- outliers_x = [outliers_x; i];
- outliers_y = [outliers_y; col];
- else
- % ## no statistics if no points
- s(:,i) = NaN;
- end
-end
-% % % % if isempty(outliers2_y)
-% % % % outliers2_y=
-% ## 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.
-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
-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
-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
-cap_x = whisker_x;
-cap_x(1,:) =cap_x(1,:)- 0.05;
-cap_x(2,:) =cap_x(2,:)+ 0.05;
-cap_y = whisker_y([1,1],:);
-
-% #quartile_x,quartile_y
-% #whisker_x,whisker_y
-% #median_x,median_y
-% #cap_x,cap_y
-%
-% ## Do the plot
-
-mm=min(min(data));
-MM=max(max(data));
-if isnan(mm), mm=0; MM=0; end
-
-if vertical
- plot (quartile_x, quartile_y, 'b', ...
- whisker_x, whisker_y, 'b--', ...
- cap_x, cap_y, 'k', ...
- median_x, median_y, 'r', ...
- outliers_x, outliers_y, [symbol(1),'r'], ...
- outliers2_x, outliers2_y, [symbol(2),'r']);
- 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