clear all;
clc;

phiinf_min = 0;
phiinf_max = 1;
phier_min = 0;
phier_max = 1; 

grid = 3;
phiinf_inc = (phiinf_max-phiinf_min)/(grid-1);
phier_inc = (phier_max-phier_min)/(grid-1);
X = zeros(14,grid^2);

% Set up the values for phiinf and ppibar

Temp = phiinf_min:phiinf_inc:phiinf_max;
Temp = repmat(Temp,1,grid);
X(1,:) = Temp;

Temp = (phier_min:phier_inc:phier_max)';
Temp = repmat(Temp,1,grid)';
Temp = Temp(:)';
X(2,:) = Temp';

first_time = 1;

for i=1:grid^2
    
      if first_time;
            
      phiinf=X(1,i);
      phier=X(2,i);
      save param_input phiinf phier;
      dynare Model10_rer_short noclearall;
      first_time = 0;
        
      else
      phiinf=X(1,i);
      phier=X(2,i);
      save param_input phiinf phier;
      stoch_simul(var_list_);
      end
        
welfare(1)=(1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,1)^2+...
    (sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,1)^2-...
    nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,1)^2+...
    (bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,1)^2-...
    gcy/(2*cy)*(1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(8,1)^2-...
    giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(26,1)^2+...
    bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,2)^2+...
    (bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,2)^2+...
    ((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,3)^2;

welfare(2)=welfare(1)+1/(1+iss)*((1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,2)^2+...
    (sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,2)^2-...
    nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,2)^2+...
    (bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,2)^2-...
    gcy/(2*cy)*(1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(8,2)^2-...
    giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(26,2)^2+...
    bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,3)^2+...
    (bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,3)^2+...
    ((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,4)^2)+...
    (bfy*(1-delta)/cy*(bfy*(1-delta)/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-alpha*sigma*nss^phi/((css-nss^phi/phi)*(1-alpha))+(delta-1/2))-nss^phi/css*((alpha/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi)-alpha/(2*(1-alpha))))*oo_.endo_simul(6,1)^2-...
    (nss^phi/css*((psi/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-psi/(2*(1-alpha))))*oo_.endo_simul(9,1)^2+...
    (bfy/(2*cy*betaf)+(inty*a2)^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))+inty/(2*cy)*(a2+a2^2)+(inty*a2)^2/(nfay*cy)-nss^phi/css*(eps*(1-alpha+alpha*eps)+sigma*nss^phi/(css-nss^phi/phi)+phi-2)*((1-gamma)/gamma)^2*theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss))))*oo_.endo_simul(12,1)^2+...
    (bfy/(cy*betaf)*(bfy*(kappa-delta)^2/(betaf*2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)))*oo_.endo_simul(25,1)^2+...
    (nfay/(2*cy)*(1-rhonfa*(1-rhonfa))+nfay^2/(cy*nxy)*((rhonfa^2+1)/2-rhonfa)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy)))*oo_.endo_simul(28,1)^2;

