clear;
clc;

tic

%% X(1,:) == ppsipi value
%% X(2,:) == ppsiy value
%% X(3,:) == rrhor value
%% X(4,:) == ppibar value
%% X(5,:) == Mean Welfare
%% X(6,:) == 2*sigmaR - R_SS <0


ppsipi_min = 3;
ppsipi_max = 30;
ppsiy_min = 1;
ppsiy_max = 3;
rrhor_min = 0.5;
rrhor_max = 1;
ppibar_min = 1;
ppibar_max = 1.005;  %% 2% annual inflation

grid_ppsipi = 301;
grid_ppsiy = 31;
grid_rrhor = 11;
grid_ppibar = 21;

ppsipi_inc = (ppsipi_max-ppsipi_min)/(grid_ppsipi-1);
ppsiy_inc = (ppsiy_max-ppsiy_min)/(grid_ppsiy-1);
rrhor_inc = (rrhor_max-rrhor_min)/(grid_rrhor-1);
ppibar_inc = (ppibar_max-ppibar_min)/(grid_ppibar-1);

%%%%% Set up the grid values for ppsipi, ppsiy & rrhor

grid1 = ppsipi_min:ppsipi_inc:ppsipi_max;
grid2 = ppsiy_min:ppsiy_inc:ppsiy_max;
grid3 = rrhor_min:rrhor_inc:rrhor_max;
grid4 = ppibar_min:ppibar_inc:ppibar_max;

X = allcomb(grid1,grid2,grid3,grid4)';

X_out=zeros(2,size(X,2));

%first run
ppsipi=X(1,1);
ppsiy=X(2,1);
rrhor=X(3,1);
ppibar=X(4,1);

dynare Loop_dynare nolog noclearall;

%set required field
if isempty(options_.qz_criterium)
    options_.qz_criterium = 1+1e-6;
end

%run loop
for ii = 2:size(X,2)
    [X_out(ii)]=wrapper_function(X(:,ii),M_,options_,oo_,var_list_);       
end


