sources of the distributed toolbox.

git-svn-id: file:///home/sebastien/dynare/gsa_dyn@23 f1850c17-3b45-254b-b221-fcb05880fee1

GSA_distrib/Morris_Measure_Groups.m

+
+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

GSA_distrib/Sampling_Function_2.m

+function [Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
+%[Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
+%	Inputs: k (1,1)                      := number of factors examined or number of groups examined.
+%                                           In case the groups are chosen the number of factors is stores in NumFact and
+%                                           sizea becomes the number of created groups. 
+%           NumFact (1,1)                := number of factors examined in the case when groups are chosen
+%	    	r (1,1)                      := sample size  
+%           p (1,1)                      := number of intervals considered in [0, 1]
+%           UB(sizea,1)                  := Upper Bound for each factor 
+%           LB(sizea,1)                  := Lower Bound for each factor 
+%           GroupNumber(1,1)             := Number of groups (eventually 0)
+%           GroupMat(NumFact,GroupNumber):= Matrix which describes the chosen groups. Each column represents a group and its elements 
+%                                           are set to 1 in correspondence of the factors that belong to the fixed group. All
+%                                           the other elements are zero.
+%   Local Variables:  
+%	    	sizeb (1,1)         := sizea+1
+%           sizec (1,1)         := 1
+%           randmult (sizea,1)  := vector of random +1 and -1  
+%           perm_e(1,sizea)     := vector of sizea random permutated indeces    
+%           fact(sizea)         := vector containing the factor varied within each traj
+% 	        DDo(sizea,sizea)    := D*       in Morris, 1991   
+%	        A(sizeb,sizea)      := Jk+1,k   in Morris, 1991
+%	        B(sizeb,sizea)      := B        in Morris, 1991
+%	        Po(sizea,sizea)     := P*       in Morris, 1991
+%           Bo(sizeb,sizea)     := B*       in Morris, 1991
+%	        Ao(sizeb,sizec)     := Jk+1,1   in Morris, 1991
+%	        xo(sizec,sizea)     := x*       in Morris, 1991 (starting point for the trajectory)
+%           In(sizeb,sizea)     := for each loop orientation matrix. It corresponds to a trajectory
+%                                  of k step in the parameter space and it provides a single elementary
+%                                  effect per factor 
+%           MyInt()
+%           Fact(sizea,1)       := for each loop vector indicating which factor or group of factors has been changed 
+%                                  in each step of the trajectory
+%           AuxMat(sizeb,sizea) := Delta*0.5*((2*B - A) * DD0 + A) in Morris, 1991. The AuxMat is used as in Morris design
+%                                  for single factor analysis, while it constitutes an intermediate step for the group analysis.
+%
+%	Output: Outmatrix(sizeb*r, sizea) := for the entire sample size computed In(i,j) matrices
+%           OutFact(sizea*r,1)        := for the entire sample size computed Fact(i,1) vectors
+%           
+%   Note: B0 is constructed as in Morris design when groups are not considered. When groups are considered the routine
+%         follows the following steps:
+%           1- Creation of P0 and DD0 matrices defined in Morris for the groups. This means that the dimensions of these
+%              2 matrices are (GroupNumber,GroupNumber).
+%           2- Creation of AuxMat matrix with (GroupNumber+1,GroupNumber) elements.
+%           3- Definition of GroupB0 starting from AuxMat, GroupMat and P0.
+%           4- The final B0 for groups is obtained as [ones(sizeb,1)*x0' + GroupB0]. The P0 permutation is present in GroupB0
+%              and it's not necessary to permute the matrix (ones(sizeb,1)*x0') because it's already randomly created. 
+%   Reference:
+%   A. Saltelli, K. Chan, E.M. Scott, "Sensitivity Analysis" on page 68 ss
+%
+%   F. Campolongo, J. Cariboni, JRC - IPSC Ispra, Varese, IT
+%   Last Update: 15 November 2005 by J.Cariboni
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% Parameters and initialisation of the output matrix
+sizea = k;
+Delta = p/(2*p-2);
+%Delta = 1/3
+NumFact = sizea;
+GroupNumber = size(GroupMat,2);
+
+if GroupNumber ~ 0;
+    sizea = size(GroupMat,2);
+end    
+
+sizeb = sizea + 1;
+sizec = 1;
+Outmatrix = [];
+OutFact = [];
+
+% For each i generate a trajectory  
+for i=1:r
+    
+    % Construct DD0 - OLD VERSION - it does not need communication toolbox
+    % RAND(N,M) is an NXM matrix with random entries, chosen from a uniform distribution on the interval (0.0,1.0).
+    % Note that DD0 tells if the factor have to be increased or ddecreased
+    % by Delta.
+    randmult = ones(k,1);           
+    v = rand(k,1);                  
+    randmult (find(v < 0.5))=-1;
+    randmult = repmat(randmult,1,k);
+    DD0 = randmult .* eye(k);
+    
+    % Construct DD0 - NEW VERSION - it needs communication toolbox
+    % randsrc(m) generates an m-by-m matrix, each of whose entries independently takes the value -1 with probability 1/2,
+    % and 1 with probability 1/2.
+    % DD0 = randsrc(NumFact) .* eye(NumFact);      
+    
+    % Construct B (lower triangular)
+    B = ones(sizeb,sizea);
+    for j = 1:sizea
+       B(1:j,j)=0;    
+    end
+    
+    % Construct A0, A
+    A0 = ones(sizeb,1);
+    A = ones(sizeb,NumFact);
+
+    % Construct the permutation matrix P0. In each column of P0 one randomly chosen element equals 1
+    % while all the others equal zero. 
+    % P0 tells the order in which order factors are changed in each
+    % trajectory. P0 is created as it follows:
+    % 1) All the elements of P0 are set equal to zero ==> P0 = zeros (sizea, sizea);
+    % 2) The function randperm create a random permutation of integer 1:sizea, without repetitions ==> perm_e; 
+    % 3) In each column of P0 the element indicated in perm_e is set equal to one.    
+    % Note that P0 is then used reading it by rows. 
+    P0 = zeros (sizea, sizea);
+    perm_e = randperm(sizea);               % RANDPERM(n) is a random permutation of the integers from 1 to n.
+    for j = 1:sizea
+        P0(perm_e(j),j) = 1;    
+    end    
+    
+    % When groups are present the random permutation is done only on B. The effect is the same since 
+    % the added part (A0*x0') is completely random. 
+    if GroupNumber ~ 0
+        B = B * (GroupMat*P0')';
+    end
+    
+    % Compute AuxMat both for single factors and groups analysis. For Single factors analysis
+    % AuxMat is added to (A0*X0) and then permutated through P0. When groups are active AuxMat is
+    % used to build GroupB0. AuxMat is created considering DD0. If the element on DD0 diagonal
+    % is 1 then AuxMat will start with zero and add Delta. If the element on DD0 diagonal is -1 
+    % then DD0 will start Delta and goes to zero.
+    AuxMat = Delta*0.5*((2*B - A) * DD0 + A);
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % a --> Define the random vector x0 for the factors. Note that x0 takes value in the hypercube
+    % [0,...,1-Delta]*[0,...,1-Delta]*[0,...,1-Delta]*[0,...,1-Delta] 
+    MyInt = repmat([0:(1/(p-1)):(1-Delta)],NumFact,1);     % Construct all possible values of the factors
+    
+    % OLD VERSION - it needs communication toolbox
+    % w = randint(NumFact,1,[1,size(MyInt,2)]);              
+    
+    % NEW VERSION - construct a version of random integers
+    % 1) create a vector of random integers
+    % 2) divide [0,1] into the needed steps
+    % 3) check in which interval the random numbers fall
+    % 4) generate the corresponding integer
+    v = repmat(rand(NumFact,1),1,size(MyInt,2)+1);     % 1)
+    IntUsed = repmat([0:1/size(MyInt,2):1],NumFact,1); % 2)
+    DiffAuxVec = IntUsed - v;                          % 3)
+    
+    for ii = 1:size(DiffAuxVec,1)
+        w(1,ii) = max(find(DiffAuxVec(ii,:)<0));       % 4)
+    end
+    x0 = MyInt(1,w)';                                  % Define x0    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
+    % b --> Compute the matrix B*, here indicated as B0. Each row in B0 is a
+    % trajectory for Morris Calculations. The dimension of B0 is (Numfactors+1,Numfactors) 
+    if GroupNumber ~ 0
+        B0 = (A0*x0' + AuxMat);
+    else
+        B0 = (A0*x0' + AuxMat)*P0;
+    end
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
+    % c --> Compute values in the original intervals
+    % B0 has values x(i,j) in [0, 1/(p -1), 2/(p -1), ... , 1].
+    % To obtain values in the original intervals [LB, UB] we compute
+    % LB(j) + x(i,j)*(UB(j)-LB(j))
+    In = repmat(LB,1,sizeb)' + B0 .* repmat((UB-LB),1,sizeb)';
+
+    % Create the Factor vector. Each component of this vector indicate which factor or group of factor
+    % has been changed in each step of the trajectory.
+    for j=1:sizea
+        Fact(1,j) = find(P0(j,:));
+    end
+    Fact(1,sizea+1) = 0;
+    
+    Outmatrix = [Outmatrix; In];
+    OutFact = [OutFact; Fact'];
+    
+end
\ No newline at end of file

GSA_distrib/cumplot.m

+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

GSA_distrib/dat_fil_.m

+function c = dat_fil_(data_file);
+%
+% Part of the Sensitivity Analysis Toolbox for DYNARE
+%
+% Written by Marco Ratto, 2006
+% Joint Research Centre, The European Commission,
+% (http://eemc.jrc.ec.europa.eu/),
+% marco.ratto@jrc.it 
+%
+% Disclaimer: This software is not subject to copyright protection and is in the public domain. 
+% It is an experimental system. The Joint Research Centre of European Commission 
+% assumes no responsibility whatsoever for its use by other parties
+% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
+% characteristic. We would appreciate acknowledgement if the software is used.
+% Reference:
+% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
+%
+
+try
+  eval(data_file);
+catch
+  load(data_file);
+end
+clear data_file;
+
+a=who;
+
+for j=1:length(a)
+  eval(['c.',a{j},'=',a{j},';']);
+end
\ No newline at end of file

GSA_distrib/dynare_MC.m

+function dynare_MC(var_list_,OutDir)
+%
+% Adapted by M. Ratto from dynare_estimation.m and posteriorsmoother.m
+% (dynare_estimation.m and posteriorsmoother.m are part of DYNARE,
+% copyright M. Juillard)
+%
+% Part of the Sensitivity Analysis Toolbox for DYNARE
+%
+% Written by Marco Ratto, 2006
+% Joint Research Centre, The European Commission,
+% (http://eemc.jrc.ec.europa.eu/),
+% marco.ratto@jrc.it 
+%
+% Disclaimer: This software is not subject to copyright protection and is in the public domain. 
+% It is an experimental system. The Joint Research Centre of European Commission 
+% assumes no responsibility whatsoever for its use by other parties
+% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
+% characteristic. We would appreciate acknowledgement if the software is used.
+% Reference:
+% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
+%
+
+global M_ options_ oo_ estim_params_ 
+global bayestopt_
+
+if options_.filtered_vars ~= 0 & options_.filter_step_ahead == 0
+  options_.filter_step_ahead = 1;
+end
+if options_.filter_step_ahead ~= 0
+  options_.nk = max(options_.filter_step_ahead);
+else
+  options_.nk = 0;
+end
+
+
+nvx = estim_params_.nvx;
+nvn = estim_params_.nvn;
+ncx = estim_params_.ncx;
+ncn = estim_params_.ncn;
+np  = estim_params_.np ;
+npar  = nvx+nvn+ncx+ncn+np;
+
+if isempty(options_.datafile)
+  error('ESTIMATION: datafile option is missing')
+end
+
+if isempty(options_.varobs)
+  error('ESTIMATION: VAROBS is missing')
+end
+
+%% Read and demean data 
+rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
+
+options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
+gend = options_.nobs;
+
+rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
+if options_.loglinear == 1 & ~options_.logdata
+  rawdata = log(rawdata);
+end
+if options_.prefilter == 1
+  bayestopt_.mean_varobs = mean(rawdata,1);
+  data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
+else
+  data = transpose(rawdata);
+end
+
+if ~isreal(rawdata)
+  error(['There are complex values in the data. Probably  a wrong' ...
+	 ' transformation'])
+end
+
+offset = npar-np;
+fname_=M_.fname;
+
+options_ = set_default_option(options_,'opt_gsa',1);
+options_gsa_ = options_.opt_gsa;
+
+if options_gsa_.pprior,
+  namfile=[fname_,'_prior'];
+else
+  namfile=[fname_,'_mc'];
+end
+load([OutDir,'\',namfile],'lpmat', 'lpmat0', 'istable')
+load(options_.mode_file)
+%%
+%%
+%%
+x=[lpmat0(istable,:) lpmat(istable,:)];
+clear lpmat lpmat0 istable %iunstable egg yys T
+B = size(x,1);
+[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(x(1,:)',gend,data);
+n1=size(atT,1);
+
+nfil=B/40;
+stock_smooth = zeros(M_.endo_nbr,gend,40);
+stock_filter = zeros(M_.endo_nbr,gend+1,40);
+stock_ys = zeros(40, M_.endo_nbr);
+logpo2=zeros(B,1);
+%%
+h = waitbar(0,'MC smoother ...');
+delete([OutDir,'\',namfile,'_*.mat'])
+ib=0;
+ifil=0;
+opt_gsa=options_.opt_gsa;
+for b=1:B
+  ib=ib+1;
+  deep = x(b,:)';
+  set_all_parameters(deep);
+  dr = resol(oo_.steady_state,0);
+  %deep(1:offset) = xparam1(1:offset);
+  logpo2(b,1) = DsgeLikelihood(deep,gend,data);
+  if opt_gsa.lik_only==0,
+  [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(deep,gend,data);
+  stock_smooth(:,:,ib)=atT(1:M_.endo_nbr,:);
+  stock_filter(:,:,ib)=filtered_state_vector(1:M_.endo_nbr,:);
+  %stock_filter(:,:,ib)=aK(options_.filter_step_ahead,1:M_.endo_nbr,:);
+  stock_ys(ib,:)=ys';
+  if ib==40,
+    ib=0;
+    ifil=ifil+1;
+    save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
+    stock_smooth = zeros(M_.endo_nbr,gend,40);
+    stock_filter = zeros(M_.endo_nbr,gend+1,40);
+    stock_ys = zeros(40, M_.endo_nbr);
+  end
+  end  
+  waitbar(b/B,h,['MC smoother ...',num2str(b),'/',num2str(B)]);
+end
+close(h)
+if opt_gsa.lik_only==0,
+if ib>0,
+    ifil=ifil+1;
+    stock_smooth = stock_smooth(:,:,1:ib);
+    stock_filter = stock_filter(:,:,1:ib);
+    stock_ys = stock_ys(1:ib,:);
+    save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys')
+end
+end
+stock_gend=gend;
+stock_data=data;
+save([OutDir,'\',namfile],'x','logpo2','stock_gend','stock_data','-append')

GSA_distrib/ghx2transition.m

+function [A,B] = ghx2transition(mm,iv,ic,aux)
+% [A,B] = ghx2transition(mm,iv,ic,aux)
+%
+% Adapted by M. Ratto from kalman_transition_matrix.m 
+% (kalman_transition_matrix.m is part of DYNARE, copyright M. Juillard)
+%
+% Part of the Sensitivity Analysis Toolbox for DYNARE
+%
+% Written by Marco Ratto, 2006
+% Joint Research Centre, The European Commission,
+% (http://eemc.jrc.ec.europa.eu/),
+% marco.ratto@jrc.it 
+%
+% Disclaimer: This software is not subject to copyright protection and is in the public domain. 
+% It is an experimental system. The Joint Research Centre of European Commission 
+% assumes no responsibility whatsoever for its use by other parties
+% and makes no guarantees, expressed or implied, about its quality, reliability, or any other
+% characteristic. We would appreciate acknowledgement if the software is used.
+% Reference:
+% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
+%
+
+global M_
+
+  [nr1, nc1] = size(mm);
+  ghx = mm(:, [1:(nc1-M_.exo_nbr)]);
+  ghu = mm(:, [(nc1-M_.exo_nbr+1):end] );
+  n_iv = length(iv);
+  n_ir1 = size(aux,1);
+  nr = n_iv + n_ir1;
+  
+  A = zeros(nr,nr);
+  B = zeros(nr,M_.exo_nbr);
+  
+  i_n_iv = 1:n_iv;
+  A(i_n_iv,ic) = ghx(iv,:);
+  if n_ir1 > 0
+    A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr);
+  end
+  
+  B(i_n_iv,:) = ghu(iv,:);

GSA_distrib/log_trans_.m

+function [yy, xdir, isig, lam]=log_trans_(y0,xdir0)