% %% pre run dynare to load the macro (this should be done only once)
clear all 
close all
clc
dynare sw4;
% 
%% setting the value for the parameter:

close all
clear all
clc

loadtheta

%% Fabio's combinations

fabio(1).subset = strvcat('pi_o','w_o ','r_o ','y_o ');
fabio(2).subset = strvcat('pi_o','i_o ','r_o ','y_o ');
fabio(3).subset = strvcat('pi_o','r_o ','y_o ','h_o ');
fabio(4).subset = strvcat('pi_o','r_o ','c_o ','h_o ');
fabio(5).subset = strvcat('pi_o','w_o ','r_o ','c_o ');
fabio(6).subset = strvcat('pi_o','y_o ','r_o ','c_o ');

delete('results_sw4_unrestricted.out');
% ouput file
diary('results_sw4_unrestricted.out');

table_comb_fabio = zeros(6,8);   
itheta           = 1;
combins          = [];
theta_0          = theta';
wb = waitbar(0, 'Computing the rank test for 6 combinations');
for co =1:6
    waitbar(co/6, wb);
%     disp(['Combination of observable ', num2str(co)]);  
%     disp(subsetn)
    % compute the ABCD matrices. they are in minimal state form solution, see
    % complementary documents in Kumunijer and Ng (2011).
    [minA,minB,minC,minD,Sigma] = sw4_abcd(theta_0,fabio(co).subset);
    % Compute the derivatives
    [Delta_0,Delta_L,Delta_T,Delta_U] = ...
        delta_as(theta_0,minA,minB,minC,minD,Sigma,fabio(co).subset);
    n_x           = size(minA,1);
    n_eps         = size(minB,2);
    n_y           = size(minC,1);
    n_theta       = size(theta_0,1);
    print_details = 0;
    tol0          = 1e-2;
    % prints the results in a table for a range of tolerance levels 
    table_comb_fabio(co,:) = rank_test_comb(Delta_0,n_theta,n_x,n_eps,n_y,tol0);    
end
close(wb);

table_n    = [n_theta n_x^2 n_eps^2  n_theta+n_x^2 n_theta+n_eps^2 n_theta+n_x^2+n_eps^2  1];
disp(sprintf('\n\n'));
disp('** Summary **');
disp(sprintf('n_theta = %d  n_x= %d  n_eps = %d',n_theta,n_x,n_eps));
disp(sprintf('Order condition: n_theta %d n_delta %d',...
	     n_theta,(n_y*n_x + n_eps*(n_x+n_y-n_eps) + n_eps*(n_eps+1)/2)));
fmt0=['%4e &' repmat('%5d &',1,cols(table_comb_fabio)-2) ' %5d' ' \\\\'];
fmt1=['%s&', repmat('%5d &',1, cols(table_n)-1) ' %5d' ' \\\\'] ;

disp(sprintf(' tol             DL     DT    DU   LT    LU    LTU    pass ' ));
for i=1:rows(table_comb_fabio);
    disp(sprintf(fmt0,table_comb_fabio(i,:)));
end;
disp(sprintf(fmt1,'  Require     ', table_n(1,:)));


%% using a subset of variables
%% all combination

close all
clear all
clc

loadtheta

all_variabs    = strvcat('y_o','c_o','i_o','w_o','pi_o','r_o','h_o');
all_nvariabs   = {'y_o','c_o','i_o','w_o','pi_o','r_o','h_o'};
[comb4]        = createcombinations(all_variabs,all_nvariabs);

disp('                            ');
disp('                            ');    
disp('                            ');  
    
theta_0    = theta';
table_comb = zeros(35,8);   
combins    = [];
wb = waitbar(0, 'Computing the rank test for 35 combinations');
for co =1:35
    waitbar(co/35, wb);
    subset  = comb4(co).var; 
    subsetn = comb4(co).nvar; 
    combins = [combins; subsetn];
    % compute the ABCD matrices. they are in minimal state form solution, see
    % complementary documents in Kumunijer and Ng (2011).
    [minA,minB,minC,minD,Sigma] = sw4_abcd(theta_0,subset);
    % Compute the derivatives
    [Delta_0,Delta_L,Delta_T,Delta_U] = ...
        delta_as(theta_0,minA,minB,minC,minD,Sigma,subset);
    n_x           = size(minA,1);
    n_eps         = size(minB,2);
    n_y           = size(minC,1);
    n_theta       = size(theta_0,1);
    tol0          = 1e-2;
    % prints the results in a table for a range of tolerance levels 
    table_comb(co,:) = rank_test_comb(Delta_0,n_theta,n_x,n_eps,n_y,tol0);
end
close(wb);


table_n    = [n_theta n_x^2 n_eps^2  n_theta+n_x^2 n_theta+n_eps^2 n_theta+n_x^2+n_eps^2  1];
disp(sprintf('\n\n'));
disp('** Summary **');
disp(sprintf('n_theta = %d  n_x= %d  n_eps = %d',n_theta,n_x,n_eps));
disp(sprintf('Order condition: n_theta %d n_delta %d',...
	     n_theta,(n_y*n_x + n_eps*(n_x+n_y-n_eps) + n_eps*(n_eps+1)/2)));
fmt0=['%4e &' repmat('%5d &',1,cols(table_comb)-2) ' %5d' ' \\\\'];
fmt1=['%s&', repmat('%5d &',1, cols(table_n)-1) ' %5d' ' \\\\'] ;

disp(sprintf(' tol             DL     DT    DU   LT    LU    LTU    pass ' ));
for i=1:rows(table_comb);
    disp(sprintf(fmt0,table_comb(i,:)));
end;
disp(sprintf(fmt1,'  Require     ', table_n(1,:)));

disp(sprintf('Combinations' ));
disp(combins);

diary off;

