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 = 0;
ppsipi_max = 30;
ppsiy_min = 0;
ppsiy_max = 3; 
rrhor_min = 0;
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 = [X;zeros(1,length(X(1,:)));zeros(1,length(X(1,:)))];

dynare Loop_optimal_switching nolog noclearall;
    
Xtemp = X(:,X(6,:)>0); %%filters out ZLB condition
[num] = max(Xtemp(5,:));
[rowloc] = find(Xtemp(5,:)==num);
ans = Xtemp(:,rowloc)
X =[ans X];

dlmwrite('output.csv',X','precision',10)
toc