clear all;
clc;

kg_min = 0;
kg_max = 1;
kp_min = 0;
kp_max = 1;


leng = 101;
jumpg = (kg_max - kg_min)/(leng-1);
jumpp = (kp_max - kp_min)/(leng-1);

E = zeros(3,leng^2);

Temp = kg_min:jumpg:kg_max;
Temp = repmat(Temp,1,leng);
E(1,:) = Temp;

Temp = (kp_min:jumpp:kp_max)';
Temp = repmat(Temp,1,leng)';
Temp = Temp(:)';
E(2,:) = Temp';

first_time = 1;

for i=1:leng^2
    
      if first_time
            
      kappag=E(1,i);
      kappap=E(2,i);
      save param kappag kappap;
      dynare det6 noclearall;
      first_time = 0;
      
      else
      set_param_value('kappag',E(1,i));
      set_param_value('kappap',E(2,i));
      save param kappag kappap;
      perfect_foresight_setup;
      perfect_foresight_solver;
      end
      
y_path = oo_.endo_simul;
cc = y_path(strmatch('c',M_.endo_names,'exact'),:);
gg = y_path(strmatch('G',M_.endo_names,'exact'),:);
ll = y_path(strmatch('L',M_.endo_names,'exact'),:);
yy = y_path(strmatch('Y',M_.endo_names,'exact'),:);

bett = beta.^([0:1:length(yy)-1])';
welfare_loss = fsolve(@(x) bett.*((kappac * cc.^((v - 1) / v) + (1 - kappac) * gg.^((v - 1) / v)).^(v/(v - 1)) - khi * ll.^(1+phil)./(1+phil) ...
                   - ((kappac * (cc(1)*(1-x/100))^((v - 1) / v) + (1 - kappac) * gg(1)^((v - 1) / v))^(v/(v - 1)) - khi * ll(1)^(1+phil)./(1+phil))),0)


            
E(3,i) = abs(welfare_loss);
end

[A,I] = min(E(3,:));
disp('Min welfare loss is '); disp(A);
disp('Corresponding kappag is '); disp(E(1,I));
disp('Corresponding kappap is '); disp(E(2,I));