welfare(3)=welfare(2)+(1/(1+iss))^2*((1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,3)^2+...
    (sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,3)^2-...
    nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,3)^2+...
    (bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,3)^2-...
    gcy/(2*cy)*(1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(8,3)^2-...
    giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(26,3)^2+...
    bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,4)^2+...
    (bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,4)^2+...
    ((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,5)^2)+...
    1/(1+iss)*((bfy*(1-delta)/cy*(bfy*(1-delta)/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-alpha*sigma*nss^phi/((css-nss^phi/phi)*(1-alpha))+(delta-1/2))-nss^phi/css*((alpha/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi)-alpha/(2*(1-alpha))))*oo_.endo_simul(6,2)^2-...
    (nss^phi/css*((psi/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-psi/(2*(1-alpha))))*oo_.endo_simul(9,2)^2+...
    (bfy/(2*cy*betaf)+(inty*a2)^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))+inty/(2*cy)*(a2+a2^2)+(inty*a2)^2/(nfay*cy)-nss^phi/css*(eps*(1-alpha+alpha*eps)+sigma*nss^phi/(css-nss^phi/phi)+phi-2)*((1-gamma)/gamma)^2*theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss))))*oo_.endo_simul(12,2)^2+...
    (bfy/(cy*betaf)*(bfy*(kappa-delta)^2/(betaf*2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)))*oo_.endo_simul(25,2)^2+...
    (nfay/(2*cy)*(1-rhonfa*(1-rhonfa))+nfay^2/(cy*nxy)*((rhonfa^2+1)/2-rhonfa)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy)))*oo_.endo_simul(28,2)^2)+...
    (bfy*(1-delta)/(betaf*cy)*(bfy*(1-delta)/(2*betaf*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta-1/2)))*oo_.endo_simul(6,1)^2;
    
    y_sq(3)=(1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,1)^2+...
        1/(1+iss)*((1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,2)^2)+...
        (1/(1+iss))^2*((1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,3)^2);
    rer_sq(3)=(sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,1)^2+...
        1/(1+iss)*((sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,2)^2)+...
        (1/(1+iss))^2*((sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,3)^2);
    infl_sq(3)=-nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,1)^2-...
        1/(1+iss)*(nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,2)^2)-...
        (1/(1+iss))^2*(nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,3)^2);
    inv_sq(3)=(bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,1)^2+...
        1/(1+iss)*((bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,2)^2)+...
        (1/(1+iss))^2*((bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,3)^2);
    gc_sq(3)=-((1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*gcy/(2*cy))*oo_.endo_simul(8,1)^2-...
        1/(1+iss)*(((1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*gcy/(2*cy))*oo_.endo_simul(8,2)^2)-...
        (1/(1+iss))^2*(((1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*gcy/(2*cy))*oo_.endo_simul(8,3)^2);
    gi_sq(3)=-(giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi))))*oo_.endo_simul(26,1)^2-...
        1/(1+iss)*((giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi))))*oo_.endo_simul(26,2)^2)-...
        (1/(1+iss))^2*((giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi))))*oo_.endo_simul(26,3)^2);
    
    rerlead_sq(3)=bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,2)^2+...
        1/(1+iss)*(bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,3)^2)+...
        (1/(1+iss))^2*(bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,4)^2);
    invlead_sq(3)=(bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,2)^2+...
        1/(1+iss)*((bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,3)^2)+...
        (1/(1+iss))^2*((bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,4)^2);
    inv2lead_sq(3)=((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,3)^2+...
        1/(1+iss)*(((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,4)^2)+...
        (1/(1+iss))^2*(((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,5)^2);
    
    knolag_sq(3)=(bfy*(1-delta)/cy*(bfy*(1-delta)/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-alpha*sigma*nss^phi/((css-nss^phi/phi)*(1-alpha))+(delta-1/2))-nss^phi/css*((alpha/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi)-alpha/(2*(1-alpha))))*oo_.endo_simul(6,1)^2+...
        1/(1+iss)*((bfy*(1-delta)/cy*(bfy*(1-delta)/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-alpha*sigma*nss^phi/((css-nss^phi/phi)*(1-alpha))+(delta-1/2))-nss^phi/css*((alpha/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi)-alpha/(2*(1-alpha))))*oo_.endo_simul(6,2)^2);
    kglag_sq(3)=-(nss^phi/css*((psi/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-psi/(2*(1-alpha))))*oo_.endo_simul(9,1)^2-...
        1/(1+iss)*((nss^phi/css*((psi/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-psi/(2*(1-alpha))))*oo_.endo_simul(9,2)^2);
    rerlag_sq(3)=(bfy/(2*cy*betaf)+(inty*a2)^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))+inty/(2*cy)*(a2+a2^2)+(inty*a2)^2/(nfay*cy)-nss^phi/css*(eps*(1-alpha+alpha*eps)+sigma*nss^phi/(css-nss^phi/phi)+phi-2)*((1-gamma)/gamma)^2*theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss))))*oo_.endo_simul(12,1)^2+...
        1/(1+iss)*((bfy/(2*cy*betaf)+(inty*a2)^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))+inty/(2*cy)*(a2+a2^2)+(inty*a2)^2/(nfay*cy)-nss^phi/css*(eps*(1-alpha+alpha*eps)+sigma*nss^phi/(css-nss^phi/phi)+phi-2)*((1-gamma)/gamma)^2*theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss))))*oo_.endo_simul(12,2)^2);
    invlag_sq(3)=(bfy/(cy*betaf)*(bfy*(kappa-delta)^2/(betaf*2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)))*oo_.endo_simul(25,1)^2+...
        1/(1+iss)*((bfy/(cy*betaf)*(bfy*(kappa-delta)^2/(betaf*2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)))*oo_.endo_simul(25,2)^2);
    nfalag_sq(3)=(nfay/(2*cy)*(1-rhonfa*(1-rhonfa))+nfay^2/(cy*nxy)*((rhonfa^2+1)/2-rhonfa)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy)))*oo_.endo_simul(28,1)^2+...
        1/(1+iss)*((nfay/(2*cy)*(1-rhonfa*(1-rhonfa))+nfay^2/(cy*nxy)*((rhonfa^2+1)/2-rhonfa)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy)))*oo_.endo_simul(28,2)^2);
    kno2lag_sq(3)=(bfy*(1-delta)/(betaf*cy)*(bfy*(1-delta)/(2*betaf*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta-1/2)))*oo_.endo_simul(6,1)^2;
    
    for k=4:404
        welfare(k)=welfare(k-1)+(1/(1+iss))^(k-1)*((1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,k)^2+...
            (sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,k)^2-...
            nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,k)^2+...
            (bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,k)^2-...
            gcy/(2*cy)*(1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(8,k)^2-...
            giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi)))*oo_.endo_simul(26,k)^2+...
            bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1)*oo_.endo_simul(12,k+1)^2+...
            (bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,k+1)^2+...
            ((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,k+2)^2)+...
            (1/(1+iss))^(k-2)*((bfy*(1-delta)/cy*(bfy*(1-delta)/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-alpha*sigma*nss^phi/((css-nss^phi/phi)*(1-alpha))+(delta-1/2))-nss^phi/css*((alpha/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi)-alpha/(2*(1-alpha))))*oo_.endo_simul(6,k-1)^2-...
            (nss^phi/css*((psi/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-psi/(2*(1-alpha))))*oo_.endo_simul(9,k-1)^2+...
            (bfy/(2*cy*betaf)+(inty*a2)^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))+inty/(2*cy)*(a2+a2^2)+(inty*a2)^2/(nfay*cy)-nss^phi/css*(eps*(1-alpha+alpha*eps)+sigma*nss^phi/(css-nss^phi/phi)+phi-2)*((1-gamma)/gamma)^2*theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss))))*oo_.endo_simul(12,k-1)^2+...
            (bfy/(cy*betaf)*(bfy*(kappa-delta)^2/(betaf*2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)))*oo_.endo_simul(25,k-1)^2+...
            (nfay/(2*cy)*(1-rhonfa*(1-rhonfa))+nfay^2/(cy*nxy)*((rhonfa^2+1)/2-rhonfa)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy)))*oo_.endo_simul(28,k-1)^2)+...
            (1/(1+iss))^(k-3)*(bfy*(1-delta)/(betaf*cy)*(bfy*(1-delta)/(2*betaf*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta-1/2)))*oo_.endo_simul(6,k-2)^2;

        y_sq(k)=y_sq(k-1)+(1/(1+iss))^(k-1)*((1/(2*cy)*(1-1/cy-sigma*css/((css-nss^phi/phi)*cy))+sigma*nss^phi/((css-nss^phi/phi)*(1-alpha)*noy)*(1/cy-nss^phi/(2*css*(1-alpha)*noy))-nss^phi/(2*css*(1-alpha)*noy)*(1-1/noy+(phi-2)/(noy*(1-alpha))))*oo_.endo_simul(1,k)^2);
        rer_sq(k)=rer_sq(k-1)+(1/(1+iss))^(k-1)*((sigma*nss^phi*oy*fdiy/((css-nss^phi/phi)*(1-alpha)*noy*cy)-bfy/(betaf*cy)-fdiy/(2*cy)*(1+fdiy/cy*(1+sigma*css/(css-nss^phi/phi)))+((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy+inty*(a1+a2))/cy*(sigma*nss^phi*oy/((css-nss^phi/phi)*(1-alpha)*noy)-fdiy/cy*(1+sigma*css*nxy*fdiy/((css-nss^phi/phi)*cy^2)))+(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/(2*cy)*((1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(bfy/betaf-fdiy-(1/betaf-1)*ofy+(1-div)*(1-tauo)*oy)/nxy-1-2*bfy*(1+1/betaf)/nxy*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi))))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*(1/betaf+1)^2+inty^2*(a1+a2)^2/(nfay*cy)-inty/cy*(a1+a2-(a1+a2)^2)/2+(1-nxy/cy-nxy*sigma*css/(cy*(css-nss^phi/phi)))*inty*(a1+a2)/cy*(((1-div)*(1-tauo)*oy-bfy-fdiy-(1/betaf-1)*ofy)/nxy+inty*(a1+a2)/(2*nxy))-nss^phi/css*(1/(2*(1-alpha)^2)*((oy/noy)^2+((1-gamma)/gamma)^2*theta/((1-theta)*(1-theta/(1+iss))))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-oy*(noy+oy)/(2*noy^2*(1-alpha))+eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*((1-gamma)/gamma)^2))*oo_.endo_simul(12,k)^2);
        infl_sq(k)=infl_sq(k-1)-(1/(1+iss))^(k-1)*(nss^phi/css*(eps*theta*(1-alpha+alpha*eps)/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))+theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss)))*(sigma*nss^phi/(css-nss^phi/phi)+phi-2))*oo_.endo_simul(10,k)^2);
        inv_sq(k)=inv_sq(k-1)+(1/(1+iss))^(k-1)*((bfy/cy*(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)-(1-cy-fdiy-gcy-giy-nxy)/(2*cy)*(1+(1-cy-fdiy-gcy-giy-nxy)/cy+sigma*css*(1-cy-fdiy-gcy-giy-nxy)/((css-nss^phi/phi)*cy))-(1-cy-fdiy-gcy-giy-nxy)/cy^2*(1+(1-cy-fdiy-gcy-giy-nxy)*nxy*sigma*css/(cy^2*(css-nss^phi/phi)))*(bfy/betaf*kappa*(1+1/(1+iss))-bfy*(delta-kappa))-bfy*kappa*(1+1/(1+iss))/(betaf*cy)*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa*(1+1/(1+iss))/betaf-(delta-kappa))^2)*oo_.endo_simul(25,k)^2);
        gc_sq(k)=gc_sq(k-1)-(1/(1+iss))^(k-1)*(((1+gcy/cy*(1+sigma*css/(css-nss^phi/phi)))*gcy/(2*cy))*oo_.endo_simul(8,k)^2);
        gi_sq(k)=gi_sq(k-1)-(1/(1+iss))^(k-1)*((giy/(2*cy)*(1+giy/cy*(1+sigma*css/(css-nss^phi/phi))))*oo_.endo_simul(26,k)^2);
        rerlead_sq(k)=rerlead_sq(k-1)+(1/(1+iss))^(k-1)*((bfy/cy*(bfy/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-1))*oo_.endo_simul(12,k+1)^2);
        invlead_sq(k)=invlead_sq(k-1)+(1/(1+iss))^(k-1)*((bfy*kappa/(cy*betaf*(1+iss))*(1/2+kappa/(1+iss))+bfy*kappa*(1+1/(1+iss))/cy*(1/2-kappa*(1+1/(1+iss)))+bfy^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))*(kappa/((1+iss)*betaf)+kappa*(1+1/(1+iss)))^2)*oo_.endo_simul(25,k+1)^2);
        inv2lead_sq(k)=inv2lead_sq(k-1)+(1/(1+iss))^(k-1)*(((bfy*kappa/(1+iss))^2/(2*nxy*cy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-bfy*kappa/(cy*(1+iss))*(1/2+kappa/(1+iss)))*oo_.endo_simul(25,k+2)^2);
        
        knolag_sq(k)=knolag_sq(k-1)+(1/(1+iss))^(k-2)*((bfy*(1-delta)/cy*(bfy*(1-delta)/(2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-alpha*sigma*nss^phi/((css-nss^phi/phi)*(1-alpha))+(delta-1/2))-nss^phi/css*((alpha/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi)-alpha/(2*(1-alpha))))*oo_.endo_simul(6,k-1)^2);
        kglag_sq(k)=kglag_sq(k-1)-(1/(1+iss))^(k-2)*((nss^phi/css*((psi/(1-alpha))^2/2*(sigma*nss^phi/(css-nss^phi/phi)+phi-2)-psi/(2*(1-alpha))))*oo_.endo_simul(9,k-1)^2);
        rerlag_sq(k)=rerlag_sq(k-1)+(1/(1+iss))^(k-2)*((bfy/(2*cy*betaf)+(inty*a2)^2/(2*cy*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))+inty/(2*cy)*(a2+a2^2)+(inty*a2)^2/(nfay*cy)-nss^phi/css*(eps*(1-alpha+alpha*eps)+sigma*nss^phi/(css-nss^phi/phi)+phi-2)*((1-gamma)/gamma)^2*theta/(2*(1-alpha)^2*(1-theta)*(1-theta/(1+iss))))*oo_.endo_simul(12,k-1)^2);
        invlag_sq(k)=invlag_sq(k-1)+(1/(1+iss))^(k-2)*((bfy/(cy*betaf)*(bfy*(kappa-delta)^2/(betaf*2*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta*(1/2-delta)-kappa*(1/2+kappa)+delta*kappa)))*oo_.endo_simul(25,k-1)^2);
        nfalag_sq(k)=nfalag_sq(k-1)+(1/(1+iss))^(k-2)*((nfay/(2*cy)*(1-rhonfa*(1-rhonfa))+nfay^2/(cy*nxy)*((rhonfa^2+1)/2-rhonfa)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy)))*oo_.endo_simul(28,k-1)^2);
        kno2lag_sq(k)=kno2lag_sq(k-1)+(1/(1+iss))^(k-3)*((bfy*(1-delta)/(betaf*cy)*(bfy*(1-delta)/(2*betaf*nxy)*(1-nxy/cy-nxy*sigma*css/((css-nss^phi/phi)*cy))-(delta-1/2)))*oo_.endo_simul(6,k-2)^2);
    end
    
    rertot=rerlead_sq(404)+rerlag_sq(404);
    knotot=knolag_sq(404)+kno2lag_sq(404);
    invtot=invlead_sq(404)+inv2lead_sq(404)+invlag_sq(404);
        
    X(3,i) = welfare(404);
    X(4,i) = y_sq(404);
    X(5,i) = rer_sq(404);
    X(6,i) = infl_sq(404);
    X(7,i) = inv_sq(404);
    X(8,i) = gc_sq(404);
    X(9,i) = gi_sq(404);
    X(10,i) = kglag_sq(404);
    X(11,i) = nfalag_sq(404);
    X(12,i) = knotot;
    X(13,i) = invtot;
    X(14,i) = rertot;    
end
      
[A,I] = max(X(3,:));
disp('Max welfare is '); disp(A);
disp('Corresponding phiinf is '); disp(X(1,I));
disp('Corresponding phier is '); disp(X(2,I));
I