% grid RHO_R
Gr_RHO_R=0.1;
% grid ALPHA_Y
Gr_ALPHA_Y=0.1;
% grid ALPHA_PIC
Gr_ALPHA_PIC=0.1;
% grid ALPHA_FIN
Gr_ALPHA_FIN=0.1;


rhoy = [0.00:Gr_ALPHA_Y:3.00];
rhoi = [0.00:Gr_RHO_R:1.00];
rhop = [0.00:Gr_ALPHA_PIC:3.00];
rhofin = [0.00:Gr_ALPHA_FIN:1.00];

% Loop sobre los parámetros
    for t=1:length(rhoi),
        for x=1:length(rhop),
            for j=1:length(rhoy),
                for m=1:length(rhofin),
                RHO_R   = rhoi(1,t);            
                ALPHA_PI= rhop(1,x);
                ALPHA_Y = rhoy(1,j);
                ALPHA_F = rhofin(1,m);
                save parameters_mat RHO_R ALPHA_PI ALPHA_Y ALPHA_F;
                dynare NKBGGGK  noclearall; 
                if info > 0
                IND_R(j,x,t,m) = RHO_R;
                IND_PI(j,x,t,m)= ALPHA_PI;
                IND_Y(j,x,t,m) = ALPHA_Y;
                IND_F(j,x,t,m) = ALPHA_F;
                WELFAREcon(j,x,t,m) =-1e12;
                WELFAREunc(j,x,t,m) =-1e12;
                DETERMINACY(j,x,t,m)=-1e12;
                else
                IND_R(j,x,t,m) = RHO_R;
                IND_PI(j,x,t,m)= ALPHA_PI;
                IND_Y(j,x,t,m) = ALPHA_Y;
                IND_F(j,x,t,m) = ALPHA_F;
                WELFAREcon(j,x,t,m) = 0.5*oo_.dr.ghs2(oo_.dr.inv_order_var(1),1)+oo_.dr.ys(oo_.dr.inv_order_var(1),1);
                WELFAREunc(j,x,t,m) = oo_.mean(oo_.dr.inv_order_var(1),1);
                DETERMINACY(j,x,t,m)= 0;
                clear parameters_mat
                end;
            end;
        end; 
    end;
end;

Wmax=max(WELFAREcon);
    for t=1:length(rhoi)
        for x=1:length(rhop)
            for j=1:length(rhoy)
                for m=1:length(rhofin)
                    if WELFAREcon(j,x,t,m)>=Wmax
                        Wmax=WELFAREcon(j,x,t,m);
                        d1max=j;
                        d2max=x;
                        d3max=t;
                        d4max=m;
                    end;
                end;
            end;
        end;
    end;

param_optimos=[IND_R(d1max,d2max,d3max,d4max), IND_PI(d1max,d2max,d3max,d4max), IND_Y(d1max,d2max,d3max,d4max), IND_F(d1max,d2max,d3max,d4max)];