function [optimized_rule, dist,exitflag] = competitive_optimized(feature1,feature2)
%[optimized_rule, dist,exitflag] = competitive_optimized(feature1,feature2) finds the coefficients of the interest-rate and fiscal  feedback rules that maximize the unconditional expectation of lifetime utility for the economy with lump-sum taxes and a fiscal rule developed in ``Optimal,  Simple, and Implementable Monetary and Fiscal Rules,'' by S. Schmitt-Grohe and Martin Uribe (2006).. 
%Output: optimized_rule is a vector of optimized coefficients in the order [GAMA1 ALFA_PAI ALFA_Y ALFA_R_BA1]
%dist is the difference between welfare under Ramsey and welfare under the optimized policy rule
%exitflag=1 if the numerical search converged and zero otherwise
%Inputs: feature1 is a string taking the values 'smoothing' (optimization over the lagged value of the interest rate) or 'no_smoothing'  the feedback rule does not feature the lagged interest rate
%feature2 is a string taking the values 'crule', 'frule', 'brule', or 'drule'. The default is  'crule'
%If only feature1 is provided then  feature2='crule'
%If no inputs are provided then feature1='no_smoothing' and feature2='crule'
%(c) Stephanie Schmitt-Grohe and Martin Uribe, June 2006

if nargin<2
feature2 = 'crule'
end

if nargin<1
feature1 = 'no_smoothing'
end


load competitive_f %produced by COMPETITIVE_F.M

load ramsey_ss %This is saved automatically, when running ramsey_ss.m

%Assign policy rule parameters simply to allow for an initial numerical evaluation of the model
crule = 0;
brule = 0;
frule = 0;
drule = 0;
ALFA_PAI = 1.5;
ALFA_Y = 0;
ALFA_R_BA1 = 0;
ALFA_R_CU = 1;
GAMA1 = 0.5;

%Order of approximation
approx = 2;

%Recall that the innovation to the state vector is given by sigma * ETAMATRIX * eps_{t+1}, where sigma is a scalar, ETAMATRIX is a nx*ne matrix and eps_t is a vector of exogenous innovations of size neX1 distributed (0,I);
ne=2; 
ETAMATRIX = zeros(size(fx,2),ne); 

%The gov't purchases and techno shocks take the last 2 elements of the state vector x (see ramsey_f.M). Thus,
ETAMATRIX(end-1,1)=STD_EPSG;
ETAMATRIX(end,2)=STD_EPSZ;

if strcmp(feature2,'crule') 
crule = 1;
elseif strcmp(feature2,'brule')
brule = 1;
elseif strcmp(feature2,'frule')
frule =1;
elseif strcmp(feature2,'drule')
drule =1;
end

%Initial guess for coefficients of policy rule
if strcmp(feature1,'no_smoothing') 
POL0 = [0.1 3 0];
elseif strcmp(feature1,'smoothing')
POL0 = [0.2 3 0 0];
end

%Numerical Evaluation of Equilibrium Conditions
num_eval

nr = find(controlvar_cu=='r_cu');

save alles.mat RSTAR PAISTAR nr ALFA_R_CU ALFA_R_BA1   ETAMATRIX brule crule frule drule vt_cu nfx nfxp nfy nfyp nfypyp nfypy nfypxp nfypx nfyyp nfyy nfyxp nfyx nfxpyp nfxpy nfxpxp nfxpx nfxyp nfxy nfxxp nfxx

OPTIONS=optimset('MaxFunEvals',10000,'MaxIter',10000,'TolFun',1e-10,'TolX',1e-10);

[optimized_rule, dist,exitflag] = fminsearch(@(x) get_gap(x), POL0, OPTIONS);

GAMA1 = optimized_rule(1);
ALFA_PAI = optimized_rule(2);
ALFA_Y = optimized_rule(3);
ALFA_R_BA1 = (length(optimized_rule)==4) * optimized_rule(end);

eval(['save optimized_' feature1 '_' feature2 '.mat'])

%GET_GAP.M
%This subfunctions computes a welfare-based difference between the competitive equilibrium and the ramsey equilibrium for a given monetary/fiscal policy

function  welfare_gap = get_gap(x);

if max(x)>3 | min(x)<0

welfare_gap = 1000;

else

load alles.mat


%Reevaluate the relevant equations.
GAMA1 = x(1);
ALFA_PAI = x(2);
ALFA_Y = x(3);

if length(x) == 4
ALFA_R_BA1 =  x(end);
end

nfx(end-1:end,:) = [[                                                                               0, crule*ALFA_R_BA1+brule*ALFA_R_BA1+frule*ALFA_R_BA1+drule*(ALFA_R_CU+ALFA_R_BA1),                                                               -drule*ALFA_R_BA1,                                                   brule*ALFA_PAI+drule*ALFA_PAI,                                                       brule*ALFA_Y+drule*ALFA_Y,                                                                   -drule*ALFA_Y,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0];
[                                                                           GAMA1,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0,                                                                               0]];
 
 
nfy(end-1:end,:) = [[                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                     crule*ALFA_Y,                                                   crule*ALFA_PAI,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0, -crule*ALFA_R_CU-brule*ALFA_R_CU-frule*ALFA_R_CU-drule*ALFA_R_CU];
[                                                                0,                                                                0,                                                               -1,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0,                                                                0]];
 
 
nfyp(end-1:end,:) = [[              0,              0,              0,              0,              0,              0,              0,              0,              0,   frule*ALFA_Y, frule*ALFA_PAI,              0,              0,              0,              0,              0,              0,              0];
[              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0,              0]];
 
 
[gx,hx,exitflag] = gx_hx(nfy,nfx,nfyp,nfxp);

if exitflag==1

%Compute standard deviation of interest rate
VARY = mom(gx,hx,ETAMATRIX*ETAMATRIX',0);
std_r_cu = VARY(nr,nr)^(1/2);

%Check that interest rate is away from zero

if log(RSTAR) > 2* std_r_cu

[gxx,hxx] = gxx_hxx_noloop(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx);

[gss,hss] = gss_hss(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx,gxx,ETAMATRIX);


[Ey,Ex] = unconditional_mean(gx, hx, gxx, hxx, gss, hss, ETAMATRIX, 1);

unconditional_welfare_competitive = Ey(1) + vt_cu;

welfare_gap = -unconditional_welfare_competitive

x

else %if log(RSTAR)>2*std_r_cu

 welfare_gap =1000; 

end %if log(RSTAR)>2*std_r_cu


else %if exitflag==1

welfare_gap = 1000;

end %if exitflag==1

end %if max(abs(x))>3