gsa: converted DOS end of lines in Unix end of lines
parent
ec0af45fc8
commit
bc8d4d8f08
1024
matlab/gsa/LPTAU.m
1024
matlab/gsa/LPTAU.m
File diff suppressed because it is too large
Load Diff
|
@ -1,130 +1,130 @@
|
|||
|
||||
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 [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
|
||||
|
|
|
@ -1,174 +1,174 @@
|
|||
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'];
|
||||
|
||||
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
|
|
@ -1,14 +1,14 @@
|
|||
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 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
|
||||
|
|
|
@ -1,30 +1,30 @@
|
|||
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},';']);
|
||||
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
|
File diff suppressed because it is too large
Load Diff
|
@ -1,54 +1,54 @@
|
|||
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 oo_ M_
|
||||
|
||||
[nr1, nc1] = size(mm);
|
||||
ghx = mm(:, [1:(nc1-M_.exo_nbr)]);
|
||||
ghu = mm(:, [(nc1-M_.exo_nbr+1):end] );
|
||||
if nargin == 1
|
||||
oo_.dr.ghx = ghx;
|
||||
oo_.dr.ghu = ghu;
|
||||
endo_nbr = M_.endo_nbr;
|
||||
nstatic = oo_.dr.nstatic;
|
||||
npred = oo_.dr.npred;
|
||||
iv = (1:endo_nbr)';
|
||||
ic = [ nstatic+(1:npred) endo_nbr+(1:size(oo_.dr.ghx,2)-npred) ]';
|
||||
aux = oo_.dr.transition_auxiliary_variables;
|
||||
k = find(aux(:,2) > npred);
|
||||
aux(:,2) = aux(:,2) + nstatic;
|
||||
aux(k,2) = aux(k,2) + oo_.dr.nfwrd;
|
||||
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 [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 oo_ M_
|
||||
|
||||
[nr1, nc1] = size(mm);
|
||||
ghx = mm(:, [1:(nc1-M_.exo_nbr)]);
|
||||
ghu = mm(:, [(nc1-M_.exo_nbr+1):end] );
|
||||
if nargin == 1
|
||||
oo_.dr.ghx = ghx;
|
||||
oo_.dr.ghu = ghu;
|
||||
endo_nbr = M_.endo_nbr;
|
||||
nstatic = oo_.dr.nstatic;
|
||||
npred = oo_.dr.npred;
|
||||
iv = (1:endo_nbr)';
|
||||
ic = [ nstatic+(1:npred) endo_nbr+(1:size(oo_.dr.ghx,2)-npred) ]';
|
||||
aux = oo_.dr.transition_auxiliary_variables;
|
||||
k = find(aux(:,2) > npred);
|
||||
aux(:,2) = aux(:,2) + nstatic;
|
||||
aux(k,2) = aux(k,2) + oo_.dr.nfwrd;
|
||||
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,:);
|
||||
|
|
|
@ -1,94 +1,94 @@
|
|||
function gsa_plotmatrix(type,varargin)
|
||||
% function gsa_plotmatrix(type,varargin)
|
||||
% extended version of the standard MATLAB plotmatrix
|
||||
|
||||
% Copyright (C) 2011-2011 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
global bayestopt_ options_ M_
|
||||
|
||||
RootDirectoryName = CheckPath('gsa');
|
||||
|
||||
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
|
||||
function gsa_plotmatrix(type,varargin)
|
||||
% function gsa_plotmatrix(type,varargin)
|
||||
% extended version of the standard MATLAB plotmatrix
|
||||
|
||||
% Copyright (C) 2011-2011 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
global bayestopt_ options_ M_
|
||||
|
||||
RootDirectoryName = CheckPath('gsa');
|
||||
|
||||
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
|
||||
|
|
|
@ -1,54 +1,54 @@
|
|||
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 [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
|
||||
|
|
|
@ -1,25 +1,25 @@
|
|||
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
|
||||
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
|
||||
|
|
File diff suppressed because it is too large
Load Diff
|
@ -1,25 +1,25 @@
|
|||
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)=triu(corr);
|
||||
dum=[];
|
||||
for i=1:options_.ar
|
||||
dum=[dum, autocorr{i}];
|
||||
end
|
||||
ac(:,:,j)=dum;
|
||||
waitbar(j/nsam,h)
|
||||
end
|
||||
close(h)
|
||||
disp(' ')
|
||||
disp('... done !')
|
||||
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)=triu(corr);
|
||||
dum=[];
|
||||
for i=1:options_.ar
|
||||
dum=[dum, autocorr{i}];
|
||||
end
|
||||
ac(:,:,j)=dum;
|
||||
waitbar(j/nsam,h)
|
||||
end
|
||||
close(h)
|
||||
disp(' ')
|
||||
disp('... done !')
|
||||
|
|
|
@ -1,158 +1,158 @@
|
|||
|
||||
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-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
|
||||
|
||||
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-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
|
|
@ -1,20 +1,20 @@
|
|||
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];
|
||||
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);
|
|
@ -1,39 +1,39 @@
|
|||
function [gend, data] = read_data
|
||||
% 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 options_ bayestopt_
|
||||
|
||||
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
|
||||
function [gend, data] = read_data
|
||||
% 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 options_ bayestopt_
|
||||
|
||||
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
|
||||
|
|
|
@ -1,351 +1,351 @@
|
|||
function redform_map(dirname)
|
||||
%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;
|
||||
%
|
||||
% 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_ 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;
|
||||
pprior = options_gsa_.pprior;
|
||||
ilog = options_gsa_.logtrans_redform;
|
||||
threshold = options_gsa_.threshold_redform;
|
||||
ksstat = options_gsa_.ksstat_redform;
|
||||
alpha2 = options_gsa_.alpha2_redform;
|
||||
|
||||
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
|
||||
if nargin==0,
|
||||
dirname='';
|
||||
end
|
||||
|
||||
if pprior
|
||||
load([dirname,'/',M_.fname,'_prior']);
|
||||
adir=[dirname '/redform_stab'];
|
||||
else
|
||||
load([dirname,'/',M_.fname,'_mc']);
|
||||
adir=[dirname '/redform_mc'];
|
||||
end
|
||||
if ~exist('T')
|
||||
stab_map_(dirname);
|
||||
if pprior
|
||||
load([dirname,'/',M_.fname,'_prior'],'T');
|
||||
else
|
||||
load([dirname,'/',M_.fname,'_mc'],'T');
|
||||
end
|
||||
end
|
||||
if isempty(dir(adir))
|
||||
mkdir(adir)
|
||||
end
|
||||
adir0=pwd;
|
||||
%cd(adir)
|
||||
|
||||
nspred=size(T,2)-M_.exo_nbr;
|
||||
x0=lpmat(istable,:);
|
||||
[kn, np]=size(x0);
|
||||
offset = length(bayestopt_.pshape)-np;
|
||||
if options_gsa_.prior_range,
|
||||
pshape=5*(ones(np,1));
|
||||
pd = [NaN(np,1) NaN(np,1) bayestopt_.lb(offset+1:end) bayestopt_.ub(offset+1:end)];
|
||||
else
|
||||
pshape = bayestopt_.pshape(offset+1:end);
|
||||
pd = [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)];
|
||||
end
|
||||
|
||||
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
|
||||
clear lpmat lpmat0 egg iunstable yys
|
||||
js=0;
|
||||
for j=1:size(anamendo,1)
|
||||
namendo=deblank(anamendo(j,:));
|
||||
iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact');
|
||||
ifig=0;
|
||||
iplo=0;
|
||||
for jx=1:size(anamexo,1)
|
||||
namexo=deblank(anamexo(jx,:));
|
||||
iexo=strmatch(namexo,M_.exo_names,'exact');
|
||||
|
||||
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,
|
||||
ifig=ifig+1;
|
||||
hfig = figure('name',[namendo,' vs. shocks ',int2str(ifig)]);
|
||||
iplo=0;
|
||||
end
|
||||
iplo=iplo+1;
|
||||
js=js+1;
|
||||
xdir0 = [adir,'/',namendo,'_vs_', namexo];
|
||||
if ilog==0,
|
||||
if isempty(threshold)
|
||||
si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0);
|
||||
else
|
||||
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
|
||||
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
|
||||
xdir = [xdir0,'_cut'];
|
||||
if ~isempty(iy),
|
||||
si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir);
|
||||
end
|
||||
if ~isempty(iy) & ~isempty(iyc)
|
||||
delete([xdir, '/*cut*.*'])
|
||||
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
|
||||
indsmirnov = find(dproba>ksstat);
|
||||
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
|
||||
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
|
||||
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
|
||||
end
|
||||
end
|
||||
else
|
||||
[yy, xdir] = log_trans_(y0,xdir0);
|
||||
silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir);
|
||||
end
|
||||
|
||||
figure(hfig)
|
||||
subplot(3,3,iplo),
|
||||
if ilog,
|
||||
[saso, iso] = sort(-silog(:,js));
|
||||
bar([silog(iso(1:min(np,10)),js)])
|
||||
logflag='log';
|
||||
else
|
||||
[saso, 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),
|
||||
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
|
||||
end
|
||||
title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
ifig=0;
|
||||
iplo=0;
|
||||
for je=1:size(anamlagendo,1)
|
||||
namlagendo=deblank(anamlagendo(je,:));
|
||||
ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
|
||||
|
||||
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,
|
||||
ifig=ifig+1;
|
||||
hfig = figure('name',[namendo,' vs. lags ',int2str(ifig)]);
|
||||
iplo=0;
|
||||
end
|
||||
iplo=iplo+1;
|
||||
js=js+1;
|
||||
xdir0 = [adir,'/',namendo,'_vs_', namlagendo];
|
||||
if ilog==0,
|
||||
if isempty(threshold)
|
||||
si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0);
|
||||
else
|
||||
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
|
||||
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
|
||||
xdir = [xdir0,'_cut'];
|
||||
if ~isempty(iy)
|
||||
si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir);
|
||||
end
|
||||
if ~isempty(iy) & ~isempty(iyc),
|
||||
delete([xdir, '/*cut*.*'])
|
||||
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
|
||||
indsmirnov = find(dproba>ksstat);
|
||||
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
|
||||
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
|
||||
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
|
||||
end
|
||||
end
|
||||
else
|
||||
[yy, xdir] = log_trans_(y0,xdir0);
|
||||
silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir);
|
||||
end
|
||||
|
||||
figure(hfig),
|
||||
subplot(3,3,iplo),
|
||||
if ilog,
|
||||
[saso, iso] = sort(-silog(:,js));
|
||||
bar([silog(iso(1:min(np,10)),js)])
|
||||
logflag='log';
|
||||
else
|
||||
[saso, 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),
|
||||
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
|
||||
end
|
||||
title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
end
|
||||
|
||||
if ilog==0,
|
||||
figure, %bar(si)
|
||||
% boxplot(si','whis',10,'symbol','r.')
|
||||
myboxplot(si',[],'.',[],10)
|
||||
xlabel(' ')
|
||||
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
|
||||
title('Reduced form GSA')
|
||||
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_gsa'])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_gsa']);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_gsa']);
|
||||
|
||||
else
|
||||
figure, %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')
|
||||
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_gsa_log'])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_gsa_log']);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_gsa_log']);
|
||||
|
||||
end
|
||||
function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir)
|
||||
global bayestopt_ options_
|
||||
|
||||
opt_gsa=options_.opt_gsa;
|
||||
np=size(x0,2);
|
||||
x00=x0;
|
||||
if opt_gsa.prior_range,
|
||||
for j=1:np,
|
||||
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));
|
||||
end
|
||||
|
||||
fname=[xdir,'/map'];
|
||||
if iload==0,
|
||||
figure, hist(y0,30), title([namy,' vs. ', namx])
|
||||
if isempty(dir(xdir))
|
||||
mkdir(xdir)
|
||||
end
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx]);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx]);
|
||||
close(gcf)
|
||||
% gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames);
|
||||
nrun=length(y0);
|
||||
nest=min(250,nrun);
|
||||
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);
|
||||
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);
|
||||
end
|
||||
save([fname,'.mat'],'gsa_')
|
||||
[sidum, iii]=sort(-gsa_.si);
|
||||
gsa_.x0=x00(1:nfit,:);
|
||||
hfig=gsa_sdp_plot(gsa_,fname,pnames,iii(1:min(12,np)));
|
||||
close(hfig);
|
||||
gsa_.x0=x0(1:nfit,:);
|
||||
% copyfile([fname,'_est.mat'],[fname,'.mat'])
|
||||
figure,
|
||||
plot(y0(1:nfit),[gsa_.fit y0(1:nfit)],'.'),
|
||||
title([namy,' vs. ', namx,' fit'])
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx,'_fit'])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx,'_fit']);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx,'_fit']);
|
||||
close(gcf)
|
||||
if nfit<nrun,
|
||||
npred=[nfit+1:nrun];
|
||||
yf = ss_anova_fcast(x0(npred,:), gsa_);
|
||||
figure,
|
||||
plot(y0(npred),[yf y0(npred)],'.'),
|
||||
title([namy,' vs. ', namx,' pred'])
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx,'_pred'])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
close(gcf)
|
||||
end
|
||||
else
|
||||
% gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames);
|
||||
gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames);
|
||||
yf = ss_anova_fcast(x0, gsa_);
|
||||
figure,
|
||||
plot(y0,[yf y0],'.'),
|
||||
title([namy,' vs. ', namx,' pred'])
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx,'_pred'])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
close(gcf)
|
||||
end
|
||||
% si = gsa_.multivariate.si;
|
||||
si = gsa_.si;
|
||||
|
||||
function redform_map(dirname)
|
||||
%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;
|
||||
%
|
||||
% 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_ 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;
|
||||
pprior = options_gsa_.pprior;
|
||||
ilog = options_gsa_.logtrans_redform;
|
||||
threshold = options_gsa_.threshold_redform;
|
||||
ksstat = options_gsa_.ksstat_redform;
|
||||
alpha2 = options_gsa_.alpha2_redform;
|
||||
|
||||
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
|
||||
if nargin==0,
|
||||
dirname='';
|
||||
end
|
||||
|
||||
if pprior
|
||||
load([dirname,'/',M_.fname,'_prior']);
|
||||
adir=[dirname '/redform_stab'];
|
||||
else
|
||||
load([dirname,'/',M_.fname,'_mc']);
|
||||
adir=[dirname '/redform_mc'];
|
||||
end
|
||||
if ~exist('T')
|
||||
stab_map_(dirname);
|
||||
if pprior
|
||||
load([dirname,'/',M_.fname,'_prior'],'T');
|
||||
else
|
||||
load([dirname,'/',M_.fname,'_mc'],'T');
|
||||
end
|
||||
end
|
||||
if isempty(dir(adir))
|
||||
mkdir(adir)
|
||||
end
|
||||
adir0=pwd;
|
||||
%cd(adir)
|
||||
|
||||
nspred=size(T,2)-M_.exo_nbr;
|
||||
x0=lpmat(istable,:);
|
||||
[kn, np]=size(x0);
|
||||
offset = length(bayestopt_.pshape)-np;
|
||||
if options_gsa_.prior_range,
|
||||
pshape=5*(ones(np,1));
|
||||
pd = [NaN(np,1) NaN(np,1) bayestopt_.lb(offset+1:end) bayestopt_.ub(offset+1:end)];
|
||||
else
|
||||
pshape = bayestopt_.pshape(offset+1:end);
|
||||
pd = [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)];
|
||||
end
|
||||
|
||||
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
|
||||
clear lpmat lpmat0 egg iunstable yys
|
||||
js=0;
|
||||
for j=1:size(anamendo,1)
|
||||
namendo=deblank(anamendo(j,:));
|
||||
iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact');
|
||||
ifig=0;
|
||||
iplo=0;
|
||||
for jx=1:size(anamexo,1)
|
||||
namexo=deblank(anamexo(jx,:));
|
||||
iexo=strmatch(namexo,M_.exo_names,'exact');
|
||||
|
||||
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,
|
||||
ifig=ifig+1;
|
||||
hfig = figure('name',[namendo,' vs. shocks ',int2str(ifig)]);
|
||||
iplo=0;
|
||||
end
|
||||
iplo=iplo+1;
|
||||
js=js+1;
|
||||
xdir0 = [adir,'/',namendo,'_vs_', namexo];
|
||||
if ilog==0,
|
||||
if isempty(threshold)
|
||||
si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0);
|
||||
else
|
||||
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
|
||||
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
|
||||
xdir = [xdir0,'_cut'];
|
||||
if ~isempty(iy),
|
||||
si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir);
|
||||
end
|
||||
if ~isempty(iy) & ~isempty(iyc)
|
||||
delete([xdir, '/*cut*.*'])
|
||||
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
|
||||
indsmirnov = find(dproba>ksstat);
|
||||
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
|
||||
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
|
||||
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
|
||||
end
|
||||
end
|
||||
else
|
||||
[yy, xdir] = log_trans_(y0,xdir0);
|
||||
silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir);
|
||||
end
|
||||
|
||||
figure(hfig)
|
||||
subplot(3,3,iplo),
|
||||
if ilog,
|
||||
[saso, iso] = sort(-silog(:,js));
|
||||
bar([silog(iso(1:min(np,10)),js)])
|
||||
logflag='log';
|
||||
else
|
||||
[saso, 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),
|
||||
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
|
||||
end
|
||||
title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
ifig=0;
|
||||
iplo=0;
|
||||
for je=1:size(anamlagendo,1)
|
||||
namlagendo=deblank(anamlagendo(je,:));
|
||||
ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
|
||||
|
||||
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,
|
||||
ifig=ifig+1;
|
||||
hfig = figure('name',[namendo,' vs. lags ',int2str(ifig)]);
|
||||
iplo=0;
|
||||
end
|
||||
iplo=iplo+1;
|
||||
js=js+1;
|
||||
xdir0 = [adir,'/',namendo,'_vs_', namlagendo];
|
||||
if ilog==0,
|
||||
if isempty(threshold)
|
||||
si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0);
|
||||
else
|
||||
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
|
||||
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
|
||||
xdir = [xdir0,'_cut'];
|
||||
if ~isempty(iy)
|
||||
si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir);
|
||||
end
|
||||
if ~isempty(iy) & ~isempty(iyc),
|
||||
delete([xdir, '/*cut*.*'])
|
||||
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
|
||||
indsmirnov = find(dproba>ksstat);
|
||||
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
|
||||
stab_map_2(x0(iy,:),alpha2,'cut',xdir)
|
||||
stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
|
||||
end
|
||||
end
|
||||
else
|
||||
[yy, xdir] = log_trans_(y0,xdir0);
|
||||
silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir);
|
||||
end
|
||||
|
||||
figure(hfig),
|
||||
subplot(3,3,iplo),
|
||||
if ilog,
|
||||
[saso, iso] = sort(-silog(:,js));
|
||||
bar([silog(iso(1:min(np,10)),js)])
|
||||
logflag='log';
|
||||
else
|
||||
[saso, 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),
|
||||
text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
|
||||
end
|
||||
title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
end
|
||||
|
||||
if ilog==0,
|
||||
figure, %bar(si)
|
||||
% boxplot(si','whis',10,'symbol','r.')
|
||||
myboxplot(si',[],'.',[],10)
|
||||
xlabel(' ')
|
||||
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
|
||||
title('Reduced form GSA')
|
||||
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_gsa'])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_gsa']);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_gsa']);
|
||||
|
||||
else
|
||||
figure, %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')
|
||||
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_gsa_log'])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_gsa_log']);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_gsa_log']);
|
||||
|
||||
end
|
||||
function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir)
|
||||
global bayestopt_ options_
|
||||
|
||||
opt_gsa=options_.opt_gsa;
|
||||
np=size(x0,2);
|
||||
x00=x0;
|
||||
if opt_gsa.prior_range,
|
||||
for j=1:np,
|
||||
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));
|
||||
end
|
||||
|
||||
fname=[xdir,'/map'];
|
||||
if iload==0,
|
||||
figure, hist(y0,30), title([namy,' vs. ', namx])
|
||||
if isempty(dir(xdir))
|
||||
mkdir(xdir)
|
||||
end
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx]);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx]);
|
||||
close(gcf)
|
||||
% gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames);
|
||||
nrun=length(y0);
|
||||
nest=min(250,nrun);
|
||||
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);
|
||||
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);
|
||||
end
|
||||
save([fname,'.mat'],'gsa_')
|
||||
[sidum, iii]=sort(-gsa_.si);
|
||||
gsa_.x0=x00(1:nfit,:);
|
||||
hfig=gsa_sdp_plot(gsa_,fname,pnames,iii(1:min(12,np)));
|
||||
close(hfig);
|
||||
gsa_.x0=x0(1:nfit,:);
|
||||
% copyfile([fname,'_est.mat'],[fname,'.mat'])
|
||||
figure,
|
||||
plot(y0(1:nfit),[gsa_.fit y0(1:nfit)],'.'),
|
||||
title([namy,' vs. ', namx,' fit'])
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx,'_fit'])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx,'_fit']);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx,'_fit']);
|
||||
close(gcf)
|
||||
if nfit<nrun,
|
||||
npred=[nfit+1:nrun];
|
||||
yf = ss_anova_fcast(x0(npred,:), gsa_);
|
||||
figure,
|
||||
plot(y0(npred),[yf y0(npred)],'.'),
|
||||
title([namy,' vs. ', namx,' pred'])
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx,'_pred'])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
close(gcf)
|
||||
end
|
||||
else
|
||||
% gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames);
|
||||
gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames);
|
||||
yf = ss_anova_fcast(x0, gsa_);
|
||||
figure,
|
||||
plot(y0,[yf y0],'.'),
|
||||
title([namy,' vs. ', namx,' pred'])
|
||||
saveas(gcf,[xdir,'/', namy,'_vs_', namx,'_pred'])
|
||||
eval(['print -depsc2 ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
eval(['print -dpdf ' xdir,'/', namy,'_vs_', namx,'_pred']);
|
||||
close(gcf)
|
||||
end
|
||||
% si = gsa_.multivariate.si;
|
||||
si = gsa_.si;
|
||||
|
||||
|
|
|
@ -1,165 +1,165 @@
|
|||
function redform_screen(dirname)
|
||||
%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;
|
||||
%
|
||||
% 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_ 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;
|
||||
nliv = options_gsa_.morris_nliv;
|
||||
|
||||
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
|
||||
if nargin==0,
|
||||
dirname='';
|
||||
end
|
||||
|
||||
load([dirname,'/',M_.fname,'_prior'],'lpmat','lpmat0','istable','T');
|
||||
|
||||
nspred=oo_.dr.nspred;
|
||||
|
||||
[kn, np]=size(lpmat);
|
||||
nshock = length(bayestopt_.pshape)-np;
|
||||
|
||||
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');
|
||||
|
||||
iplo=0;
|
||||
ifig=0;
|
||||
for jx=1:size(anamexo,1)
|
||||
namexo=deblank(anamexo(jx,:));
|
||||
iexo=strmatch(namexo,M_.exo_names,'exact');
|
||||
|
||||
if ~isempty(iexo),
|
||||
y0=teff(T(iendo,iexo+nspred,:),kn,istable);
|
||||
if ~isempty(y0),
|
||||
if mod(iplo,9)==0,
|
||||
ifig=ifig+1;
|
||||
figure('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);
|
||||
SAM = squeeze(SAMorris(nshock+1:end,1));
|
||||
SA(:,js)=SAM./(max(SAM)+eps);
|
||||
[saso, 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')
|
||||
end
|
||||
title([namendo,' vs. ',namexo],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
iplo=0;
|
||||
ifig=0;
|
||||
for je=1:size(anamlagendo,1)
|
||||
namlagendo=deblank(anamlagendo(je,:));
|
||||
ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
|
||||
|
||||
if ~isempty(ilagendo),
|
||||
y0=teff(T(iendo,ilagendo,:),kn,istable);
|
||||
if ~isempty(y0),
|
||||
if mod(iplo,9)==0,
|
||||
ifig=ifig+1;
|
||||
figure('name',[namendo,' vs. lagged endogenous ',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);
|
||||
SAM = squeeze(SAMorris(nshock+1:end,1));
|
||||
SA(:,js)=SAM./(max(SAM)+eps);
|
||||
[saso, 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')
|
||||
end
|
||||
|
||||
title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
end
|
||||
|
||||
figure,
|
||||
%bar(SA)
|
||||
% boxplot(SA','whis',10,'symbol','r.')
|
||||
myboxplot(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')
|
||||
end
|
||||
xlabel(' ')
|
||||
ylabel('Elementary Effects')
|
||||
title('Reduced form screening')
|
||||
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_screen'])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_screen']);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_screen']);
|
||||
|
||||
function redform_screen(dirname)
|
||||
%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;
|
||||
%
|
||||
% 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_ 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;
|
||||
nliv = options_gsa_.morris_nliv;
|
||||
|
||||
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
|
||||
if nargin==0,
|
||||
dirname='';
|
||||
end
|
||||
|
||||
load([dirname,'/',M_.fname,'_prior'],'lpmat','lpmat0','istable','T');
|
||||
|
||||
nspred=oo_.dr.nspred;
|
||||
|
||||
[kn, np]=size(lpmat);
|
||||
nshock = length(bayestopt_.pshape)-np;
|
||||
|
||||
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');
|
||||
|
||||
iplo=0;
|
||||
ifig=0;
|
||||
for jx=1:size(anamexo,1)
|
||||
namexo=deblank(anamexo(jx,:));
|
||||
iexo=strmatch(namexo,M_.exo_names,'exact');
|
||||
|
||||
if ~isempty(iexo),
|
||||
y0=teff(T(iendo,iexo+nspred,:),kn,istable);
|
||||
if ~isempty(y0),
|
||||
if mod(iplo,9)==0,
|
||||
ifig=ifig+1;
|
||||
figure('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);
|
||||
SAM = squeeze(SAMorris(nshock+1:end,1));
|
||||
SA(:,js)=SAM./(max(SAM)+eps);
|
||||
[saso, 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')
|
||||
end
|
||||
title([namendo,' vs. ',namexo],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_shock_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
|
||||
iplo=0;
|
||||
ifig=0;
|
||||
for je=1:size(anamlagendo,1)
|
||||
namlagendo=deblank(anamlagendo(je,:));
|
||||
ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
|
||||
|
||||
if ~isempty(ilagendo),
|
||||
y0=teff(T(iendo,ilagendo,:),kn,istable);
|
||||
if ~isempty(y0),
|
||||
if mod(iplo,9)==0,
|
||||
ifig=ifig+1;
|
||||
figure('name',[namendo,' vs. lagged endogenous ',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);
|
||||
SAM = squeeze(SAMorris(nshock+1:end,1));
|
||||
SA(:,js)=SAM./(max(SAM)+eps);
|
||||
[saso, 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')
|
||||
end
|
||||
|
||||
title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
|
||||
if iplo==9,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
if iplo<9 & iplo>0 & ifig,
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]);
|
||||
close(gcf)
|
||||
end
|
||||
end
|
||||
|
||||
figure,
|
||||
%bar(SA)
|
||||
% boxplot(SA','whis',10,'symbol','r.')
|
||||
myboxplot(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')
|
||||
end
|
||||
xlabel(' ')
|
||||
ylabel('Elementary Effects')
|
||||
title('Reduced form screening')
|
||||
|
||||
saveas(gcf,[dirname,'/',M_.fname,'_redform_screen'])
|
||||
eval(['print -depsc2 ' dirname,'/',M_.fname,'_redform_screen']);
|
||||
eval(['print -dpdf ' dirname,'/',M_.fname,'_redform_screen']);
|
||||
|
||||
|
|
|
@ -1,30 +1,30 @@
|
|||
function set_shocks_param(xparam1)
|
||||
global estim_params_ M_
|
||||
|
||||
nvx = estim_params_.nvx;
|
||||
ncx = estim_params_.ncx;
|
||||
np = estim_params_.np;
|
||||
Sigma_e = M_.Sigma_e;
|
||||
offset = 0;
|
||||
if nvx
|
||||
offset = offset + nvx;
|
||||
var_exo = estim_params_.var_exo;
|
||||
for i=1:nvx
|
||||
k = var_exo(i,1);
|
||||
Sigma_e(k,k) = xparam1(i)^2;
|
||||
end
|
||||
end
|
||||
|
||||
if ncx
|
||||
offset = offset + estim_params_.nvn;
|
||||
corrx = estim_params_.corrx;
|
||||
for i=1:ncx
|
||||
k1 = corrx(i,1);
|
||||
k2 = corrx(i,2);
|
||||
Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e_(k1,k1)*Sigma_e_(k2,k2));
|
||||
Sigma_e(k2,k1) = Sigma_e_(k1,k2);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
function set_shocks_param(xparam1)
|
||||
global estim_params_ M_
|
||||
|
||||
nvx = estim_params_.nvx;
|
||||
ncx = estim_params_.ncx;
|
||||
np = estim_params_.np;
|
||||
Sigma_e = M_.Sigma_e;
|
||||
offset = 0;
|
||||
if nvx
|
||||
offset = offset + nvx;
|
||||
var_exo = estim_params_.var_exo;
|
||||
for i=1:nvx
|
||||
k = var_exo(i,1);
|
||||
Sigma_e(k,k) = xparam1(i)^2;
|
||||
end
|
||||
end
|
||||
|
||||
if ncx
|
||||
offset = offset + estim_params_.nvn;
|
||||
corrx = estim_params_.corrx;
|
||||
for i=1:ncx
|
||||
k1 = corrx(i,1);
|
||||
k2 = corrx(i,2);
|
||||
Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e_(k1,k1)*Sigma_e_(k2,k2));
|
||||
Sigma_e(k2,k1) = Sigma_e_(k1,k2);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
M_.Sigma_e = Sigma_e;
|
|
@ -1,7 +1,7 @@
|
|||
function s=skewness(y),
|
||||
|
||||
% y=stand_(y);
|
||||
% s=mean(y.^3);
|
||||
m2=mean((y-mean(y)).^2);
|
||||
m3=mean((y-mean(y)).^3);
|
||||
function s=skewness(y),
|
||||
|
||||
% y=stand_(y);
|
||||
% s=mean(y.^3);
|
||||
m2=mean((y-mean(y)).^2);
|
||||
m3=mean((y-mean(y)).^3);
|
||||
s=m3/m2^1.5;
|
|
@ -1,73 +1,73 @@
|
|||
function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
|
||||
% Smirnov test for 2 distributions
|
||||
% [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
|
||||
|
||||
if nargin<3
|
||||
alpha = 0.05;
|
||||
end
|
||||
if nargin<4,
|
||||
iflag=0;
|
||||
end
|
||||
|
||||
% empirical cdfs.
|
||||
xmix= [x1;x2];
|
||||
bin = [-inf ; sort(xmix) ; inf];
|
||||
|
||||
ncount1 = histc (x1 , bin);
|
||||
ncount1 = ncount1(:);
|
||||
ncount2 = histc (x2 , bin);
|
||||
ncount2 = ncount2(:);
|
||||
|
||||
cum1 = cumsum(ncount1)./sum(ncount1);
|
||||
cum1 = cum1(1:end-1);
|
||||
|
||||
cum2 = cumsum(ncount2)./sum(ncount2);
|
||||
cum2 = cum2(1:end-1);
|
||||
|
||||
n1= length(x1);
|
||||
n2= length(x2);
|
||||
n = n1*n2 /(n1+n2);
|
||||
|
||||
% Compute the d(n1,n2) statistics.
|
||||
|
||||
if iflag==0,
|
||||
d = max(abs(cum1 - cum2));
|
||||
elseif iflag==-1
|
||||
d = max(cum2 - cum1);
|
||||
elseif iflag==1
|
||||
d = max(cum1 - cum2);
|
||||
end
|
||||
%
|
||||
% Compute P-value check H0 hypothesis
|
||||
%
|
||||
|
||||
lam = max((sqrt(n) + 0.12 + 0.11/sqrt(n)) * d , 0);
|
||||
if iflag == 0
|
||||
j = [1:101]';
|
||||
prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2));
|
||||
|
||||
prob=max(prob,0);
|
||||
prob=min(1,prob);
|
||||
else
|
||||
prob = exp(-2*lam*lam);
|
||||
end
|
||||
|
||||
H = (alpha >= prob);
|
||||
function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
|
||||
% Smirnov test for 2 distributions
|
||||
% [H,prob,d] = smirnov(x1 , x2 , alpha, iflag )
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
|
||||
|
||||
if nargin<3
|
||||
alpha = 0.05;
|
||||
end
|
||||
if nargin<4,
|
||||
iflag=0;
|
||||
end
|
||||
|
||||
% empirical cdfs.
|
||||
xmix= [x1;x2];
|
||||
bin = [-inf ; sort(xmix) ; inf];
|
||||
|
||||
ncount1 = histc (x1 , bin);
|
||||
ncount1 = ncount1(:);
|
||||
ncount2 = histc (x2 , bin);
|
||||
ncount2 = ncount2(:);
|
||||
|
||||
cum1 = cumsum(ncount1)./sum(ncount1);
|
||||
cum1 = cum1(1:end-1);
|
||||
|
||||
cum2 = cumsum(ncount2)./sum(ncount2);
|
||||
cum2 = cum2(1:end-1);
|
||||
|
||||
n1= length(x1);
|
||||
n2= length(x2);
|
||||
n = n1*n2 /(n1+n2);
|
||||
|
||||
% Compute the d(n1,n2) statistics.
|
||||
|
||||
if iflag==0,
|
||||
d = max(abs(cum1 - cum2));
|
||||
elseif iflag==-1
|
||||
d = max(cum2 - cum1);
|
||||
elseif iflag==1
|
||||
d = max(cum1 - cum2);
|
||||
end
|
||||
%
|
||||
% Compute P-value check H0 hypothesis
|
||||
%
|
||||
|
||||
lam = max((sqrt(n) + 0.12 + 0.11/sqrt(n)) * d , 0);
|
||||
if iflag == 0
|
||||
j = [1:101]';
|
||||
prob = 2 * sum((-1).^(j-1).*exp(-2*lam*lam*j.^2));
|
||||
|
||||
prob=max(prob,0);
|
||||
prob=min(1,prob);
|
||||
else
|
||||
prob = exp(-2*lam*lam);
|
||||
end
|
||||
|
||||
H = (alpha >= prob);
|
||||
|
|
|
@ -1,52 +1,52 @@
|
|||
function [tadj, iff] = speed(A,B,mf,p),
|
||||
% [tadj, iff] = speed(A,B,mf,p),
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
nvar=length(mf);
|
||||
nstate= size(A,1);
|
||||
nshock = size(B,2);
|
||||
nrun=size(B,3);
|
||||
|
||||
iff=zeros(nvar,nshock,nrun);
|
||||
tadj=iff;
|
||||
disp('Computing speed of adjustement ...')
|
||||
h = waitbar(0,'Speed of adjustement...');
|
||||
|
||||
for i=1:nrun,
|
||||
irf=zeros(nvar,nshock);
|
||||
a=squeeze(A(:,:,i));
|
||||
b=squeeze(B(:,:,i));
|
||||
IFF=inv(eye(nstate)-a)*b;
|
||||
iff(:,:,i)=IFF(mf,:);
|
||||
IF=IFF-b;
|
||||
|
||||
t=0;
|
||||
while any(any(irf<0.5))
|
||||
t=t+1;
|
||||
IFT=((eye(nstate)-a^(t+1))*inv(eye(nstate)-a))*b-b;
|
||||
irf=IFT(mf,:)./(IF(mf,:)+eps);
|
||||
irf = irf.*(abs(IF(mf,:))>1.e-7)+(abs(IF(mf,:))<=1.e-7);
|
||||
%irf=ft(mf,:);
|
||||
tt=(irf>0.5).*t;
|
||||
tadj(:,:,i)=((tt-tadj(:,:,i))==tt).*tt+tadj(:,:,i);
|
||||
end
|
||||
waitbar(i/nrun,h)
|
||||
end
|
||||
disp(' ')
|
||||
disp('.. done !')
|
||||
close(h)
|
||||
function [tadj, iff] = speed(A,B,mf,p),
|
||||
% [tadj, iff] = speed(A,B,mf,p),
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
nvar=length(mf);
|
||||
nstate= size(A,1);
|
||||
nshock = size(B,2);
|
||||
nrun=size(B,3);
|
||||
|
||||
iff=zeros(nvar,nshock,nrun);
|
||||
tadj=iff;
|
||||
disp('Computing speed of adjustement ...')
|
||||
h = waitbar(0,'Speed of adjustement...');
|
||||
|
||||
for i=1:nrun,
|
||||
irf=zeros(nvar,nshock);
|
||||
a=squeeze(A(:,:,i));
|
||||
b=squeeze(B(:,:,i));
|
||||
IFF=inv(eye(nstate)-a)*b;
|
||||
iff(:,:,i)=IFF(mf,:);
|
||||
IF=IFF-b;
|
||||
|
||||
t=0;
|
||||
while any(any(irf<0.5))
|
||||
t=t+1;
|
||||
IFT=((eye(nstate)-a^(t+1))*inv(eye(nstate)-a))*b-b;
|
||||
irf=IFT(mf,:)./(IF(mf,:)+eps);
|
||||
irf = irf.*(abs(IF(mf,:))>1.e-7)+(abs(IF(mf,:))<=1.e-7);
|
||||
%irf=ft(mf,:);
|
||||
tt=(irf>0.5).*t;
|
||||
tadj(:,:,i)=((tt-tadj(:,:,i))==tt).*tt+tadj(:,:,i);
|
||||
end
|
||||
waitbar(i/nrun,h)
|
||||
end
|
||||
disp(' ')
|
||||
disp('.. done !')
|
||||
close(h)
|
||||
|
|
File diff suppressed because it is too large
Load Diff
|
@ -1,98 +1,98 @@
|
|||
function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname, pcrit)
|
||||
%function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname, pcrit)
|
||||
%
|
||||
% lpmat = Monte Carlo matrix
|
||||
% ibehaviour = index of behavioural runs
|
||||
% inonbehaviour = index of non-behavioural runs
|
||||
% aname = label of the analysis
|
||||
% iplot = 1 plot cumulative distributions (default)
|
||||
% iplot = 0 no plots
|
||||
% ipar = index array of parameters to plot
|
||||
% dirname = (OPTIONAL) path of the directory where to save
|
||||
% (default: current directory)
|
||||
% pcrit = (OPTIONAL) critical value of the pvalue below which show the plots
|
||||
%
|
||||
% Plots: dotted lines for BEHAVIOURAL
|
||||
% solid lines for NON BEHAVIOURAL
|
||||
% USES smirnov
|
||||
%
|
||||
% 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 estim_params_ bayestopt_ M_ options_
|
||||
|
||||
if nargin<5,
|
||||
iplot=1;
|
||||
end
|
||||
fname_ = M_.fname;
|
||||
if nargin<7,
|
||||
dirname='';;
|
||||
end
|
||||
|
||||
nshock = estim_params_.nvx;
|
||||
nshock = nshock + estim_params_.nvn;
|
||||
nshock = nshock + estim_params_.ncx;
|
||||
nshock = nshock + estim_params_.ncn;
|
||||
|
||||
npar=size(lpmat,2);
|
||||
ishock= npar>estim_params_.np;
|
||||
|
||||
if nargin<6,
|
||||
ipar=[];
|
||||
end
|
||||
if nargin<8 || isempty(pcrit),
|
||||
pcrit=1;
|
||||
end
|
||||
|
||||
% Smirnov test for Blanchard;
|
||||
for j=1:npar,
|
||||
[H,P,KSSTAT] = smirnov(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
|
||||
proba(j)=P;
|
||||
dproba(j)=KSSTAT;
|
||||
end
|
||||
if isempty(ipar),
|
||||
% ipar=find(dproba>dcrit);
|
||||
ipar=find(proba<pcrit);
|
||||
end
|
||||
nparplot=length(ipar);
|
||||
if iplot
|
||||
lpmat=lpmat(:,ipar);
|
||||
ftit=bayestopt_.name(ipar+nshock*(1-ishock));
|
||||
|
||||
for i=1:ceil(nparplot/12),
|
||||
hh=figure('name',aname);
|
||||
for j=1+12*(i-1):min(nparplot,12*i),
|
||||
subplot(3,4,j-12*(i-1))
|
||||
if ~isempty(ibehaviour),
|
||||
h=cumplot(lpmat(ibehaviour,j));
|
||||
set(h,'color',[0 0 0], 'linestyle',':')
|
||||
end
|
||||
hold on,
|
||||
if ~isempty(inonbehaviour),
|
||||
h=cumplot(lpmat(inonbehaviour,j));
|
||||
set(h,'color',[0 0 0])
|
||||
end
|
||||
% title([ftit{j},'. D-stat ', num2str(dproba(ipar(j)),2)],'interpreter','none')
|
||||
title([ftit{j},'. p-value ', num2str(proba(ipar(j)),2)],'interpreter','none')
|
||||
end
|
||||
eval(['print -depsc2 ' dirname '/' fname_ '_' aname '_SA_' int2str(i) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION'),
|
||||
saveas(hh,[dirname,'/',fname_,'_',aname,'_SA_',int2str(i)])
|
||||
eval(['print -dpdf ' dirname '/' fname_ '_' aname '_SA_' int2str(i)]);
|
||||
end
|
||||
if options_.nograph, close(hh), end
|
||||
end
|
||||
end
|
||||
function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname, pcrit)
|
||||
%function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, iplot, ipar, dirname, pcrit)
|
||||
%
|
||||
% lpmat = Monte Carlo matrix
|
||||
% ibehaviour = index of behavioural runs
|
||||
% inonbehaviour = index of non-behavioural runs
|
||||
% aname = label of the analysis
|
||||
% iplot = 1 plot cumulative distributions (default)
|
||||
% iplot = 0 no plots
|
||||
% ipar = index array of parameters to plot
|
||||
% dirname = (OPTIONAL) path of the directory where to save
|
||||
% (default: current directory)
|
||||
% pcrit = (OPTIONAL) critical value of the pvalue below which show the plots
|
||||
%
|
||||
% Plots: dotted lines for BEHAVIOURAL
|
||||
% solid lines for NON BEHAVIOURAL
|
||||
% USES smirnov
|
||||
%
|
||||
% 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 estim_params_ bayestopt_ M_ options_
|
||||
|
||||
if nargin<5,
|
||||
iplot=1;
|
||||
end
|
||||
fname_ = M_.fname;
|
||||
if nargin<7,
|
||||
dirname='';;
|
||||
end
|
||||
|
||||
nshock = estim_params_.nvx;
|
||||
nshock = nshock + estim_params_.nvn;
|
||||
nshock = nshock + estim_params_.ncx;
|
||||
nshock = nshock + estim_params_.ncn;
|
||||
|
||||
npar=size(lpmat,2);
|
||||
ishock= npar>estim_params_.np;
|
||||
|
||||
if nargin<6,
|
||||
ipar=[];
|
||||
end
|
||||
if nargin<8 || isempty(pcrit),
|
||||
pcrit=1;
|
||||
end
|
||||
|
||||
% Smirnov test for Blanchard;
|
||||
for j=1:npar,
|
||||
[H,P,KSSTAT] = smirnov(lpmat(ibehaviour,j),lpmat(inonbehaviour,j));
|
||||
proba(j)=P;
|
||||
dproba(j)=KSSTAT;
|
||||
end
|
||||
if isempty(ipar),
|
||||
% ipar=find(dproba>dcrit);
|
||||
ipar=find(proba<pcrit);
|
||||
end
|
||||
nparplot=length(ipar);
|
||||
if iplot
|
||||
lpmat=lpmat(:,ipar);
|
||||
ftit=bayestopt_.name(ipar+nshock*(1-ishock));
|
||||
|
||||
for i=1:ceil(nparplot/12),
|
||||
hh=figure('name',aname);
|
||||
for j=1+12*(i-1):min(nparplot,12*i),
|
||||
subplot(3,4,j-12*(i-1))
|
||||
if ~isempty(ibehaviour),
|
||||
h=cumplot(lpmat(ibehaviour,j));
|
||||
set(h,'color',[0 0 0], 'linestyle',':')
|
||||
end
|
||||
hold on,
|
||||
if ~isempty(inonbehaviour),
|
||||
h=cumplot(lpmat(inonbehaviour,j));
|
||||
set(h,'color',[0 0 0])
|
||||
end
|
||||
% title([ftit{j},'. D-stat ', num2str(dproba(ipar(j)),2)],'interpreter','none')
|
||||
title([ftit{j},'. p-value ', num2str(proba(ipar(j)),2)],'interpreter','none')
|
||||
end
|
||||
eval(['print -depsc2 ' dirname '/' fname_ '_' aname '_SA_' int2str(i) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION'),
|
||||
saveas(hh,[dirname,'/',fname_,'_',aname,'_SA_',int2str(i)])
|
||||
eval(['print -dpdf ' dirname '/' fname_ '_' aname '_SA_' int2str(i)]);
|
||||
end
|
||||
if options_.nograph, close(hh), end
|
||||
end
|
||||
end
|
||||
|
|
|
@ -1,105 +1,105 @@
|
|||
function stab_map_2(x,alpha2, pvalue, fnam, dirname)
|
||||
% function stab_map_2(x, alpha2, pvalue, fnam, dirname)
|
||||
%
|
||||
% 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 bayestopt_ estim_params_ dr_ options_ ys_ fname_
|
||||
global bayestopt_ estim_params_ options_ oo_ M_
|
||||
|
||||
npar=size(x,2);
|
||||
nsam=size(x,1);
|
||||
ishock= npar>estim_params_.np;
|
||||
if nargin<4,
|
||||
fnam='';
|
||||
end
|
||||
if nargin<5,
|
||||
dirname='';
|
||||
end
|
||||
|
||||
ys_ = oo_.dr.ys;
|
||||
dr_ = oo_.dr;
|
||||
fname_ = M_.fname;
|
||||
nshock = estim_params_.nvx;
|
||||
nshock = nshock + estim_params_.nvn;
|
||||
nshock = nshock + estim_params_.ncx;
|
||||
nshock = nshock + estim_params_.ncn;
|
||||
|
||||
c0=corrcoef(x);
|
||||
c00=tril(c0,-1);
|
||||
fig_nam_=[fname_,'_',fnam,'_corr_'];
|
||||
|
||||
ifig=0;
|
||||
j2=0;
|
||||
if ishock==0
|
||||
npar=estim_params_.np;
|
||||
else
|
||||
npar=estim_params_.np+nshock;
|
||||
end
|
||||
for j=1:npar,
|
||||
i2=find(abs(c00(:,j))>alpha2);
|
||||
if length(i2)>0,
|
||||
for jx=1:length(i2),
|
||||
tval = abs(c00(i2(jx),j)*sqrt( (nsam-2)/(1-c00(i2(jx),j)^2) ));
|
||||
tcr = tcrit(nsam-2,pvalue);
|
||||
if tval>tcr,
|
||||
j2=j2+1;
|
||||
if mod(j2,12)==1,
|
||||
ifig=ifig+1;
|
||||
hh=figure('name',['Correlations in the ',fnam,' sample ', num2str(ifig)]);
|
||||
end
|
||||
subplot(3,4,j2-(ifig-1)*12)
|
||||
% bar(c0(i2,j)),
|
||||
% set(gca,'xticklabel',bayestopt_.name(i2)),
|
||||
% set(gca,'xtick',[1:length(i2)])
|
||||
%plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
|
||||
%hold on,
|
||||
plot(x(:,j),x(:,i2(jx)),'.')
|
||||
% xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
|
||||
% ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
|
||||
if ishock,
|
||||
xlabel(bayestopt_.name{j},'interpreter','none'),
|
||||
ylabel(bayestopt_.name{i2(jx)},'interpreter','none'),
|
||||
else
|
||||
xlabel(bayestopt_.name{j+nshock},'interpreter','none'),
|
||||
ylabel(bayestopt_.name{i2(jx)+nshock},'interpreter','none'),
|
||||
end
|
||||
title(['cc = ',num2str(c0(i2(jx),j))])
|
||||
if (mod(j2,12)==0) && j2>0,
|
||||
eval(['print -depsc2 ' dirname '/' fig_nam_ int2str(ifig) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION'),
|
||||
saveas(hh,[dirname,'/',fig_nam_,int2str(ifig)])
|
||||
eval(['print -dpdf ' dirname '/' fig_nam_ int2str(ifig)]);
|
||||
end
|
||||
if options_.nograph, close(hh), end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
if (j==(npar)) && j2>0,
|
||||
eval(['print -depsc2 ' dirname '/' fig_nam_ int2str(ifig) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION'),
|
||||
saveas(gcf,[dirname,'/',fig_nam_,int2str(ifig)])
|
||||
eval(['print -dpdf ' dirname '/' fig_nam_ int2str(ifig)]);
|
||||
end
|
||||
if options_.nograph, close(gcf), end
|
||||
end
|
||||
|
||||
end
|
||||
if ifig==0,
|
||||
disp(['No correlation term >', num2str(alpha2),' found for ',fnam])
|
||||
end
|
||||
%close all
|
||||
function stab_map_2(x,alpha2, pvalue, fnam, dirname)
|
||||
% function stab_map_2(x, alpha2, pvalue, fnam, dirname)
|
||||
%
|
||||
% 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 bayestopt_ estim_params_ dr_ options_ ys_ fname_
|
||||
global bayestopt_ estim_params_ options_ oo_ M_
|
||||
|
||||
npar=size(x,2);
|
||||
nsam=size(x,1);
|
||||
ishock= npar>estim_params_.np;
|
||||
if nargin<4,
|
||||
fnam='';
|
||||
end
|
||||
if nargin<5,
|
||||
dirname='';
|
||||
end
|
||||
|
||||
ys_ = oo_.dr.ys;
|
||||
dr_ = oo_.dr;
|
||||
fname_ = M_.fname;
|
||||
nshock = estim_params_.nvx;
|
||||
nshock = nshock + estim_params_.nvn;
|
||||
nshock = nshock + estim_params_.ncx;
|
||||
nshock = nshock + estim_params_.ncn;
|
||||
|
||||
c0=corrcoef(x);
|
||||
c00=tril(c0,-1);
|
||||
fig_nam_=[fname_,'_',fnam,'_corr_'];
|
||||
|
||||
ifig=0;
|
||||
j2=0;
|
||||
if ishock==0
|
||||
npar=estim_params_.np;
|
||||
else
|
||||
npar=estim_params_.np+nshock;
|
||||
end
|
||||
for j=1:npar,
|
||||
i2=find(abs(c00(:,j))>alpha2);
|
||||
if length(i2)>0,
|
||||
for jx=1:length(i2),
|
||||
tval = abs(c00(i2(jx),j)*sqrt( (nsam-2)/(1-c00(i2(jx),j)^2) ));
|
||||
tcr = tcrit(nsam-2,pvalue);
|
||||
if tval>tcr,
|
||||
j2=j2+1;
|
||||
if mod(j2,12)==1,
|
||||
ifig=ifig+1;
|
||||
hh=figure('name',['Correlations in the ',fnam,' sample ', num2str(ifig)]);
|
||||
end
|
||||
subplot(3,4,j2-(ifig-1)*12)
|
||||
% bar(c0(i2,j)),
|
||||
% set(gca,'xticklabel',bayestopt_.name(i2)),
|
||||
% set(gca,'xtick',[1:length(i2)])
|
||||
%plot(stock_par(ixx(nfilt+1:end,i),j),stock_par(ixx(nfilt+1:end,i),i2(jx)),'.k')
|
||||
%hold on,
|
||||
plot(x(:,j),x(:,i2(jx)),'.')
|
||||
% xlabel(deblank(estim_params_.param_names(j,:)),'interpreter','none'),
|
||||
% ylabel(deblank(estim_params_.param_names(i2(jx),:)),'interpreter','none'),
|
||||
if ishock,
|
||||
xlabel(bayestopt_.name{j},'interpreter','none'),
|
||||
ylabel(bayestopt_.name{i2(jx)},'interpreter','none'),
|
||||
else
|
||||
xlabel(bayestopt_.name{j+nshock},'interpreter','none'),
|
||||
ylabel(bayestopt_.name{i2(jx)+nshock},'interpreter','none'),
|
||||
end
|
||||
title(['cc = ',num2str(c0(i2(jx),j))])
|
||||
if (mod(j2,12)==0) && j2>0,
|
||||
eval(['print -depsc2 ' dirname '/' fig_nam_ int2str(ifig) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION'),
|
||||
saveas(hh,[dirname,'/',fig_nam_,int2str(ifig)])
|
||||
eval(['print -dpdf ' dirname '/' fig_nam_ int2str(ifig)]);
|
||||
end
|
||||
if options_.nograph, close(hh), end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
if (j==(npar)) && j2>0,
|
||||
eval(['print -depsc2 ' dirname '/' fig_nam_ int2str(ifig) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION'),
|
||||
saveas(gcf,[dirname,'/',fig_nam_,int2str(ifig)])
|
||||
eval(['print -dpdf ' dirname '/' fig_nam_ int2str(ifig)]);
|
||||
end
|
||||
if options_.nograph, close(gcf), end
|
||||
end
|
||||
|
||||
end
|
||||
if ifig==0,
|
||||
disp(['No correlation term >', num2str(alpha2),' found for ',fnam])
|
||||
end
|
||||
%close all
|
||||
|
|
|
@ -1,25 +1,25 @@
|
|||
function [y, meany, stdy] = stand_(x)
|
||||
% STAND_ Standardise a matrix by columns
|
||||
%
|
||||
% [x,my,sy]=stand(y)
|
||||
%
|
||||
% y: 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
|
||||
%
|
||||
% Copyright (c) 2006 by JRC, European Commission, United Kingdom
|
||||
% Author : Marco Ratto
|
||||
|
||||
|
||||
if nargin==0,
|
||||
return;
|
||||
end
|
||||
|
||||
for j=1:size(x,2);
|
||||
meany(j)=mean(x(find(~isnan(x(:,j))),j));
|
||||
stdy(j)=std(x(find(~isnan(x(:,j))),j));
|
||||
y(:,j)=(x(:,j)-meany(j))./stdy(j);
|
||||
end
|
||||
function [y, meany, stdy] = stand_(x)
|
||||
% STAND_ Standardise a matrix by columns
|
||||
%
|
||||
% [x,my,sy]=stand(y)
|
||||
%
|
||||
% y: 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
|
||||
%
|
||||
% Copyright (c) 2006 by JRC, European Commission, United Kingdom
|
||||
% Author : Marco Ratto
|
||||
|
||||
|
||||
if nargin==0,
|
||||
return;
|
||||
end
|
||||
|
||||
for j=1:size(x,2);
|
||||
meany(j)=mean(x(find(~isnan(x(:,j))),j));
|
||||
stdy(j)=std(x(find(~isnan(x(:,j))),j));
|
||||
y(:,j)=(x(:,j)-meany(j))./stdy(j);
|
||||
end
|
||||
% end of m-file
|
|
@ -1,149 +1,149 @@
|
|||
function t_crit = tcrit(n,pval0)
|
||||
|
||||
% function t_crit = tcrit(n,pval0)
|
||||
%
|
||||
% given the p-value pval0, the function givese the
|
||||
% critical value t_crit of the t-distribution with n degress of freedom
|
||||
%
|
||||
% Part of the Sensitivity Analysis Toolbox for DYNARE
|
||||
%
|
||||
% Written by Marco Ratto
|
||||
% Joint Research Centre, The European Commission,
|
||||
% (http://eemc.jrc.ec.europa.eu/),
|
||||
% marco.ratto@jrc.it
|
||||
%
|
||||
|
||||
% Copyright (C) 2011-2011 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
if nargin==1 || isempty(pval0),
|
||||
pval0=0.05;
|
||||
end
|
||||
pval = [ 0.10 0.05 0.025 0.01 0.005 0.001];
|
||||
pval0=max(pval0,min(pval));
|
||||
ncol=min(find(pval<=pval0))+1;
|
||||
|
||||
t_crit=[
|
||||
1 3.078 6.314 12.706 31.821 63.657 318.313
|
||||
2 1.886 2.920 4.303 6.965 9.925 22.327
|
||||
3 1.638 2.353 3.182 4.541 5.841 10.215
|
||||
4 1.533 2.132 2.776 3.747 4.604 7.173
|
||||
5 1.476 2.015 2.571 3.365 4.032 5.893
|
||||
6 1.440 1.943 2.447 3.143 3.707 5.208
|
||||
7 1.415 1.895 2.365 2.998 3.499 4.782
|
||||
8 1.397 1.860 2.306 2.896 3.355 4.499
|
||||
9 1.383 1.833 2.262 2.821 3.250 4.296
|
||||
10 1.372 1.812 2.228 2.764 3.169 4.143
|
||||
11 1.363 1.796 2.201 2.718 3.106 4.024
|
||||
12 1.356 1.782 2.179 2.681 3.055 3.929
|
||||
13 1.350 1.771 2.160 2.650 3.012 3.852
|
||||
14 1.345 1.761 2.145 2.624 2.977 3.787
|
||||
15 1.341 1.753 2.131 2.602 2.947 3.733
|
||||
16 1.337 1.746 2.120 2.583 2.921 3.686
|
||||
17 1.333 1.740 2.110 2.567 2.898 3.646
|
||||
18 1.330 1.734 2.101 2.552 2.878 3.610
|
||||
19 1.328 1.729 2.093 2.539 2.861 3.579
|
||||
20 1.325 1.725 2.086 2.528 2.845 3.552
|
||||
21 1.323 1.721 2.080 2.518 2.831 3.527
|
||||
22 1.321 1.717 2.074 2.508 2.819 3.505
|
||||
23 1.319 1.714 2.069 2.500 2.807 3.485
|
||||
24 1.318 1.711 2.064 2.492 2.797 3.467
|
||||
25 1.316 1.708 2.060 2.485 2.787 3.450
|
||||
26 1.315 1.706 2.056 2.479 2.779 3.435
|
||||
27 1.314 1.703 2.052 2.473 2.771 3.421
|
||||
28 1.313 1.701 2.048 2.467 2.763 3.408
|
||||
29 1.311 1.699 2.045 2.462 2.756 3.396
|
||||
30 1.310 1.697 2.042 2.457 2.750 3.385
|
||||
31 1.309 1.696 2.040 2.453 2.744 3.375
|
||||
32 1.309 1.694 2.037 2.449 2.738 3.365
|
||||
33 1.308 1.692 2.035 2.445 2.733 3.356
|
||||
34 1.307 1.691 2.032 2.441 2.728 3.348
|
||||
35 1.306 1.690 2.030 2.438 2.724 3.340
|
||||
36 1.306 1.688 2.028 2.434 2.719 3.333
|
||||
37 1.305 1.687 2.026 2.431 2.715 3.326
|
||||
38 1.304 1.686 2.024 2.429 2.712 3.319
|
||||
39 1.304 1.685 2.023 2.426 2.708 3.313
|
||||
40 1.303 1.684 2.021 2.423 2.704 3.307
|
||||
41 1.303 1.683 2.020 2.421 2.701 3.301
|
||||
42 1.302 1.682 2.018 2.418 2.698 3.296
|
||||
43 1.302 1.681 2.017 2.416 2.695 3.291
|
||||
44 1.301 1.680 2.015 2.414 2.692 3.286
|
||||
45 1.301 1.679 2.014 2.412 2.690 3.281
|
||||
46 1.300 1.679 2.013 2.410 2.687 3.277
|
||||
47 1.300 1.678 2.012 2.408 2.685 3.273
|
||||
48 1.299 1.677 2.011 2.407 2.682 3.269
|
||||
49 1.299 1.677 2.010 2.405 2.680 3.265
|
||||
50 1.299 1.676 2.009 2.403 2.678 3.261
|
||||
51 1.298 1.675 2.008 2.402 2.676 3.258
|
||||
52 1.298 1.675 2.007 2.400 2.674 3.255
|
||||
53 1.298 1.674 2.006 2.399 2.672 3.251
|
||||
54 1.297 1.674 2.005 2.397 2.670 3.248
|
||||
55 1.297 1.673 2.004 2.396 2.668 3.245
|
||||
56 1.297 1.673 2.003 2.395 2.667 3.242
|
||||
57 1.297 1.672 2.002 2.394 2.665 3.239
|
||||
58 1.296 1.672 2.002 2.392 2.663 3.237
|
||||
59 1.296 1.671 2.001 2.391 2.662 3.234
|
||||
60 1.296 1.671 2.000 2.390 2.660 3.232
|
||||
61 1.296 1.670 2.000 2.389 2.659 3.229
|
||||
62 1.295 1.670 1.999 2.388 2.657 3.227
|
||||
63 1.295 1.669 1.998 2.387 2.656 3.225
|
||||
64 1.295 1.669 1.998 2.386 2.655 3.223
|
||||
65 1.295 1.669 1.997 2.385 2.654 3.220
|
||||
66 1.295 1.668 1.997 2.384 2.652 3.218
|
||||
67 1.294 1.668 1.996 2.383 2.651 3.216
|
||||
68 1.294 1.668 1.995 2.382 2.650 3.214
|
||||
69 1.294 1.667 1.995 2.382 2.649 3.213
|
||||
70 1.294 1.667 1.994 2.381 2.648 3.211
|
||||
71 1.294 1.667 1.994 2.380 2.647 3.209
|
||||
72 1.293 1.666 1.993 2.379 2.646 3.207
|
||||
73 1.293 1.666 1.993 2.379 2.645 3.206
|
||||
74 1.293 1.666 1.993 2.378 2.644 3.204
|
||||
75 1.293 1.665 1.992 2.377 2.643 3.202
|
||||
76 1.293 1.665 1.992 2.376 2.642 3.201
|
||||
77 1.293 1.665 1.991 2.376 2.641 3.199
|
||||
78 1.292 1.665 1.991 2.375 2.640 3.198
|
||||
79 1.292 1.664 1.990 2.374 2.640 3.197
|
||||
80 1.292 1.664 1.990 2.374 2.639 3.195
|
||||
81 1.292 1.664 1.990 2.373 2.638 3.194
|
||||
82 1.292 1.664 1.989 2.373 2.637 3.193
|
||||
83 1.292 1.663 1.989 2.372 2.636 3.191
|
||||
84 1.292 1.663 1.989 2.372 2.636 3.190
|
||||
85 1.292 1.663 1.988 2.371 2.635 3.189
|
||||
86 1.291 1.663 1.988 2.370 2.634 3.188
|
||||
87 1.291 1.663 1.988 2.370 2.634 3.187
|
||||
88 1.291 1.662 1.987 2.369 2.633 3.185
|
||||
89 1.291 1.662 1.987 2.369 2.632 3.184
|
||||
90 1.291 1.662 1.987 2.368 2.632 3.183
|
||||
91 1.291 1.662 1.986 2.368 2.631 3.182
|
||||
92 1.291 1.662 1.986 2.368 2.630 3.181
|
||||
93 1.291 1.661 1.986 2.367 2.630 3.180
|
||||
94 1.291 1.661 1.986 2.367 2.629 3.179
|
||||
95 1.291 1.661 1.985 2.366 2.629 3.178
|
||||
96 1.290 1.661 1.985 2.366 2.628 3.177
|
||||
97 1.290 1.661 1.985 2.365 2.627 3.176
|
||||
98 1.290 1.661 1.984 2.365 2.627 3.175
|
||||
99 1.290 1.660 1.984 2.365 2.626 3.175
|
||||
100 1.290 1.660 1.984 2.364 2.626 3.174
|
||||
inf 1.282 1.645 1.960 2.326 2.576 3.090
|
||||
];
|
||||
|
||||
if n<=100,
|
||||
t_crit=t_crit(n,ncol);
|
||||
else
|
||||
t_crit=t_crit(end,ncol);
|
||||
function t_crit = tcrit(n,pval0)
|
||||
|
||||
% function t_crit = tcrit(n,pval0)
|
||||
%
|
||||
% given the p-value pval0, the function givese the
|
||||
% critical value t_crit of the t-distribution with n degress of freedom
|
||||
%
|
||||
% Part of the Sensitivity Analysis Toolbox for DYNARE
|
||||
%
|
||||
% Written by Marco Ratto
|
||||
% Joint Research Centre, The European Commission,
|
||||
% (http://eemc.jrc.ec.europa.eu/),
|
||||
% marco.ratto@jrc.it
|
||||
%
|
||||
|
||||
% Copyright (C) 2011-2011 Dynare Team
|
||||
%
|
||||
% This file is part of Dynare.
|
||||
%
|
||||
% Dynare is free software: you can redistribute it and/or modify
|
||||
% it under the terms of the GNU General Public License as published by
|
||||
% the Free Software Foundation, either version 3 of the License, or
|
||||
% (at your option) any later version.
|
||||
%
|
||||
% Dynare is distributed in the hope that it will be useful,
|
||||
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
% GNU General Public License for more details.
|
||||
%
|
||||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
if nargin==1 || isempty(pval0),
|
||||
pval0=0.05;
|
||||
end
|
||||
pval = [ 0.10 0.05 0.025 0.01 0.005 0.001];
|
||||
pval0=max(pval0,min(pval));
|
||||
ncol=min(find(pval<=pval0))+1;
|
||||
|
||||
t_crit=[
|
||||
1 3.078 6.314 12.706 31.821 63.657 318.313
|
||||
2 1.886 2.920 4.303 6.965 9.925 22.327
|
||||
3 1.638 2.353 3.182 4.541 5.841 10.215
|
||||
4 1.533 2.132 2.776 3.747 4.604 7.173
|
||||
5 1.476 2.015 2.571 3.365 4.032 5.893
|
||||
6 1.440 1.943 2.447 3.143 3.707 5.208
|
||||
7 1.415 1.895 2.365 2.998 3.499 4.782
|
||||
8 1.397 1.860 2.306 2.896 3.355 4.499
|
||||
9 1.383 1.833 2.262 2.821 3.250 4.296
|
||||
10 1.372 1.812 2.228 2.764 3.169 4.143
|
||||
11 1.363 1.796 2.201 2.718 3.106 4.024
|
||||
12 1.356 1.782 2.179 2.681 3.055 3.929
|
||||
13 1.350 1.771 2.160 2.650 3.012 3.852
|
||||
14 1.345 1.761 2.145 2.624 2.977 3.787
|
||||
15 1.341 1.753 2.131 2.602 2.947 3.733
|
||||
16 1.337 1.746 2.120 2.583 2.921 3.686
|
||||
17 1.333 1.740 2.110 2.567 2.898 3.646
|
||||
18 1.330 1.734 2.101 2.552 2.878 3.610
|
||||
19 1.328 1.729 2.093 2.539 2.861 3.579
|
||||
20 1.325 1.725 2.086 2.528 2.845 3.552
|
||||
21 1.323 1.721 2.080 2.518 2.831 3.527
|
||||
22 1.321 1.717 2.074 2.508 2.819 3.505
|
||||
23 1.319 1.714 2.069 2.500 2.807 3.485
|
||||
24 1.318 1.711 2.064 2.492 2.797 3.467
|
||||
25 1.316 1.708 2.060 2.485 2.787 3.450
|
||||
26 1.315 1.706 2.056 2.479 2.779 3.435
|
||||
27 1.314 1.703 2.052 2.473 2.771 3.421
|
||||
28 1.313 1.701 2.048 2.467 2.763 3.408
|
||||
29 1.311 1.699 2.045 2.462 2.756 3.396
|
||||
30 1.310 1.697 2.042 2.457 2.750 3.385
|
||||
31 1.309 1.696 2.040 2.453 2.744 3.375
|
||||
32 1.309 1.694 2.037 2.449 2.738 3.365
|
||||
33 1.308 1.692 2.035 2.445 2.733 3.356
|
||||
34 1.307 1.691 2.032 2.441 2.728 3.348
|
||||
35 1.306 1.690 2.030 2.438 2.724 3.340
|
||||
36 1.306 1.688 2.028 2.434 2.719 3.333
|
||||
37 1.305 1.687 2.026 2.431 2.715 3.326
|
||||
38 1.304 1.686 2.024 2.429 2.712 3.319
|
||||
39 1.304 1.685 2.023 2.426 2.708 3.313
|
||||
40 1.303 1.684 2.021 2.423 2.704 3.307
|
||||
41 1.303 1.683 2.020 2.421 2.701 3.301
|
||||
42 1.302 1.682 2.018 2.418 2.698 3.296
|
||||
43 1.302 1.681 2.017 2.416 2.695 3.291
|
||||
44 1.301 1.680 2.015 2.414 2.692 3.286
|
||||
45 1.301 1.679 2.014 2.412 2.690 3.281
|
||||
46 1.300 1.679 2.013 2.410 2.687 3.277
|
||||
47 1.300 1.678 2.012 2.408 2.685 3.273
|
||||
48 1.299 1.677 2.011 2.407 2.682 3.269
|
||||
49 1.299 1.677 2.010 2.405 2.680 3.265
|
||||
50 1.299 1.676 2.009 2.403 2.678 3.261
|
||||
51 1.298 1.675 2.008 2.402 2.676 3.258
|
||||
52 1.298 1.675 2.007 2.400 2.674 3.255
|
||||
53 1.298 1.674 2.006 2.399 2.672 3.251
|
||||
54 1.297 1.674 2.005 2.397 2.670 3.248
|
||||
55 1.297 1.673 2.004 2.396 2.668 3.245
|
||||
56 1.297 1.673 2.003 2.395 2.667 3.242
|
||||
57 1.297 1.672 2.002 2.394 2.665 3.239
|
||||
58 1.296 1.672 2.002 2.392 2.663 3.237
|
||||
59 1.296 1.671 2.001 2.391 2.662 3.234
|
||||
60 1.296 1.671 2.000 2.390 2.660 3.232
|
||||
61 1.296 1.670 2.000 2.389 2.659 3.229
|
||||
62 1.295 1.670 1.999 2.388 2.657 3.227
|
||||
63 1.295 1.669 1.998 2.387 2.656 3.225
|
||||
64 1.295 1.669 1.998 2.386 2.655 3.223
|
||||
65 1.295 1.669 1.997 2.385 2.654 3.220
|
||||
66 1.295 1.668 1.997 2.384 2.652 3.218
|
||||
67 1.294 1.668 1.996 2.383 2.651 3.216
|
||||
68 1.294 1.668 1.995 2.382 2.650 3.214
|
||||
69 1.294 1.667 1.995 2.382 2.649 3.213
|
||||
70 1.294 1.667 1.994 2.381 2.648 3.211
|
||||
71 1.294 1.667 1.994 2.380 2.647 3.209
|
||||
72 1.293 1.666 1.993 2.379 2.646 3.207
|
||||
73 1.293 1.666 1.993 2.379 2.645 3.206
|
||||
74 1.293 1.666 1.993 2.378 2.644 3.204
|
||||
75 1.293 1.665 1.992 2.377 2.643 3.202
|
||||
76 1.293 1.665 1.992 2.376 2.642 3.201
|
||||
77 1.293 1.665 1.991 2.376 2.641 3.199
|
||||
78 1.292 1.665 1.991 2.375 2.640 3.198
|
||||
79 1.292 1.664 1.990 2.374 2.640 3.197
|
||||
80 1.292 1.664 1.990 2.374 2.639 3.195
|
||||
81 1.292 1.664 1.990 2.373 2.638 3.194
|
||||
82 1.292 1.664 1.989 2.373 2.637 3.193
|
||||
83 1.292 1.663 1.989 2.372 2.636 3.191
|
||||
84 1.292 1.663 1.989 2.372 2.636 3.190
|
||||
85 1.292 1.663 1.988 2.371 2.635 3.189
|
||||
86 1.291 1.663 1.988 2.370 2.634 3.188
|
||||
87 1.291 1.663 1.988 2.370 2.634 3.187
|
||||
88 1.291 1.662 1.987 2.369 2.633 3.185
|
||||
89 1.291 1.662 1.987 2.369 2.632 3.184
|
||||
90 1.291 1.662 1.987 2.368 2.632 3.183
|
||||
91 1.291 1.662 1.986 2.368 2.631 3.182
|
||||
92 1.291 1.662 1.986 2.368 2.630 3.181
|
||||
93 1.291 1.661 1.986 2.367 2.630 3.180
|
||||
94 1.291 1.661 1.986 2.367 2.629 3.179
|
||||
95 1.291 1.661 1.985 2.366 2.629 3.178
|
||||
96 1.290 1.661 1.985 2.366 2.628 3.177
|
||||
97 1.290 1.661 1.985 2.365 2.627 3.176
|
||||
98 1.290 1.661 1.984 2.365 2.627 3.175
|
||||
99 1.290 1.660 1.984 2.365 2.626 3.175
|
||||
100 1.290 1.660 1.984 2.364 2.626 3.174
|
||||
inf 1.282 1.645 1.960 2.326 2.576 3.090
|
||||
];
|
||||
|
||||
if n<=100,
|
||||
t_crit=t_crit(n,ncol);
|
||||
else
|
||||
t_crit=t_crit(end,ncol);
|
||||
end
|
|
@ -1,51 +1,51 @@
|
|||
function [yt, j0, ir, ic]=teff(T,Nsam,istable)
|
||||
%
|
||||
% [yt, j0, ir, ic]=teff(T,Nsam,istable)
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
ndim = (length(size(T)));
|
||||
if ndim==3,
|
||||
if nargin==1,
|
||||
Nsam=size(T,3);
|
||||
istable = [1:Nsam]';
|
||||
end
|
||||
tmax=max(T,[],3);
|
||||
tmin=min(T,[],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);
|
||||
ir=(find( (tmax-tmin)>1.e-8));
|
||||
j0 = length(ir);
|
||||
yt=NaN(Nsam, j0);
|
||||
yt(istable,:)=T(ir,:)';
|
||||
|
||||
|
||||
end
|
||||
%clear y0 y1;
|
||||
function [yt, j0, ir, ic]=teff(T,Nsam,istable)
|
||||
%
|
||||
% [yt, j0, ir, ic]=teff(T,Nsam,istable)
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
ndim = (length(size(T)));
|
||||
if ndim==3,
|
||||
if nargin==1,
|
||||
Nsam=size(T,3);
|
||||
istable = [1:Nsam]';
|
||||
end
|
||||
tmax=max(T,[],3);
|
||||
tmin=min(T,[],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);
|
||||
ir=(find( (tmax-tmin)>1.e-8));
|
||||
j0 = length(ir);
|
||||
yt=NaN(Nsam, j0);
|
||||
yt(istable,:)=T(ir,:)';
|
||||
|
||||
|
||||
end
|
||||
%clear y0 y1;
|
||||
|
|
|
@ -1,44 +1,44 @@
|
|||
function tau = thet2tau(params, indx, indexo, flagmoments,mf,nlags,useautocorr)
|
||||
global M_ oo_ options_
|
||||
|
||||
if nargin==1,
|
||||
indx = [1:M_.param_nbr];
|
||||
indexo = [];
|
||||
end
|
||||
|
||||
if nargin<4,
|
||||
flagmoments=0;
|
||||
end
|
||||
if nargin<7 | isempty(useautocorr),
|
||||
useautocorr=0;
|
||||
end
|
||||
|
||||
M_.params(indx) = params(length(indexo)+1:end);
|
||||
if ~isempty(indexo)
|
||||
M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
|
||||
end
|
||||
[A,B,plouf,plouf,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
|
||||
if flagmoments==0,
|
||||
tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];
|
||||
else
|
||||
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
|
||||
k = find(abs(GAM) < 1e-12);
|
||||
GAM(k) = 0;
|
||||
if useautocorr,
|
||||
sy = sqrt(diag(GAM));
|
||||
sy = sy*sy';
|
||||
sy0 = sy-diag(diag(sy))+eye(length(sy));
|
||||
dum = GAM./sy0;
|
||||
tau = dyn_vech(dum(mf,mf));
|
||||
else
|
||||
tau = dyn_vech(GAM(mf,mf));
|
||||
end
|
||||
for ii = 1:nlags
|
||||
dum = A^(ii)*GAM;
|
||||
if useautocorr,
|
||||
dum = dum./sy;
|
||||
end
|
||||
tau = [tau;vec(dum(mf,mf))];
|
||||
end
|
||||
tau = [ oo_.dr.ys(oo_.dr.order_var(mf)); tau];
|
||||
function tau = thet2tau(params, indx, indexo, flagmoments,mf,nlags,useautocorr)
|
||||
global M_ oo_ options_
|
||||
|
||||
if nargin==1,
|
||||
indx = [1:M_.param_nbr];
|
||||
indexo = [];
|
||||
end
|
||||
|
||||
if nargin<4,
|
||||
flagmoments=0;
|
||||
end
|
||||
if nargin<7 | isempty(useautocorr),
|
||||
useautocorr=0;
|
||||
end
|
||||
|
||||
M_.params(indx) = params(length(indexo)+1:end);
|
||||
if ~isempty(indexo)
|
||||
M_.Sigma_e(indexo,indexo) = diag(params(1:length(indexo)).^2);
|
||||
end
|
||||
[A,B,plouf,plouf,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
|
||||
if flagmoments==0,
|
||||
tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];
|
||||
else
|
||||
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
|
||||
k = find(abs(GAM) < 1e-12);
|
||||
GAM(k) = 0;
|
||||
if useautocorr,
|
||||
sy = sqrt(diag(GAM));
|
||||
sy = sy*sy';
|
||||
sy0 = sy-diag(diag(sy))+eye(length(sy));
|
||||
dum = GAM./sy0;
|
||||
tau = dyn_vech(dum(mf,mf));
|
||||
else
|
||||
tau = dyn_vech(GAM(mf,mf));
|
||||
end
|
||||
for ii = 1:nlags
|
||||
dum = A^(ii)*GAM;
|
||||
if useautocorr,
|
||||
dum = dum./sy;
|
||||
end
|
||||
tau = [tau;vec(dum(mf,mf))];
|
||||
end
|
||||
tau = [ oo_.dr.ys(oo_.dr.order_var(mf)); tau];
|
||||
end
|
|
@ -1,25 +1,25 @@
|
|||
function yr = trank(y);
|
||||
% yr = trank(y);
|
||||
% yr is the rank transformation of y
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
[nr, nc] = size(y);
|
||||
for j=1:nc,
|
||||
[dum, is]=sort(y(:,j));
|
||||
yr(is,j)=[1:nr]'./nr;
|
||||
end
|
||||
function yr = trank(y);
|
||||
% yr = trank(y);
|
||||
% yr is the rank transformation of y
|
||||
%
|
||||
% 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.
|
||||
%
|
||||
|
||||
[nr, nc] = size(y);
|
||||
for j=1:nc,
|
||||
[dum, is]=sort(y(:,j));
|
||||
yr(is,j)=[1:nr]'./nr;
|
||||
end
|
||||
|
|
|
@ -1,100 +1,100 @@
|
|||
data = [0.928467646476 11.8716889412 20 0.418037507392 0.227382377518 ...
|
||||
-0.705994063083 11.7522582094 21.25 1.09254424511 -1.29488274994 ...
|
||||
-0.511895351926 9.68144025625 17.25 -1.66150408407 0.331508393098 ...
|
||||
-0.990955971267 10.0890781236 17 1.43016275252 -2.43589670141 ...
|
||||
-0.981233061806 12.1094840679 18.25 2.91293288733 -0.790246576864 ...
|
||||
-0.882182844512 8.54559460406 15 0.419579139481 0.358729719566 ...
|
||||
-0.930893002836 6.19238374422 12.5 -1.48847457959 0.739779938797 ...
|
||||
1.53158206947 2.76544271886 11.5 -0.336216769682 0.455559918769 ...
|
||||
2.2659052834 5.47418162513 11 0.306436789767 -0.0707985731221 ...
|
||||
1.05419803797 6.35698426189 11 0.140700250477 0.620401487202 ...
|
||||
1.20161076793 3.4253301593 11 0.461296492351 0.14354323987 ...
|
||||
1.73934077971 4.70926070322 11.5 1.35798282982 0.38564694435 ...
|
||||
1.71735262584 3.54232079749 12.5 2.9097529155 -0.804308583301 ...
|
||||
0.426343657844 3.32719108897 13 1.64214862652 -1.18214664701 ...
|
||||
1.67751812324 2.93444727338 11.25 0.344434910651 -1.6529373719 ...
|
||||
1.37013301099 4.72303361923 11.75 2.61511526582 0.327684243041 ...
|
||||
0.281231073781 4.4893853071 10.5 1.17043449257 1.12855106649 ...
|
||||
1.53638992834 3.7325309699 10.25 -0.683947046728 0.11943538737 ...
|
||||
1.68081431462 3.34729969129 10 1.41159342106 -1.59065680853 ...
|
||||
-0.343321601133 5.05563513564 12 1.75117366498 -2.40127764642 ...
|
||||
0.873415608666 3.2779996255 10.25 -1.39895866711 0.0971444398216 ...
|
||||
0.26399696544 4.78229419828 9.75 0.0914692438124 0.299310457612 ...
|
||||
-0.562233624818 3.88598638237 9.75 -0.0505384765105 0.332826708151 ...
|
||||
2.15161914936 3.84859710132 8.75 -3.44811080489 0.789138678784 ...
|
||||
1.2345093726 5.62225030942 9.5 -0.366945407434 2.32974981198 ...
|
||||
1.62554967459 4.24667132831 10 -0.800958371402 0.0293183770935 ...
|
||||
1.33035402527 2.75248979249 9.75 -0.855723113225 0.852493939813 ...
|
||||
1.52078814077 3.53415985826 9.75 -3.37963469203 -1.05133958119 ...
|
||||
1.16704983697 4.92754079464 10.75 -3.0142303324 0.459907431978 ...
|
||||
0.277213572101 4.55532133037 11.75 -0.851995599415 2.03242034852 ...
|
||||
0.842215068977 3.11164509647 12.25 -1.08290421696 0.014323281961 ...
|
||||
1.05325028606 4.92882647578 13.5 -1.1953883867 0.706764750654 ...
|
||||
0.453051253568 6.82998950103 13.5 0.111803656462 0.088462593153 ...
|
||||
0.199885995525 5.82643354662 13.5 -0.920501518421 -0.26504958666 ...
|
||||
0.137907999624 2.66076369132 13.5 -1.17122929812 -0.995642430514 ...
|
||||
0.721949686709 5.70497876823 14.25 1.19378169018 -1.10644839651 ...
|
||||
-0.418465249225 3.75861110232 14.75 -1.03131674824 0.188507675831 ...
|
||||
-0.644028342116 4.15104788154 13.75 -1.48911756546 0.204560913792 ...
|
||||
-0.848213852668 5.65580324027 12.75 0.677011703877 -0.849628054542 ...
|
||||
-1.51954076928 11.4866911266 11.25 -0.446024680774 -0.456342350765 ...
|
||||
0.265275055215 2.85472749592 9.75 -0.598778202436 -0.907311640831 ...
|
||||
0.356162529063 2.29614015658 9.5 -0.46820788432 -1.22130883441 ...
|
||||
0.368308864363 -0.539083504685 8 -0.781333991956 0.374007246518 ...
|
||||
-0.145751412732 1.61507621789 8.25 3.68291932628 1.32438399845 ...
|
||||
0.285457283664 2.14334055993 7 1.42819405379 -0.00818660844123 ...
|
||||
0.372390129412 1.60000213334 6.25 0.626106424052 -0.10136772765 ...
|
||||
0.382720203063 1.72614243263 7.25 4.89631941021 -1.10060711916 ...
|
||||
0.737957515573 2.90430582851 6 -0.0422721010314 0.4178952497 ...
|
||||
0.649532581668 0.657135682543 6 0.692066153971 0.422299120276 ...
|
||||
0.627159201987 1.70352689913 5.75 2.62066711305 -1.29237304034 ...
|
||||
0.905441299817 1.95663197267 5.5 1.5949697565 -0.27115830703 ...
|
||||
1.49322577898 -2.08741765309 6.25 1.23027694802 0.418336889527 ...
|
||||
1.48750731567 -1.57274121871 8 3.01660550994 -0.893958254365 ...
|
||||
1.39783858087 2.22623066426 7 -0.80842319214 1.47625453886 ...
|
||||
0.89274836317 1.30378081742 8 -0.249485058661 0.159871204185 ...
|
||||
0.920652246088 4.1437741965 9.75 2.8204453623 0.178149239655 ...
|
||||
-0.00264276644799 3.07989972052 8.75 -2.56342461535 2.105998353 ...
|
||||
0.0198190461681 0.766283759256 8 -1.15838865989 1.56888883418 ...
|
||||
0.440050515311 0.127570085801 7.5 0.0400753569995 0.028914333532 ...
|
||||
0.129536637901 1.78174141526 6.75 0.959943962785 0.307781224401 ...
|
||||
0.398549827172 3.03606770667 6.5 -0.340209794742 0.100979469478 ...
|
||||
1.17174775425 0.629625188037 5.75 0.403003686814 0.902394579377 ...
|
||||
0.991163981251 2.50862910684 4.75 -1.44963996982 1.16150986945 ...
|
||||
0.967603566096 2.12003739013 4.75 0.610846030775 -0.889994896068 ...
|
||||
1.14689383604 1.24185011459 4.75 2.01098091308 -1.73846431001 ...
|
||||
1.32593824054 0.990713820685 4.75 -0.0955142989332 -0.0369257308362 ...
|
||||
0.861135002644 -0.24744943605 6 1.72793107135 -0.691506789639 ...
|
||||
1.26870850151 2.09844764887 6.5 1.50720217572 -1.31399187077 ...
|
||||
0.260364987715 1.10650139716 6.5 1.13659047496 0.0720441664643 ...
|
||||
1.09731242214 0.490796381346 7.25 4.59123894147 -2.14073070763 ...
|
||||
1.63792841781 0.612652594286 6.75 1.79604605035 -0.644363995357 ...
|
||||
1.48465576034 0.978295808687 6.75 -2.00753620902 1.39437534964 ...
|
||||
1.0987608663 4.25212569087 6.25 -2.58901196498 2.56054320803 ...
|
||||
1.42592178132 2.76984518311 6.25 0.888195752358 1.03114549274 ...
|
||||
1.52958239462 1.31795955491 6.5 -0.902907564082 -0.0952198893776 ...
|
||||
1.0170168994 2.14733589918 7 -1.3054866978 2.68803738466 ...
|
||||
0.723253652257 3.43552889347 7.5 1.8213700853 0.592593586195 ...
|
||||
1.24720806008 3.87383806577 7.5 0.0522300654168 0.988871238698 ...
|
||||
0.482531471239 2.67793287032 7.5 2.9693944293 -0.108591166081 ...
|
||||
0.154056100439 0.927269031704 6.75 0.119222057561 3.30489209451 ...
|
||||
0.0694865769274 6.65916526788 6.25 0.889014476084 -2.83976849035 ...
|
||||
-0.121267434867 0.341442615696 5.25 0.323053239216 -3.49289229012 ...
|
||||
0.726473690375 -3.5423730964 4 2.19149290449 -3.20855054004 ...
|
||||
1.39271709108 2.63121085718 3.75 0.88406577736 0.75622580197 ...
|
||||
1.07502077727 5.88578836799 4.25 -2.55088273352 2.89018116374 ...
|
||||
0.759049251607 4.24703604223 4.5 0.575687665685 -0.388292506167 ...
|
||||
];
|
||||
|
||||
data = reshape(data,5,86)';
|
||||
y_obs = data(:,1);
|
||||
pie_obs = data(:,2);
|
||||
R_obs = data(:,3);
|
||||
de = data(:,4);
|
||||
dq = data(:,5);
|
||||
|
||||
%Country: Canada
|
||||
%Sample Range: 1981:2 to 2002:3
|
||||
%Observations: 86
|
||||
%Variables: Real GDP Growth [%], Inflation [annualized %], Nom Rate [%],
|
||||
% Exchange Rate Change [%], Terms of Trade Change [%]
|
||||
data = [0.928467646476 11.8716889412 20 0.418037507392 0.227382377518 ...
|
||||
-0.705994063083 11.7522582094 21.25 1.09254424511 -1.29488274994 ...
|
||||
-0.511895351926 9.68144025625 17.25 -1.66150408407 0.331508393098 ...
|
||||
-0.990955971267 10.0890781236 17 1.43016275252 -2.43589670141 ...
|
||||
-0.981233061806 12.1094840679 18.25 2.91293288733 -0.790246576864 ...
|
||||
-0.882182844512 8.54559460406 15 0.419579139481 0.358729719566 ...
|
||||
-0.930893002836 6.19238374422 12.5 -1.48847457959 0.739779938797 ...
|
||||
1.53158206947 2.76544271886 11.5 -0.336216769682 0.455559918769 ...
|
||||
2.2659052834 5.47418162513 11 0.306436789767 -0.0707985731221 ...
|
||||
1.05419803797 6.35698426189 11 0.140700250477 0.620401487202 ...
|
||||
1.20161076793 3.4253301593 11 0.461296492351 0.14354323987 ...
|
||||
1.73934077971 4.70926070322 11.5 1.35798282982 0.38564694435 ...
|
||||
1.71735262584 3.54232079749 12.5 2.9097529155 -0.804308583301 ...
|
||||
0.426343657844 3.32719108897 13 1.64214862652 -1.18214664701 ...
|
||||
1.67751812324 2.93444727338 11.25 0.344434910651 -1.6529373719 ...
|
||||
1.37013301099 4.72303361923 11.75 2.61511526582 0.327684243041 ...
|
||||
0.281231073781 4.4893853071 10.5 1.17043449257 1.12855106649 ...
|
||||
1.53638992834 3.7325309699 10.25 -0.683947046728 0.11943538737 ...
|
||||
1.68081431462 3.34729969129 10 1.41159342106 -1.59065680853 ...
|
||||
-0.343321601133 5.05563513564 12 1.75117366498 -2.40127764642 ...
|
||||
0.873415608666 3.2779996255 10.25 -1.39895866711 0.0971444398216 ...
|
||||
0.26399696544 4.78229419828 9.75 0.0914692438124 0.299310457612 ...
|
||||
-0.562233624818 3.88598638237 9.75 -0.0505384765105 0.332826708151 ...
|
||||
2.15161914936 3.84859710132 8.75 -3.44811080489 0.789138678784 ...
|
||||
1.2345093726 5.62225030942 9.5 -0.366945407434 2.32974981198 ...
|
||||
1.62554967459 4.24667132831 10 -0.800958371402 0.0293183770935 ...
|
||||
1.33035402527 2.75248979249 9.75 -0.855723113225 0.852493939813 ...
|
||||
1.52078814077 3.53415985826 9.75 -3.37963469203 -1.05133958119 ...
|
||||
1.16704983697 4.92754079464 10.75 -3.0142303324 0.459907431978 ...
|
||||
0.277213572101 4.55532133037 11.75 -0.851995599415 2.03242034852 ...
|
||||
0.842215068977 3.11164509647 12.25 -1.08290421696 0.014323281961 ...
|
||||
1.05325028606 4.92882647578 13.5 -1.1953883867 0.706764750654 ...
|
||||
0.453051253568 6.82998950103 13.5 0.111803656462 0.088462593153 ...
|
||||
0.199885995525 5.82643354662 13.5 -0.920501518421 -0.26504958666 ...
|
||||
0.137907999624 2.66076369132 13.5 -1.17122929812 -0.995642430514 ...
|
||||
0.721949686709 5.70497876823 14.25 1.19378169018 -1.10644839651 ...
|
||||
-0.418465249225 3.75861110232 14.75 -1.03131674824 0.188507675831 ...
|
||||
-0.644028342116 4.15104788154 13.75 -1.48911756546 0.204560913792 ...
|
||||
-0.848213852668 5.65580324027 12.75 0.677011703877 -0.849628054542 ...
|
||||
-1.51954076928 11.4866911266 11.25 -0.446024680774 -0.456342350765 ...
|
||||
0.265275055215 2.85472749592 9.75 -0.598778202436 -0.907311640831 ...
|
||||
0.356162529063 2.29614015658 9.5 -0.46820788432 -1.22130883441 ...
|
||||
0.368308864363 -0.539083504685 8 -0.781333991956 0.374007246518 ...
|
||||
-0.145751412732 1.61507621789 8.25 3.68291932628 1.32438399845 ...
|
||||
0.285457283664 2.14334055993 7 1.42819405379 -0.00818660844123 ...
|
||||
0.372390129412 1.60000213334 6.25 0.626106424052 -0.10136772765 ...
|
||||
0.382720203063 1.72614243263 7.25 4.89631941021 -1.10060711916 ...
|
||||
0.737957515573 2.90430582851 6 -0.0422721010314 0.4178952497 ...
|
||||
0.649532581668 0.657135682543 6 0.692066153971 0.422299120276 ...
|
||||
0.627159201987 1.70352689913 5.75 2.62066711305 -1.29237304034 ...
|
||||
0.905441299817 1.95663197267 5.5 1.5949697565 -0.27115830703 ...
|
||||
1.49322577898 -2.08741765309 6.25 1.23027694802 0.418336889527 ...
|
||||
1.48750731567 -1.57274121871 8 3.01660550994 -0.893958254365 ...
|
||||
1.39783858087 2.22623066426 7 -0.80842319214 1.47625453886 ...
|
||||
0.89274836317 1.30378081742 8 -0.249485058661 0.159871204185 ...
|
||||
0.920652246088 4.1437741965 9.75 2.8204453623 0.178149239655 ...
|
||||
-0.00264276644799 3.07989972052 8.75 -2.56342461535 2.105998353 ...
|
||||
0.0198190461681 0.766283759256 8 -1.15838865989 1.56888883418 ...
|
||||
0.440050515311 0.127570085801 7.5 0.0400753569995 0.028914333532 ...
|
||||
0.129536637901 1.78174141526 6.75 0.959943962785 0.307781224401 ...
|
||||
0.398549827172 3.03606770667 6.5 -0.340209794742 0.100979469478 ...
|
||||
1.17174775425 0.629625188037 5.75 0.403003686814 0.902394579377 ...
|
||||
0.991163981251 2.50862910684 4.75 -1.44963996982 1.16150986945 ...
|
||||
0.967603566096 2.12003739013 4.75 0.610846030775 -0.889994896068 ...
|
||||
1.14689383604 1.24185011459 4.75 2.01098091308 -1.73846431001 ...
|
||||
1.32593824054 0.990713820685 4.75 -0.0955142989332 -0.0369257308362 ...
|
||||
0.861135002644 -0.24744943605 6 1.72793107135 -0.691506789639 ...
|
||||
1.26870850151 2.09844764887 6.5 1.50720217572 -1.31399187077 ...
|
||||
0.260364987715 1.10650139716 6.5 1.13659047496 0.0720441664643 ...
|
||||
1.09731242214 0.490796381346 7.25 4.59123894147 -2.14073070763 ...
|
||||
1.63792841781 0.612652594286 6.75 1.79604605035 -0.644363995357 ...
|
||||
1.48465576034 0.978295808687 6.75 -2.00753620902 1.39437534964 ...
|
||||
1.0987608663 4.25212569087 6.25 -2.58901196498 2.56054320803 ...
|
||||
1.42592178132 2.76984518311 6.25 0.888195752358 1.03114549274 ...
|
||||
1.52958239462 1.31795955491 6.5 -0.902907564082 -0.0952198893776 ...
|
||||
1.0170168994 2.14733589918 7 -1.3054866978 2.68803738466 ...
|
||||
0.723253652257 3.43552889347 7.5 1.8213700853 0.592593586195 ...
|
||||
1.24720806008 3.87383806577 7.5 0.0522300654168 0.988871238698 ...
|
||||
0.482531471239 2.67793287032 7.5 2.9693944293 -0.108591166081 ...
|
||||
0.154056100439 0.927269031704 6.75 0.119222057561 3.30489209451 ...
|
||||
0.0694865769274 6.65916526788 6.25 0.889014476084 -2.83976849035 ...
|
||||
-0.121267434867 0.341442615696 5.25 0.323053239216 -3.49289229012 ...
|
||||
0.726473690375 -3.5423730964 4 2.19149290449 -3.20855054004 ...
|
||||
1.39271709108 2.63121085718 3.75 0.88406577736 0.75622580197 ...
|
||||
1.07502077727 5.88578836799 4.25 -2.55088273352 2.89018116374 ...
|
||||
0.759049251607 4.24703604223 4.5 0.575687665685 -0.388292506167 ...
|
||||
];
|
||||
|
||||
data = reshape(data,5,86)';
|
||||
y_obs = data(:,1);
|
||||
pie_obs = data(:,2);
|
||||
R_obs = data(:,3);
|
||||
de = data(:,4);
|
||||
dq = data(:,5);
|
||||
|
||||
%Country: Canada
|
||||
%Sample Range: 1981:2 to 2002:3
|
||||
%Observations: 86
|
||||
%Variables: Real GDP Growth [%], Inflation [annualized %], Nom Rate [%],
|
||||
% Exchange Rate Change [%], Terms of Trade Change [%]
|
||||
|
|
|
@ -1,226 +1,226 @@
|
|||
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
|
||||
varexo e_R e_q e_ys e_pies e_A;
|
||||
|
||||
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
|
||||
|
||||
psi1 = 1.54;
|
||||
psi2 = 0.25;
|
||||
psi3 = 0.25;
|
||||
rho_R = 0.5;
|
||||
alpha = 0.3;
|
||||
rr = 2.51;
|
||||
k = 0.5;
|
||||
tau = 0.5;
|
||||
rho_q = 0.4;
|
||||
rho_A = 0.2;
|
||||
rho_ys = 0.9;
|
||||
rho_pies = 0.7;
|
||||
|
||||
|
||||
model(linear);
|
||||
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
|
||||
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
|
||||
pie = de+(1-alpha)*dq+pie_s;
|
||||
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
|
||||
dq = rho_q*dq(-1)+e_q;
|
||||
y_s = rho_ys*y_s(-1)+e_ys;
|
||||
pie_s = rho_pies*pie_s(-1)+e_pies;
|
||||
A = rho_A*A(-1)+e_A;
|
||||
y_obs = y-y(-1)+A;
|
||||
pie_obs = 4*pie;
|
||||
R_obs = 4*R;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_R = 1.25^2;
|
||||
var e_q = 2.5^2;
|
||||
var e_A = 1.89;
|
||||
var e_ys = 1.89;
|
||||
var e_pies = 1.89;
|
||||
end;
|
||||
|
||||
varobs y_obs R_obs pie_obs dq de;
|
||||
|
||||
estimated_params;
|
||||
psi1 , gamma_pdf,1.5,0.5;
|
||||
psi2 , gamma_pdf,0.25,0.125;
|
||||
psi3 , gamma_pdf,0.25,0.125;
|
||||
rho_R ,beta_pdf,0.5,0.2;
|
||||
alpha ,beta_pdf,0.3,0.1;
|
||||
rr ,gamma_pdf,2.5,1;
|
||||
k , gamma_pdf,0.5,0.25;
|
||||
tau ,gamma_pdf,0.5,0.2;
|
||||
rho_q ,beta_pdf,0.4,0.2;
|
||||
rho_A ,beta_pdf,0.5,0.2;
|
||||
rho_ys ,beta_pdf,0.8,0.1;
|
||||
rho_pies,beta_pdf,0.7,0.15;
|
||||
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
|
||||
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
|
||||
end;
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('NOW I DO STABILITY MAPPING and prepare sample for Reduced form Mapping');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(redform=1, //create sample of reduced form coefficients
|
||||
alpha2_stab=0.4,
|
||||
ksstat=0);
|
||||
// NOTE: since namendo is emppty by default,
|
||||
// this call does not perform the mapping of reduced form coefficient: just prepares the sample
|
||||
|
||||
/*
|
||||
disp(' ');
|
||||
disp('ANALYSIS OF REDUCED FORM COEFFICIENTS');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(load_stab=1, // loead previously generated sample analysed for stability
|
||||
redform=1, // do the reduced form mapping
|
||||
logtrans_redform=1, // estimate log-transformed reduced form coefficients (default=0)
|
||||
namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
|
||||
namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
|
||||
namlagendo=(R), // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
|
||||
stab=0 // don't repeat again the stability mapping
|
||||
);
|
||||
*/
|
||||
|
||||
disp(' ');
|
||||
disp('THE PREVIOUS TWO CALLS COULD BE DONE TOGETHER');
|
||||
disp('BY USING THE COMBINED CALL');
|
||||
disp(' ');
|
||||
disp('dynare_sensitivity(alpha2_stab=0.4, ksstat=0, redform=1,')
|
||||
disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R));')
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
//dynare_sensitivity(
|
||||
//alpha2_stab=0.4,
|
||||
//ksstat=0,
|
||||
//redform=1, //create sample of reduced form coefficients
|
||||
//logtrans_redform=1, // estimate log-transformed reduced form coefficients (default=0)
|
||||
//namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
|
||||
//namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
|
||||
//namlagendo=(R) // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
|
||||
//);
|
||||
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('MC FILTERING(rmse=1), TO MAP THE FIT FROM PRIORS');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible
|
||||
load_stab=1, // load prior sample
|
||||
istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error)
|
||||
stab=0, // don't plot again stability analysis results
|
||||
rmse=1, // do rmse analysis
|
||||
pfilt_rmse=0.1, // filtering criterion, i.e. filter the best 10 percent rmse's
|
||||
alpha2_rmse=0.3, // critical value for correlations in the rmse filterting analysis:
|
||||
// if ==1, means no corrleation analysis done
|
||||
alpha_rmse=1 // critical value for smirnov statistics of filtered samples
|
||||
// if ==1, no smornov analysis is done
|
||||
);
|
||||
|
||||
disp(' ');
|
||||
disp('THE PREVIOUS THREE CALLS COULD BE DONE TOGETHER');
|
||||
disp('BY USING THE COMBINED CALL');
|
||||
disp(' ');
|
||||
disp('dynare_sensitivity(alpha2_stab=0.4, ksstat=0, redform=1,')
|
||||
disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R),')
|
||||
disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,')
|
||||
disp('istart_rmse=2, rmse=1, pfilt_rmse=0.1, alpha2_rmse=0.3, alpha_rmse=1);')
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
//dynare_sensitivity(
|
||||
//alpha2_stab=0.4,
|
||||
//ksstat=0,
|
||||
//redform=1, //create sample of reduced form coefficients
|
||||
//logtrans_redform=1, // estimate log-transformed reduced form coefficients (default=0)
|
||||
//namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
|
||||
//namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
|
||||
//namlagendo=(R), // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
|
||||
//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
|
||||
//istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error)
|
||||
//rmse=1, // do rmse analysis
|
||||
//pfilt_rmse=0.1, // filtering criterion, i.e. filter the best 10 percent rmse's
|
||||
//alpha2_rmse=0.3, // critical value for correlations in the rmse filterting analysis:
|
||||
// // if ==1, means no corrleation analysis done
|
||||
//alpha_rmse=1 // critical value for smirnov statistics of filtered samples
|
||||
// // if ==1, no smirnov sensitivity analysis is done
|
||||
//);
|
||||
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('I ESTIMATE THE MODEL');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
// run this to generate posterior mode and Metropolis files if not yet done
|
||||
estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,
|
||||
prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=4, nograph, mh_drop=0.6,
|
||||
bayesian_irf, filtered_vars, smoother) y_obs R_obs pie_obs dq de;
|
||||
|
||||
|
||||
// run this to produce posterior samples of filtered, smoothed and irf variables, if not yet done
|
||||
//estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3,
|
||||
// mh_replic=0, mode_file=ls2003_mode, mode_compute=0, nograph, load_mh_file, bayesian_irf,
|
||||
// filtered_vars, smoother, mh_drop=0.6);
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('WE DO STABILITY MAPPING AGAIN, BUT FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE (or ML) and Hessian (pprior=0 & ppost=0)');
|
||||
disp('Typical for ML estimation, also feasible for posterior mode');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,
|
||||
mode_file=ls2003_mode // specifies the mode file where the mode and Hessian are stored
|
||||
);
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('RMSE ANALYSIS FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
dynare_sensitivity(mode_file=ls2003_mode,
|
||||
datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
|
||||
pprior=0,
|
||||
stab=0,
|
||||
rmse=1,
|
||||
alpha2_rmse=1, // no correlation analysis
|
||||
alpha_rmse=1 // no Smirnov sensitivity analysis
|
||||
);
|
||||
|
||||
disp(' ');
|
||||
disp('THE LAST TWO CALLS COULD BE DONE TOGETHER');
|
||||
disp('BY USING THE COMBINED CALL');
|
||||
disp(' ');
|
||||
disp('dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,')
|
||||
disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,')
|
||||
disp('rmse=1, alpha2_rmse=1, alpha_rmse=1);')
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
//dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,
|
||||
//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
|
||||
//rmse=1,
|
||||
//alpha2_rmse=1, // no correlation analysis
|
||||
//alpha_rmse=1 // no Smirnov sensitivity analysis
|
||||
//);
|
||||
|
||||
/*
|
||||
disp(' ');
|
||||
disp('RMSE ANALYSIS FOR POSTERIOR MCMC sample (ppost=1)');
|
||||
disp('Needs a call to dynare_estimation to load all MH environment');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2, mode_file=ls2003_mode, load_mh_file,
|
||||
prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=0, nograph, mh_drop=0.6);
|
||||
|
||||
dynare_sensitivity(stab=0, // no need for stability analysis since the posterior sample is surely OK
|
||||
rmse=1,ppost=1,alpha2_rmse=1,alpha_rmse=1);
|
||||
*/
|
||||
|
||||
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
|
||||
varexo e_R e_q e_ys e_pies e_A;
|
||||
|
||||
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
|
||||
|
||||
psi1 = 1.54;
|
||||
psi2 = 0.25;
|
||||
psi3 = 0.25;
|
||||
rho_R = 0.5;
|
||||
alpha = 0.3;
|
||||
rr = 2.51;
|
||||
k = 0.5;
|
||||
tau = 0.5;
|
||||
rho_q = 0.4;
|
||||
rho_A = 0.2;
|
||||
rho_ys = 0.9;
|
||||
rho_pies = 0.7;
|
||||
|
||||
|
||||
model(linear);
|
||||
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
|
||||
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
|
||||
pie = de+(1-alpha)*dq+pie_s;
|
||||
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
|
||||
dq = rho_q*dq(-1)+e_q;
|
||||
y_s = rho_ys*y_s(-1)+e_ys;
|
||||
pie_s = rho_pies*pie_s(-1)+e_pies;
|
||||
A = rho_A*A(-1)+e_A;
|
||||
y_obs = y-y(-1)+A;
|
||||
pie_obs = 4*pie;
|
||||
R_obs = 4*R;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_R = 1.25^2;
|
||||
var e_q = 2.5^2;
|
||||
var e_A = 1.89;
|
||||
var e_ys = 1.89;
|
||||
var e_pies = 1.89;
|
||||
end;
|
||||
|
||||
varobs y_obs R_obs pie_obs dq de;
|
||||
|
||||
estimated_params;
|
||||
psi1 , gamma_pdf,1.5,0.5;
|
||||
psi2 , gamma_pdf,0.25,0.125;
|
||||
psi3 , gamma_pdf,0.25,0.125;
|
||||
rho_R ,beta_pdf,0.5,0.2;
|
||||
alpha ,beta_pdf,0.3,0.1;
|
||||
rr ,gamma_pdf,2.5,1;
|
||||
k , gamma_pdf,0.5,0.25;
|
||||
tau ,gamma_pdf,0.5,0.2;
|
||||
rho_q ,beta_pdf,0.4,0.2;
|
||||
rho_A ,beta_pdf,0.5,0.2;
|
||||
rho_ys ,beta_pdf,0.8,0.1;
|
||||
rho_pies,beta_pdf,0.7,0.15;
|
||||
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
|
||||
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
|
||||
end;
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('NOW I DO STABILITY MAPPING and prepare sample for Reduced form Mapping');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(redform=1, //create sample of reduced form coefficients
|
||||
alpha2_stab=0.4,
|
||||
ksstat=0);
|
||||
// NOTE: since namendo is emppty by default,
|
||||
// this call does not perform the mapping of reduced form coefficient: just prepares the sample
|
||||
|
||||
/*
|
||||
disp(' ');
|
||||
disp('ANALYSIS OF REDUCED FORM COEFFICIENTS');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(load_stab=1, // loead previously generated sample analysed for stability
|
||||
redform=1, // do the reduced form mapping
|
||||
logtrans_redform=1, // estimate log-transformed reduced form coefficients (default=0)
|
||||
namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
|
||||
namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
|
||||
namlagendo=(R), // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
|
||||
stab=0 // don't repeat again the stability mapping
|
||||
);
|
||||
*/
|
||||
|
||||
disp(' ');
|
||||
disp('THE PREVIOUS TWO CALLS COULD BE DONE TOGETHER');
|
||||
disp('BY USING THE COMBINED CALL');
|
||||
disp(' ');
|
||||
disp('dynare_sensitivity(alpha2_stab=0.4, ksstat=0, redform=1,')
|
||||
disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R));')
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
//dynare_sensitivity(
|
||||
//alpha2_stab=0.4,
|
||||
//ksstat=0,
|
||||
//redform=1, //create sample of reduced form coefficients
|
||||
//logtrans_redform=1, // estimate log-transformed reduced form coefficients (default=0)
|
||||
//namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
|
||||
//namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
|
||||
//namlagendo=(R) // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
|
||||
//);
|
||||
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('MC FILTERING(rmse=1), TO MAP THE FIT FROM PRIORS');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(datafile=data_ca1,first_obs=8,nobs=79,prefilter=1, // also presample=2,loglinear, are admissible
|
||||
load_stab=1, // load prior sample
|
||||
istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error)
|
||||
stab=0, // don't plot again stability analysis results
|
||||
rmse=1, // do rmse analysis
|
||||
pfilt_rmse=0.1, // filtering criterion, i.e. filter the best 10 percent rmse's
|
||||
alpha2_rmse=0.3, // critical value for correlations in the rmse filterting analysis:
|
||||
// if ==1, means no corrleation analysis done
|
||||
alpha_rmse=1 // critical value for smirnov statistics of filtered samples
|
||||
// if ==1, no smornov analysis is done
|
||||
);
|
||||
|
||||
disp(' ');
|
||||
disp('THE PREVIOUS THREE CALLS COULD BE DONE TOGETHER');
|
||||
disp('BY USING THE COMBINED CALL');
|
||||
disp(' ');
|
||||
disp('dynare_sensitivity(alpha2_stab=0.4, ksstat=0, redform=1,')
|
||||
disp('logtrans_redform=1, namendo=(pie,R), namexo=(e_R), namlagendo=(R),')
|
||||
disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,')
|
||||
disp('istart_rmse=2, rmse=1, pfilt_rmse=0.1, alpha2_rmse=0.3, alpha_rmse=1);')
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
//dynare_sensitivity(
|
||||
//alpha2_stab=0.4,
|
||||
//ksstat=0,
|
||||
//redform=1, //create sample of reduced form coefficients
|
||||
//logtrans_redform=1, // estimate log-transformed reduced form coefficients (default=0)
|
||||
//namendo=(pie,R), // evaluate relationships for pie and R (namendo=(:) for all variables)
|
||||
//namexo=(e_R), // evaluate relationships with exogenous e_R (use namexo=(:) for all shocks)
|
||||
//namlagendo=(R), // evaluate relationships with lagged R (use namlagendo=(:) for all lagged endogenous)
|
||||
//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
|
||||
//istart_rmse=2, //start computing rmse from second observation (i.e. rmse does not inlude initial big error)
|
||||
//rmse=1, // do rmse analysis
|
||||
//pfilt_rmse=0.1, // filtering criterion, i.e. filter the best 10 percent rmse's
|
||||
//alpha2_rmse=0.3, // critical value for correlations in the rmse filterting analysis:
|
||||
// // if ==1, means no corrleation analysis done
|
||||
//alpha_rmse=1 // critical value for smirnov statistics of filtered samples
|
||||
// // if ==1, no smirnov sensitivity analysis is done
|
||||
//);
|
||||
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('I ESTIMATE THE MODEL');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
// run this to generate posterior mode and Metropolis files if not yet done
|
||||
estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,
|
||||
prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=4, nograph, mh_drop=0.6,
|
||||
bayesian_irf, filtered_vars, smoother) y_obs R_obs pie_obs dq de;
|
||||
|
||||
|
||||
// run this to produce posterior samples of filtered, smoothed and irf variables, if not yet done
|
||||
//estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2,prefilter=1,mh_jscale=0.3,
|
||||
// mh_replic=0, mode_file=ls2003_mode, mode_compute=0, nograph, load_mh_file, bayesian_irf,
|
||||
// filtered_vars, smoother, mh_drop=0.6);
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('WE DO STABILITY MAPPING AGAIN, BUT FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE (or ML) and Hessian (pprior=0 & ppost=0)');
|
||||
disp('Typical for ML estimation, also feasible for posterior mode');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
|
||||
dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,
|
||||
mode_file=ls2003_mode // specifies the mode file where the mode and Hessian are stored
|
||||
);
|
||||
|
||||
|
||||
disp(' ');
|
||||
disp('RMSE ANALYSIS FOR MULTIVARIATE SAMPLE AT THE POSTERIOR MODE');
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
dynare_sensitivity(mode_file=ls2003_mode,
|
||||
datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
|
||||
pprior=0,
|
||||
stab=0,
|
||||
rmse=1,
|
||||
alpha2_rmse=1, // no correlation analysis
|
||||
alpha_rmse=1 // no Smirnov sensitivity analysis
|
||||
);
|
||||
|
||||
disp(' ');
|
||||
disp('THE LAST TWO CALLS COULD BE DONE TOGETHER');
|
||||
disp('BY USING THE COMBINED CALL');
|
||||
disp(' ');
|
||||
disp('dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,')
|
||||
disp('datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,')
|
||||
disp('rmse=1, alpha2_rmse=1, alpha_rmse=1);')
|
||||
disp(' ');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
//dynare_sensitivity(pprior=0,Nsam=2048,alpha2_stab=0.4,mode_file=ls2003_mode,
|
||||
//datafile=data_ca1,first_obs=8,nobs=79,prefilter=1,
|
||||
//rmse=1,
|
||||
//alpha2_rmse=1, // no correlation analysis
|
||||
//alpha_rmse=1 // no Smirnov sensitivity analysis
|
||||
//);
|
||||
|
||||
/*
|
||||
disp(' ');
|
||||
disp('RMSE ANALYSIS FOR POSTERIOR MCMC sample (ppost=1)');
|
||||
disp('Needs a call to dynare_estimation to load all MH environment');
|
||||
disp('Press ENTER to continue'); pause(5);
|
||||
estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=2, mode_file=ls2003_mode, load_mh_file,
|
||||
prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=0, nograph, mh_drop=0.6);
|
||||
|
||||
dynare_sensitivity(stab=0, // no need for stability analysis since the posterior sample is surely OK
|
||||
rmse=1,ppost=1,alpha2_rmse=1,alpha_rmse=1);
|
||||
*/
|
||||
|
||||
|
|
|
@ -1,88 +1,88 @@
|
|||
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
|
||||
varexo e_R e_q e_ys e_pies e_A;
|
||||
|
||||
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
|
||||
|
||||
psi1 = 1.54;
|
||||
psi2 = 0.25;
|
||||
psi3 = 0.25;
|
||||
rho_R = 0.5;
|
||||
alpha = 0.3;
|
||||
rr = 2.51;
|
||||
k = 0.5;
|
||||
tau = 0.5;
|
||||
rho_q = 0.4;
|
||||
rho_A = 0.2;
|
||||
rho_ys = 0.9;
|
||||
rho_pies = 0.7;
|
||||
|
||||
|
||||
model(linear);
|
||||
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
|
||||
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
|
||||
pie = de+(1-alpha)*dq+pie_s;
|
||||
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
|
||||
dq = rho_q*dq(-1)+e_q;
|
||||
y_s = rho_ys*y_s(-1)+e_ys;
|
||||
pie_s = rho_pies*pie_s(-1)+e_pies;
|
||||
A = rho_A*A(-1)+e_A;
|
||||
y_obs = y-y(-1)+A;
|
||||
pie_obs = 4*pie;
|
||||
R_obs = 4*R;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_R = 1.25^2;
|
||||
var e_q = 2.5^2;
|
||||
var e_A = 1.89;
|
||||
var e_ys = 1.89;
|
||||
var e_pies = 1.89;
|
||||
end;
|
||||
|
||||
varobs y_obs R_obs pie_obs dq de;
|
||||
|
||||
estimated_params;
|
||||
psi1 , gamma_pdf,1.5,0.5;
|
||||
psi2 , gamma_pdf,0.25,0.125;
|
||||
psi3 , gamma_pdf,0.25,0.125;
|
||||
rho_R ,beta_pdf,0.5,0.2;
|
||||
alpha ,beta_pdf,0.3,0.1;
|
||||
rr ,gamma_pdf,2.5,1;
|
||||
k , gamma_pdf,0.5,0.25;
|
||||
tau ,gamma_pdf,0.5,0.2;
|
||||
rho_q ,beta_pdf,0.4,0.2;
|
||||
rho_A ,beta_pdf,0.5,0.2;
|
||||
rho_ys ,beta_pdf,0.8,0.1;
|
||||
rho_pies,beta_pdf,0.7,0.15;
|
||||
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
|
||||
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
|
||||
end;
|
||||
|
||||
|
||||
disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM A SCREENING FOR IDENTIFICATION ANALYSIS');
|
||||
disp('TYPE II ERRORS')
|
||||
disp(' ')
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
|
||||
dynare_sensitivity(identification=1, morris_nliv=6, morris_ntra=50);
|
||||
|
||||
disp('CREATE MC SAMPLE, CHECK FOR STABILITY AND PERFORM IDENTIFICATION ANALYSIS');
|
||||
disp('TYPE I (main effects)')
|
||||
disp(' ')
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
dynare_sensitivity(identification=1, morris=0);
|
||||
|
||||
disp('USE PREVIOUS MC SAMPLE AND PERFORM IDENTIFICATION ANALYSIS');
|
||||
disp('WIth analytic derivatives')
|
||||
disp(' ')
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
dynare_sensitivity(identification=1, load_stab=1, stab=0, morris=2);
|
||||
|
||||
|
||||
stoch_simul(order=1,irf=40);
|
||||
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
|
||||
varexo e_R e_q e_ys e_pies e_A;
|
||||
|
||||
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
|
||||
|
||||
psi1 = 1.54;
|
||||
psi2 = 0.25;
|
||||
psi3 = 0.25;
|
||||
rho_R = 0.5;
|
||||
alpha = 0.3;
|
||||
rr = 2.51;
|
||||
k = 0.5;
|
||||
tau = 0.5;
|
||||
rho_q = 0.4;
|
||||
rho_A = 0.2;
|
||||
rho_ys = 0.9;
|
||||
rho_pies = 0.7;
|
||||
|
||||
|
||||
model(linear);
|
||||
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
|
||||
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
|
||||
pie = de+(1-alpha)*dq+pie_s;
|
||||
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
|
||||
dq = rho_q*dq(-1)+e_q;
|
||||
y_s = rho_ys*y_s(-1)+e_ys;
|
||||
pie_s = rho_pies*pie_s(-1)+e_pies;
|
||||
A = rho_A*A(-1)+e_A;
|
||||
y_obs = y-y(-1)+A;
|
||||
pie_obs = 4*pie;
|
||||
R_obs = 4*R;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_R = 1.25^2;
|
||||
var e_q = 2.5^2;
|
||||
var e_A = 1.89;
|
||||
var e_ys = 1.89;
|
||||
var e_pies = 1.89;
|
||||
end;
|
||||
|
||||
varobs y_obs R_obs pie_obs dq de;
|
||||
|
||||
estimated_params;
|
||||
psi1 , gamma_pdf,1.5,0.5;
|
||||
psi2 , gamma_pdf,0.25,0.125;
|
||||
psi3 , gamma_pdf,0.25,0.125;
|
||||
rho_R ,beta_pdf,0.5,0.2;
|
||||
alpha ,beta_pdf,0.3,0.1;
|
||||
rr ,gamma_pdf,2.5,1;
|
||||
k , gamma_pdf,0.5,0.25;
|
||||
tau ,gamma_pdf,0.5,0.2;
|
||||
rho_q ,beta_pdf,0.4,0.2;
|
||||
rho_A ,beta_pdf,0.5,0.2;
|
||||
rho_ys ,beta_pdf,0.8,0.1;
|
||||
rho_pies,beta_pdf,0.7,0.15;
|
||||
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
|
||||
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
|
||||
end;
|
||||
|
||||
|
||||
disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM A SCREENING FOR IDENTIFICATION ANALYSIS');
|
||||
disp('TYPE II ERRORS')
|
||||
disp(' ')
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
|
||||
dynare_sensitivity(identification=1, morris_nliv=6, morris_ntra=50);
|
||||
|
||||
disp('CREATE MC SAMPLE, CHECK FOR STABILITY AND PERFORM IDENTIFICATION ANALYSIS');
|
||||
disp('TYPE I (main effects)')
|
||||
disp(' ')
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
dynare_sensitivity(identification=1, morris=0);
|
||||
|
||||
disp('USE PREVIOUS MC SAMPLE AND PERFORM IDENTIFICATION ANALYSIS');
|
||||
disp('WIth analytic derivatives')
|
||||
disp(' ')
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
dynare_sensitivity(identification=1, load_stab=1, stab=0, morris=2);
|
||||
|
||||
|
||||
stoch_simul(order=1,irf=40);
|
||||
|
|
|
@ -1,72 +1,72 @@
|
|||
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
|
||||
varexo e_R e_q e_ys e_pies e_A;
|
||||
|
||||
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
|
||||
|
||||
psi1 = 1.54;
|
||||
psi2 = 0.25;
|
||||
psi3 = 0.25;
|
||||
rho_R = 0.5;
|
||||
alpha = 0.3;
|
||||
rr = 2.51;
|
||||
k = 0.5;
|
||||
tau = 0.5;
|
||||
rho_q = 0.4;
|
||||
rho_A = 0.2;
|
||||
rho_ys = 0.9;
|
||||
rho_pies = 0.7;
|
||||
|
||||
|
||||
model(linear);
|
||||
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
|
||||
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
|
||||
pie = de+(1-alpha)*dq+pie_s;
|
||||
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
|
||||
dq = rho_q*dq(-1)+e_q;
|
||||
y_s = rho_ys*y_s(-1)+e_ys;
|
||||
pie_s = rho_pies*pie_s(-1)+e_pies;
|
||||
A = rho_A*A(-1)+e_A;
|
||||
y_obs = y-y(-1)+A;
|
||||
pie_obs = 4*pie;
|
||||
R_obs = 4*R;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_R = 1.25^2;
|
||||
var e_q = 2.5^2;
|
||||
var e_A = 1.89;
|
||||
var e_ys = 1.89;
|
||||
var e_pies = 1.89;
|
||||
end;
|
||||
|
||||
varobs y_obs R_obs pie_obs dq de;
|
||||
|
||||
estimated_params;
|
||||
psi1 , gamma_pdf,1.5,0.5;
|
||||
psi2 , gamma_pdf,0.25,0.125;
|
||||
psi3 , gamma_pdf,0.25,0.125;
|
||||
rho_R ,beta_pdf,0.5,0.2;
|
||||
alpha ,beta_pdf,0.3,0.1;
|
||||
rr ,gamma_pdf,2.5,1;
|
||||
k , gamma_pdf,0.5,0.25;
|
||||
tau ,gamma_pdf,0.5,0.2;
|
||||
rho_q ,beta_pdf,0.4,0.2;
|
||||
rho_A ,beta_pdf,0.5,0.2;
|
||||
rho_ys ,beta_pdf,0.8,0.1;
|
||||
rho_pies,beta_pdf,0.7,0.15;
|
||||
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
|
||||
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
|
||||
end;
|
||||
|
||||
|
||||
disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM SENSITIVITY ANALYSIS');
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
|
||||
dynare_sensitivity(morris=1, morris_nliv=6, morris_ntra=20, redform=1,
|
||||
namendo=(:), namexo=(:), namlagendo=(:));
|
||||
|
||||
stoch_simul(order=1,irf=40);
|
||||
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
|
||||
varexo e_R e_q e_ys e_pies e_A;
|
||||
|
||||
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
|
||||
|
||||
psi1 = 1.54;
|
||||
psi2 = 0.25;
|
||||
psi3 = 0.25;
|
||||
rho_R = 0.5;
|
||||
alpha = 0.3;
|
||||
rr = 2.51;
|
||||
k = 0.5;
|
||||
tau = 0.5;
|
||||
rho_q = 0.4;
|
||||
rho_A = 0.2;
|
||||
rho_ys = 0.9;
|
||||
rho_pies = 0.7;
|
||||
|
||||
|
||||
model(linear);
|
||||
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
|
||||
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
|
||||
pie = de+(1-alpha)*dq+pie_s;
|
||||
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
|
||||
dq = rho_q*dq(-1)+e_q;
|
||||
y_s = rho_ys*y_s(-1)+e_ys;
|
||||
pie_s = rho_pies*pie_s(-1)+e_pies;
|
||||
A = rho_A*A(-1)+e_A;
|
||||
y_obs = y-y(-1)+A;
|
||||
pie_obs = 4*pie;
|
||||
R_obs = 4*R;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var e_R = 1.25^2;
|
||||
var e_q = 2.5^2;
|
||||
var e_A = 1.89;
|
||||
var e_ys = 1.89;
|
||||
var e_pies = 1.89;
|
||||
end;
|
||||
|
||||
varobs y_obs R_obs pie_obs dq de;
|
||||
|
||||
estimated_params;
|
||||
psi1 , gamma_pdf,1.5,0.5;
|
||||
psi2 , gamma_pdf,0.25,0.125;
|
||||
psi3 , gamma_pdf,0.25,0.125;
|
||||
rho_R ,beta_pdf,0.5,0.2;
|
||||
alpha ,beta_pdf,0.3,0.1;
|
||||
rr ,gamma_pdf,2.5,1;
|
||||
k , gamma_pdf,0.5,0.25;
|
||||
tau ,gamma_pdf,0.5,0.2;
|
||||
rho_q ,beta_pdf,0.4,0.2;
|
||||
rho_A ,beta_pdf,0.5,0.2;
|
||||
rho_ys ,beta_pdf,0.8,0.1;
|
||||
rho_pies,beta_pdf,0.7,0.15;
|
||||
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
|
||||
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
|
||||
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
|
||||
end;
|
||||
|
||||
|
||||
disp('CREATE SCREENING SAMPLE, CHECK FOR STABILITY AND PERFORM SENSITIVITY ANALYSIS');
|
||||
disp('PRESS ENTER TO CONTUNUE');
|
||||
pause;
|
||||
|
||||
dynare_sensitivity(morris=1, morris_nliv=6, morris_ntra=20, redform=1,
|
||||
namendo=(:), namexo=(:), namlagendo=(:));
|
||||
|
||||
stoch_simul(order=1,irf=40);
|
||||
|
|
Loading…
Reference in New Issue