Skip to content
Snippets Groups Projects
Commit 753aa48a authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

GSA: remove tags for older Dynare versions

parent fe4dbbfe
Branches
Tags
No related merge requests found
Showing
with 0 additions and 3414 deletions
File deleted
This diff is collapsed.
function [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)
% [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)
%
% Given the Morris sample matrix, the output values and the group matrix compute the Morris measures
% -------------------------------------------------------------------------
% INPUTS
% -------------------------------------------------------------------------
% Group [NumFactor, NumGroups] := Matrix describing the groups.
% Each column represents one group.
% The element of each column are zero if the factor is not in the
% group. Otherwise it is 1.
% Sample := Matrix of the Morris sampled trajectories
% Output := Matrix of the output(s) values in correspondence of each point
% of each trajectory
% k = Number of factors
% -------------------------------------------------------------------------
% OUTPUTS
% OutMatrix (NumFactor*NumOutputs, 3)= [Mu*, Mu, StDev]
% for each output it gives the three measures of each factor
% -------------------------------------------------------------------------
if nargin==0,
disp(' ')
disp('[SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group);')
return
end
OutMatrix=[];
if nargin < 5, Group=[]; end
NumGroups = size(Group,2);
if nargin < 4 | isempty(p),
p = 4;
end
Delt = p/(2*p-2);
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
% 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
Single_Sample = Sample(i+(i-1)*sizea:i+(i-1)*sizea+sizea,:);
Single_OutValues = OutValues(i+(i-1)*sizea:i+(i-1)*sizea+sizea,:);
A = (Single_Sample(2:sizea+1,:)-Single_Sample(1:sizea,:))';
Delta = A(find(A));
% 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.
for j=1:sizea % For each point in the trajectory i.e for each factor
% matrix of factor which changes
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
% opposite direction and the definition of the standard
% 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,:));
% If no groups --> we compute both the original and
% modified measure
if Delta(j) > 0 %=> +Delta
SAmeas(Change_factor(j,i),i) = (Single_OutValues(j+1) - Single_OutValues(j) )/Delt; %(2/3);
else %=> -Delta
SAmeas(Change_factor(j,i),i) = (Single_OutValues(j) - Single_OutValues(j+1) )/Delt; %(2/3);
end
end
end %for j=1:sizea
end %for i=1:r
if NumGroups ~ 0
SAmeas = SAmeas';
end
% Compute Mu AbsMu and StDev
if any(any(isnan(SAmeas)))
for j=1:NumFact,
SAm = SAmeas(j,:);
SAm = SAm(find(~isnan(SAm)));
rr=length(SAm);
AbsMu(j,1) = sum(abs(SAm),2)/rr;
if NumGroups == 0
Mu(j,1) = sum(SAm,2)/rr;
StDev(j,1) = sum((SAm - repmat(Mu(j),1,rr)).^2/(rr*(rr-1)),2).^0.5;
end
end
else
AbsMu = sum(abs(SAmeas),2)/r;
if NumGroups == 0
Mu = sum(SAmeas,2)/r;
StDev = sum((SAmeas - repmat(Mu,1,r)).^2/(r*(r-1)),2).^0.5;
end
end
% 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];
else
OutMatrix = [OutMatrix; 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.
% 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
% r (1,1) := sample size
% p (1,1) := number of intervals considered in [0, 1]
% UB(sizea,1) := Upper Bound for each factor
% LB(sizea,1) := Lower Bound for each factor
% GroupNumber(1,1) := Number of groups (eventually 0)
% GroupMat(NumFact,GroupNumber):= Matrix which describes the chosen groups. Each column represents a group and its elements
% are set to 1 in correspondence of the factors that belong to the fixed group. All
% the other elements are zero.
% Local Variables:
% sizeb (1,1) := sizea+1
% sizec (1,1) := 1
% randmult (sizea,1) := vector of random +1 and -1
% perm_e(1,sizea) := vector of sizea random permutated indeces
% fact(sizea) := vector containing the factor varied within each traj
% DDo(sizea,sizea) := D* in Morris, 1991
% A(sizeb,sizea) := Jk+1,k in Morris, 1991
% B(sizeb,sizea) := B in Morris, 1991
% Po(sizea,sizea) := P* in Morris, 1991
% Bo(sizeb,sizea) := B* in Morris, 1991
% Ao(sizeb,sizec) := Jk+1,1 in Morris, 1991
% xo(sizec,sizea) := x* in Morris, 1991 (starting point for the trajectory)
% In(sizeb,sizea) := for each loop orientation matrix. It corresponds to a trajectory
% of k step in the parameter space and it provides a single elementary
% effect per factor
% MyInt()
% Fact(sizea,1) := for each loop vector indicating which factor or group of factors has been changed
% in each step of the trajectory
% 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
% 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
% follows the following steps:
% 1- Creation of P0 and DD0 matrices defined in Morris for the groups. This means that the dimensions of these
% 2 matrices are (GroupNumber,GroupNumber).
% 2- Creation of AuxMat matrix with (GroupNumber+1,GroupNumber) elements.
% 3- Definition of GroupB0 starting from AuxMat, GroupMat and P0.
% 4- The final B0 for groups is obtained as [ones(sizeb,1)*x0' + GroupB0]. The P0 permutation is present in GroupB0
% and it's not necessary to permute the matrix (ones(sizeb,1)*x0') because it's already randomly created.
% Reference:
% A. Saltelli, K. Chan, E.M. Scott, "Sensitivity Analysis" on page 68 ss
%
% F. Campolongo, J. Cariboni, JRC - IPSC Ispra, Varese, IT
% Last Update: 15 November 2005 by J.Cariboni
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters and initialisation of the output matrix
sizea = k;
Delta = p/(2*p-2);
%Delta = 1/3
NumFact = sizea;
GroupNumber = size(GroupMat,2);
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 = repmat(randmult,1,k);
DD0 = randmult .* eye(k);
% Construct DD0 - NEW VERSION - it needs communication toolbox
% randsrc(m) generates an m-by-m matrix, each of whose entries independently takes the value -1 with probability 1/2,
% and 1 with probability 1/2.
% DD0 = randsrc(NumFact) .* eye(NumFact);
% Construct B (lower triangular)
B = ones(sizeb,sizea);
for j = 1:sizea
B(1:j,j)=0;
end
% Construct A0, A
A0 = ones(sizeb,1);
A = ones(sizeb,NumFact);
% Construct the permutation matrix P0. In each column of P0 one randomly chosen element equals 1
% while all the others equal zero.
% P0 tells the order in which order factors are changed in each
% trajectory. P0 is created as it follows:
% 1) All the elements of P0 are set equal to zero ==> P0 = zeros (sizea, sizea);
% 2) The function randperm create a random permutation of integer 1:sizea, without repetitions ==> perm_e;
% 3) In each column of P0 the element indicated in perm_e is set equal to one.
% Note that P0 is then used reading it by rows.
P0 = zeros (sizea, sizea);
perm_e = randperm(sizea); % RANDPERM(n) is a random permutation of the integers from 1 to n.
for j = 1:sizea
P0(perm_e(j),j) = 1;
end
% 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
B = B * (GroupMat*P0')';
end
% Compute AuxMat both for single factors and groups analysis. For Single factors analysis
% AuxMat is added to (A0*X0) and then permutated through P0. When groups are active AuxMat is
% used to build GroupB0. AuxMat is created considering DD0. If the element on DD0 diagonal
% is 1 then AuxMat will start with zero and add Delta. If the element on DD0 diagonal is -1
% then DD0 will start Delta and goes to zero.
AuxMat = Delta*0.5*((2*B - A) * DD0 + A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
% OLD VERSION - it needs communication toolbox
% w = randint(NumFact,1,[1,size(MyInt,2)]);
% NEW VERSION - construct a version of random integers
% 1) create a vector of random integers
% 2) divide [0,1] into the needed steps
% 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)
DiffAuxVec = IntUsed - v; % 3)
for ii = 1:size(DiffAuxVec,1)
w(1,ii) = max(find(DiffAuxVec(ii,:)<0)); % 4)
end
x0 = MyInt(1,w)'; % Define x0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
B0 = (A0*x0' + AuxMat);
else
B0 = (A0*x0' + AuxMat)*P0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% c --> Compute values in the original intervals
% B0 has values x(i,j) in [0, 1/(p -1), 2/(p -1), ... , 1].
% To obtain values in the original intervals [LB, UB] we compute
% LB(j) + x(i,j)*(UB(j)-LB(j))
In = repmat(LB,1,sizeb)' + B0 .* repmat((UB-LB),1,sizeb)';
% 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.
for j=1:sizea
Fact(1,j) = find(P0(j,:));
end
Fact(1,sizea+1) = 0;
Outmatrix = [Outmatrix; In];
OutFact = [OutFact; Fact'];
end
\ No newline at end of file
function x = beta_inv(p, a, b)
% PURPOSE: inverse of the cdf (quantile) of the beta(a,b) distribution
%--------------------------------------------------------------
% USAGE: x = beta_inv(p,a,b)
% where: p = vector of probabilities
% a = beta distribution parameter, a = scalar
% b = beta distribution parameter b = scalar
% NOTE: mean [beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
%--------------------------------------------------------------
% RETURNS: x at each element of p for the beta(a,b) distribution
%--------------------------------------------------------------
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
%--------------------------------------------------------------
% Anders Holtsberg, 18-11-93
% Copyright (c) Anders Holtsberg
% documentation modified by LeSage to
% match the format of the econometrics toolbox
if (nargin ~= 3)
error('Wrong # of arguments to beta_inv');
end
if any(any((a<=0)|(b<=0)))
error('beta_inv parameter a or b is nonpositive');
end
if any(any(abs(2*p-1)>1))
error('beta_inv: A probability should be 0<=p<=1');
end
x = a ./ (a+b);
dx = 1;
while any(any(abs(dx)>256*eps*max(x,1)))
dx = (betainc(x,a,b) - p) ./ beta_pdf(x,a,b);
x = x - dx;
x = x + (dx - x) / 2 .* (x<0);
end
\ No newline at end of file
function pdf = beta_pdf(x, a, b)
% PURPOSE: pdf of the beta(a,b) distribution
%--------------------------------------------------------------
% USAGE: pdf = beta_pdf(x,a,b)
% where: x = vector of components
% a = beta distribution parameter, a = scalar
% b = beta distribution parameter b = scalar
% NOTE: mean[(beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
%--------------------------------------------------------------
% RETURNS: pdf at each element of x of the beta(a,b) distribution
%--------------------------------------------------------------
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
%--------------------------------------------------------------
% Anders Holtsberg, 18-11-93
% Copyright (c) Anders Holtsberg
% documentation modified by LeSage to
% match the format of the econometrics toolbox
if (nargin ~=3)
error('Wrong # of arguments to beta_pdf');
end
if any(any((a<=0)|(b<=0)))
error('Parameter a or b is nonpositive');
end
I = find((x<0)|(x>1));
pdf = x.^(a-1) .* (1-x).^(b-1) ./ beta(a,b);
pdf(I) = 0*I;
function h = cumplot(x);
%function h =cumplot(x)
% Copyright (C) 2005 Marco Ratto
n=length(x);
x=[-inf; sort(x); Inf];
y=[0:n n]./n;
h0 = stairs(x,y);
grid on,
if nargout,
h=h0;
end
function c = dat_fil_(data_file);
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
try
eval(data_file);
catch
load(data_file);
end
clear data_file;
a=who;
for j=1:length(a)
eval(['c.',a{j},'=',a{j},';']);
end
\ No newline at end of file
function dynare_MC(var_list_,OutDir)
%
% Adapted by M. Ratto from dynare_estimation.m and posteriorsmoother.m
% (dynare_estimation.m and posteriorsmoother.m are part of DYNARE,
% copyright M. Juillard)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_ options_ oo_ estim_params_
global bayestopt_
if options_.filtered_vars ~= 0 & options_.filter_step_ahead == 0
options_.filter_step_ahead = 1;
end
if options_.filter_step_ahead ~= 0
options_.nk = max(options_.filter_step_ahead);
else
options_.nk = 0;
end
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
if isempty(options_.datafile)
error('ESTIMATION: datafile option is missing')
end
if isempty(options_.varobs)
error('ESTIMATION: VAROBS is missing')
end
%% Read and demean data
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
bayestopt_.mean_varobs = mean(rawdata,1);
data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
else
data = transpose(rawdata);
end
if ~isreal(rawdata)
error(['There are complex values in the data. Probably a wrong' ...
' transformation'])
end
offset = npar-np;
fname_=M_.fname;
options_ = set_default_option(options_,'opt_gsa',1);
options_gsa_ = options_.opt_gsa;
if options_gsa_.pprior,
namfile=[fname_,'_prior'];
else
namfile=[fname_,'_mc'];
end
load([OutDir,'\',namfile],'lpmat', 'lpmat0', 'istable')
% load(options_.mode_file)
%%
%%
%%
x=[lpmat0(istable,:) lpmat(istable,:)];
clear lpmat lpmat0 istable %iunstable egg yys T
B = size(x,1);
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(x(1,:)',gend,data);
n1=size(atT,1);
nfil=B/40;
stock_smooth = zeros(M_.endo_nbr,gend,40);
stock_filter = zeros(M_.endo_nbr,gend+1,40);
stock_ys = zeros(40, M_.endo_nbr);
logpo2=zeros(B,1);
%%
h = waitbar(0,'MC smoother ...');
delete([OutDir,'\',namfile,'_*.mat'])
ib=0;
ifil=0;
opt_gsa=options_.opt_gsa;
for b=1:B
ib=ib+1;
deep = x(b,:)';
set_all_parameters(deep);
dr = resol(oo_.steady_state,0);
%deep(1:offset) = xparam1(1:offset);
logpo2(b,1) = DsgeLikelihood(deep,gend,data);
if opt_gsa.lik_only==0,
[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(deep,gend,data);
stock_smooth(:,:,ib)=atT(1:M_.endo_nbr,:);
% stock_filter(:,:,ib)=filtered_state_vector(1:M_.endo_nbr,:);
stock_filter(:,:,ib)=aK(1,1:M_.endo_nbr,:);
stock_ys(ib,:)=ys';
if ib==40,
ib=0;
ifil=ifil+1;
save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
stock_smooth = zeros(M_.endo_nbr,gend,40);
stock_filter = zeros(M_.endo_nbr,gend+1,40);
stock_ys = zeros(40, M_.endo_nbr);
end
end
waitbar(b/B,h,['MC smoother ...',num2str(b),'/',num2str(B)]);
end
close(h)
if opt_gsa.lik_only==0,
if ib>0,
ifil=ifil+1;
stock_smooth = stock_smooth(:,:,1:ib);
stock_filter = stock_filter(:,:,1:ib);
stock_ys = stock_ys(1:ib,:);
save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
end
end
stock_gend=gend;
stock_data=data;
save([OutDir,'\',namfile],'x','logpo2','stock_gend','stock_data','-append')
This diff is collapsed.
function cdf = gamm_cdf (x, a)
% PURPOSE: returns the cdf at x of the gamma(a) distribution
%---------------------------------------------------
% USAGE: cdf = gamm_cdf(x,a)
% where: x = a vector
% a = a scalar gamma(a)
%---------------------------------------------------
% RETURNS:
% a vector of cdf at each element of x of the gamma(a) distribution
% --------------------------------------------------
% SEE ALSO: gamm_d, gamm_pdf, gamm_rnd, gamm_inv
%---------------------------------------------------
% Anders Holtsberg, 18-11-93
% Copyright (c) Anders Holtsberg
if nargin ~= 2
error('Wrong # of arguments to gamm_cdf');
end;
if any(any(a<=0))
error('gamm_cdf: parameter a is wrong')
end
cdf = gammainc(x,a);
I0 = find(x<0);
cdf(I0) = zeros(size(I0));
function x = gamm_inv(p,a,b)
% PURPOSE: returns the inverse of the cdf at p of the gamma(a,b) distribution
%---------------------------------------------------
% USAGE: x = gamm_inv(p,a)
% where: p = a vector of probabilities
% a = a scalar parameter gamma(a)
% b = scaling factor for gamma
%---------------------------------------------------
% RETURNS:
% a vector x of the quantile at each element of p of the gamma(a) distribution
% --------------------------------------------------
% SEE ALSO: gamm_d, gamm_pdf, gamm_rnd, gamm_cdf
%---------------------------------------------------
% Anders Holtsberg, 18-11-93
% Copyright (c) Anders Holtsberg
% documentation modified by LeSage to fit the format
% of the econometrics toolbox
if (nargin < 2 | isempty(b))
b=1;
end
if (nargin > 3)
error('Wrong # of arguments to gamm_inv');
end
if any(any(abs(2*p-1)>1))
error('gamm_inv: a probability should be 0<=p<=1')
end
if any(any(a<=0))
error('gamma_inv: parameter a is wrong')
end
x = max(a-1,0.1);
dx = 1;
while any(any(abs(dx)>256*eps*max(x,1)))
dx = (gamm_cdf(x,a) - p) ./ gamm_pdf(x,a);
x = x - dx;
x = x + (dx - x) / 2 .* (x<0);
end
I0 = find(p==0);
x(I0) = zeros(size(I0));
I1 = find(p==1);
x(I1) = zeros(size(I0)) + Inf;
x=x.*b;
function f = gamm_pdf (x, a)
% PURPOSE: returns the pdf at x of the gamma(a) distribution
%---------------------------------------------------
% USAGE: pdf = gamm_pdf(x,a)
% where: x = a vector
% a = a scalar for gamma(a)
%---------------------------------------------------
% RETURNS:
% a vector of pdf at each element of x of the gamma(a) distribution
% --------------------------------------------------
% SEE ALSO: gamm_cdf, gamm_rnd, gamm_inv
%---------------------------------------------------
% Anders Holtsberg, 18-11-93
% Copyright (c) Anders Holtsberg
if nargin ~= 2
error('Wrong # of arguments to gamm_cdf');
end;
if any(any(a<=0))
error('gamm_pdf: parameter a is wrong')
end
f = x .^ (a-1) .* exp(-x) ./ gamma(a);
I0 = find(x<0);
f(I0) = zeros(size(I0));
function [A,B] = ghx2transition(mm,iv,ic,aux)
% [A,B] = ghx2transition(mm,iv,ic,aux)
%
% Adapted by M. Ratto from kalman_transition_matrix.m
% (kalman_transition_matrix.m is part of DYNARE, copyright M. Juillard)
%
% Part of the Sensitivity Analysis Toolbox for DYNARE
%
% Written by Marco Ratto, 2006
% Joint Research Centre, The European Commission,
% (http://eemc.jrc.ec.europa.eu/),
% marco.ratto@jrc.it
%
% Disclaimer: This software is not subject to copyright protection and is in the public domain.
% It is an experimental system. The Joint Research Centre of European Commission
% assumes no responsibility whatsoever for its use by other parties
% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
% characteristic. We would appreciate acknowledgement if the software is used.
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
%
global M_
[nr1, nc1] = size(mm);
ghx = mm(:, [1:(nc1-M_.exo_nbr)]);
ghu = mm(:, [(nc1-M_.exo_nbr+1):end] );
n_iv = length(iv);
n_ir1 = size(aux,1);
nr = n_iv + n_ir1;
A = zeros(nr,nr);
B = zeros(nr,M_.exo_nbr);
i_n_iv = 1:n_iv;
A(i_n_iv,ic) = ghx(iv,:);
if n_ir1 > 0
A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr);
end
B(i_n_iv,:) = ghu(iv,:);
function [yy, xdir, isig, lam]=log_trans_(y0,xdir0)
if nargin==1,
xdir0='';
end
f=inline('skewness(log(y+lam))','lam','y');
isig=1;
if ~(max(y0)<0 | min(y0)>0)
if skewness(y0)<0,
isig=-1;
y0=-y0;
end
n=hist(y0,10);
if n(1)>20*n(end),
try lam=fzero(f,[-min(y0)+10*eps -min(y0)+abs(median(y0))],[],y0);
catch
yl(1)=f(-min(y0)+10*eps,y0);
yl(2)=f(-min(y0)+abs(median(y0)),y0);
if abs(yl(1))<abs(yl(2))
lam=-min(y0)+eps;
else
lam = -min(y0)+abs(median(y0)); %abs(100*(1+min(y0)));
end
end
yy = log(y0+lam);
xdir=[xdir0,'_logskew'];
else
isig=0;
lam=0;
yy = log(y0.^2);
xdir=[xdir0,'_logsquared'];
end
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 lam=fzero(f,[-min(y0)+10*eps -min(y0)+median(y0)],[],y0);
catch
yl(1)=f(-min(y0)+10*eps,y0);
yl(2)=f(-min(y0)+abs(median(y0)),y0);
if abs(yl(1))<abs(yl(2))
lam=-min(y0)+eps;
else
lam = -min(y0)+abs(median(y0)); %abs(100*(1+min(y0)));
end
end
yy = log(y0+lam);
end
function [lpmat] = lptauSEQ(Nsam,Nvar)
% [lpmat] = lptauSEQ(Nsam,Nvar)
%
% function call LPTAU and generates a sample of dimension Nsam for a
% number of parameters Nvar
%
% Copyright (C) 2005 Marco Ratto
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
% Marco Ratto,
% Unit of Econometrics and Statistics AF
% (http://www.jrc.cec.eu.int/uasa/),
% IPSC, Joint Research Centre
% The European Commission,
% TP 361, 21020 ISPRA(VA), ITALY
% marco.ratto@jrc.it
%
clear lptau
lpmat = zeros(Nsam, Nvar);
for j=1:Nsam,
lpmat(j,:)=LPTAU(j,Nvar);
end
return
This diff is collapsed.
function [vdec, cc, ac] = mc_moments(mm, ss, dr)
global options_ M_
[nr1, nc1, nsam] = size(mm);
disp('Computing theoretical moments ...')
h = waitbar(0,'Theoretical moments ...');
for j=1:nsam,
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,:));
end
[vdec(:,:,j), corr, autocorr, z, zz] = th_moments(dr,options_.varobs);
cc(:,:,j)=tril(corr,-1);
ac(:,:,j)=autocorr{1};
waitbar(j/nsam,h)
end
close(h)
disp(' ')
disp('... done !')
function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
% % % % 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 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, MM+(MM-mm)*0.05]);
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment