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 = [X;zeros(1,length(X(1,:)));zeros(1,length(X(1,:)))];

for ii = 1:length(X(1,:))

    set_param_value('ppsipi', X(1,ii));
    set_param_value('ppsiy', X(2,ii));
    set_param_value('rrhor', X(3,ii));
    set_param_value('ppibar', X(4,ii));
		
		try %if indeterminate instead of stopping move on to catch
		
			dynare Loop_dynare nolog noclearall;
		
		  X(5,ii) = oo_.dr.ys(10) + oo_.dr.ghs2(oo_.dr.inv_order_var(10))/2; %%my definition of welfare instead of using oo_.mean(10) which is stochastic SS of Welf.
      X(6,ii) = oo_.mean(3)- 2*sqrt(oo_.var(3,3));   %%I say this condition on R must be satisfied for welfare to count. 
      
   catch
     	X(5,ii) = -Inf;  %%set welfare to NAN if indeterminate and move onto next loop
     	
    end

    ii

end
    
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