%RAMSEY_SS.M
%Ramsey steady state  in ``Optimal Inflation Stabilization in a Medium-Scale Macroeconomic Model'' by S. Schmitt-Grohe and Martin Uribe (2005)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2005 
function [r_cu, pai_cu, dist, exitflag] = ramsey_ss;

%r_cu = gross quarterly nominal interest rate

%pai_cu = gross quarterly inflation rate

%dist = accuracy of numerical solution to the Ramsey steady state

%exitflag =1 if numerical minimization was successful

allsyms

%Utility function
muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu;
muzstar_fu1 = muupsilon_fu1^(THETA/(1-THETA))*muz_fu1;

objective_cu =  (1-PHI4)*log((c_cu - B * c_ba1/muzstar_cu )) +  PHI4 * log((1-stil_fu1*hd_cu));
objective_fu1 = (1-PHI4)*log((c_fu1 - B * c_cu/muzstar_fu1 )) +  PHI4 * log( (1-stil_fu1p*hd_fu1));

[constraints_cu, constraints_ba1, constraints_fu1] = constraints;

% Lc are the constraints that would be mulitplied by the various xi, the
% Lagrange mulipliers
L = BETTA^(-1) * ( constraints_ba1) + ( constraints_cu) + BETTA * ( constraints_fu1);

endogvar_cu = [u_cu qq_cu iv_cu w_cu hd_cu output_cu c_cu la_cu pai_cu f2_cu x2_cu s_cup stil_cup k_cup r_cu]; 

endogvar_cu=endogvar_cu(:);

endogvar_ba1p = [u_ba1p qq_ba1p iv_ba1p w_ba1p hd_ba1p output_ba1p c_ba1p la_ba1p pai_ba1p f2_ba1p x2_ba1p s_fu1 stil_fu1  k_fu1 r_ba1p]; 

endogvar_ba1p=endogvar_ba1p(:);

dLdx =jacobian(L, endogvar_cu)+jacobian(L, endogvar_ba1p);

O = objective_cu  + BETTA * objective_fu1;

dOdx =jacobian(O, endogvar_cu)+jacobian(O, endogvar_ba1p);

%Initial guess for gross nominal interest rate
X0 = 1.000001;

OPTIONS=optimset('MaxFunEvals',100000,'MaxIter',100000,'TolFun',1e-11,'TolX',1e-11);

[oppolicy,dist,exitflag] = fminsearch(@(x) get_gap(x,dLdx,dOdx),X0,OPTIONS);

[gap,XI] = get_gap(oppolicy,dLdx,dOdx); 

xi1_cu = XI(1);
xi2_cu = XI(2);
xi3_cu = XI(3);
xi4_cu = XI(4);
xi5_cu = XI(5);
xi6_cu = XI(6);
xi7_cu = XI(7);
xi8_cu = XI(8);
xi9_cu = XI(9);
xi10_cu = XI(10);
xi11_cu = XI(11);
xi12_cu = XI(12);
xi13_cu = XI(13);
xi14_cu = XI(14);

xi1_cup = XI(1);
xi2_cup = XI(2);
xi3_cup = XI(3);
xi4_cup = XI(4);
xi5_cup = XI(5);
xi6_cup = XI(6);
xi7_cup = XI(7);
xi8_cup = XI(8);
xi9_cup = XI(9);
xi10_cup = XI(10);
xi11_cup = XI(11);
xi12_cup = XI(12);
xi13_cup = XI(13);
xi14_cup = XI(14);

xi1_ba1 = XI(1);
xi2_ba1 = XI(2);
xi3_ba1 = XI(3);
xi4_ba1 = XI(4);
xi5_ba1 = XI(5);
xi6_ba1 = XI(6);
xi7_ba1 = XI(7);
xi8_ba1 = XI(8);
xi9_ba1 = XI(9);
xi10_ba1 = XI(10);
xi11_ba1 = XI(11);
xi12_ba1 = XI(12);
xi13_ba1 = XI(13);
xi14_ba1 = XI(14);

xi1_ba1p = XI(1);
xi2_ba1p = XI(2);
xi3_ba1p = XI(3);
xi4_ba1p = XI(4);
xi5_ba1p = XI(5);
xi6_ba1p = XI(6);
xi7_ba1p = XI(7);
xi8_ba1p = XI(8);
xi9_ba1p = XI(9);
xi10_ba1p = XI(10);
xi11_ba1p = XI(11);
xi12_ba1p = XI(12);
xi13_ba1p = XI(13);
xi14_ba1p = XI(14);

xi1_fu1 = XI(1);
xi2_fu1 = XI(2);
xi3_fu1 = XI(3);
xi4_fu1 = XI(4);
xi5_fu1 = XI(5);
xi6_fu1 = XI(6);
xi7_fu1 = XI(7);
xi8_fu1 = XI(8);
xi9_fu1 = XI(9);
xi10_fu1 = XI(10);
xi11_fu1 = XI(11);
xi12_fu1 = XI(12);
xi13_fu1 = XI(13);
xi14_fu1 = XI(14);

xi1_fu1p = XI(1);
xi2_fu1p = XI(2);
xi3_fu1p = XI(3);
xi4_fu1p = XI(4);
xi5_fu1p = XI(5);
xi6_fu1p = XI(6);
xi7_fu1p = XI(7);
xi8_fu1p = XI(8);
xi9_fu1p = XI(9);
xi10_fu1p = XI(10);
xi11_fu1p = XI(11);
xi12_fu1p = XI(12);
xi13_fu1p = XI(13);
xi14_fu1p = XI(14);

[ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, u_ba1, v_ba1,  profits_ba1,  qq_ba1, r_ba1, ptil_ba1, mutil_ba1, mc_ba1, s_ba1, htil_ba1, hd_ba1, stil_ba1, k_ba1, h_ba1, iv_ba1, w_ba1, wtil_ba1, output_ba1, g_ba1, mh_ba1,  m_ba1,   c_ba1, la_ba1, pai_ba1, f1_ba1, f2_ba1, x1_ba1, x2_ba1, u_ba1p, v_ba1p, profits_ba1p,  qq_ba1p, r_ba1p, ptil_ba1p, mutil_ba1p, mc_ba1p, s_ba1p, htil_ba1p, hd_ba1p, stil_ba1p, k_ba1p, h_ba1p, iv_ba1p, w_ba1p, wtil_ba1p, output_ba1p, g_ba1p, mh_ba1p,  m_ba1p,  c_ba1p, la_ba1p, pai_ba1p, f1_ba1p, f2_ba1p, x1_ba1p, x2_ba1p, u_ba2, v_ba2, profits_ba2,  qq_ba2, r_ba2, ptil_ba2, mutil_ba2, mc_ba2, s_ba2, htil_ba2, hd_ba2, stil_ba2, k_ba2, h_ba2, iv_ba2, w_ba2, wtil_ba2, output_ba2, g_ba2, mh_ba2,  m_ba2, c_ba2, la_ba2, pai_ba2, f1_ba2, f2_ba2, x1_ba2, x2_ba2, u_ba2p, v_ba2p,  profits_ba2p,  qq_ba2p, r_ba2p, ptil_ba2p, mutil_ba2p, mc_ba2p, s_ba2p, htil_ba2p, hd_ba2p, stil_ba2p, k_ba2p, h_ba2p, iv_ba2p, w_ba2p, wtil_ba2p, output_ba2p, g_ba2p, mh_ba2p, m_ba2p, c_ba2p, la_ba2p, pai_ba2p, f1_ba2p, f2_ba2p, x1_ba2p, x2_ba2p, u_cu, v_cu, profits_cu,  qq_cu, r_cu, ptil_cu, mutil_cu, mc_cu, s_cu, htil_cu, hd_cu, stil_cu, k_cu, h_cu, iv_cu, w_cu, wtil_cu, output_cu, g_cu, mh_cu, m_cu, c_cu, la_cu, pai_cu, f1_cu, f2_cu, x1_cu, x2_cu, u_cup, v_cup, profits_cup,  qq_cup, r_cup, ptil_cup, mutil_cup, mc_cup, s_cup, htil_cup, hd_cup, stil_cup, k_cup, h_cup, iv_cup, w_cup, wtil_cup, output_cup, g_cup, mh_cup, m_cup, c_cup, la_cup, pai_cup, f1_cup, f2_cup, x1_cup, x2_cup, u_fu1, v_fu1, profits_fu1,  qq_fu1, r_fu1, ptil_fu1, mutil_fu1, mc_fu1, s_fu1, htil_fu1, hd_fu1, stil_fu1, k_fu1, h_fu1, iv_fu1, w_fu1, wtil_fu1, output_fu1, g_fu1, mh_fu1, m_fu1, c_fu1, la_fu1, pai_fu1, f1_fu1, f2_fu1, x1_fu1, x2_fu1, u_fu1p, v_fu1p,  profits_fu1p,  qq_fu1p, r_fu1p, ptil_fu1p, mutil_fu1p, mc_fu1p, s_fu1p, htil_fu1p, hd_fu1p, stil_fu1p, k_fu1p, h_fu1p, iv_fu1p, w_fu1p, wtil_fu1p, output_fu1p, g_fu1p, mh_fu1p, m_fu1p, c_fu1p, la_fu1p, pai_fu1p, f1_fu1p, f2_fu1p, x1_fu1p, x2_fu1p, muz_cu, muz_cup, muz_ba1, muz_ba1p, muz_fu1, muz_fu1p,muupsilon_cu, muupsilon_cup, muupsilon_ba1, muupsilon_ba1p, muupsilon_fu1, muupsilon_fu1p,MC, RK,RHOMUZ, RHOMUUPSILON, RHOG, STD_EPSMUZ, STD_EPSMUUPSILON, STD_EPSG, G, MUI, MUZ, MUUPSILON, MUZSTAR, MULAMBDA] = ss(oppolicy);

save ramsey_ss

%SS.M
function [ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, u_ba1, v_ba1, profits_ba1,  qq_ba1, r_ba1, ptil_ba1, mutil_ba1, mc_ba1, s_ba1, htil_ba1, hd_ba1, stil_ba1, k_ba1, h_ba1, iv_ba1, w_ba1, wtil_ba1, output_ba1, g_ba1, mh_ba1,  m_ba1,  c_ba1, la_ba1, pai_ba1, f1_ba1, f2_ba1, x1_ba1, x2_ba1, u_ba1p, v_ba1p, profits_ba1p,  qq_ba1p, r_ba1p, ptil_ba1p, mutil_ba1p, mc_ba1p, s_ba1p, htil_ba1p, hd_ba1p, stil_ba1p, k_ba1p, h_ba1p, iv_ba1p, w_ba1p, wtil_ba1p, output_ba1p, g_ba1p, mh_ba1p,  m_ba1p,  c_ba1p, la_ba1p, pai_ba1p, f1_ba1p, f2_ba1p, x1_ba1p, x2_ba1p, u_ba2, v_ba2, profits_ba2,  qq_ba2, r_ba2, ptil_ba2, mutil_ba2, mc_ba2, s_ba2, htil_ba2, hd_ba2, stil_ba2, k_ba2, h_ba2, iv_ba2, w_ba2, wtil_ba2, output_ba2, g_ba2, mh_ba2,  m_ba2, c_ba2, la_ba2, pai_ba2, f1_ba2, f2_ba2, x1_ba2, x2_ba2, u_ba2p, v_ba2p,  profits_ba2p,  qq_ba2p, r_ba2p, ptil_ba2p, mutil_ba2p, mc_ba2p, s_ba2p, htil_ba2p, hd_ba2p, stil_ba2p, k_ba2p, h_ba2p, iv_ba2p, w_ba2p, wtil_ba2p, output_ba2p, g_ba2p, mh_ba2p, m_ba2p, c_ba2p, la_ba2p, pai_ba2p, f1_ba2p, f2_ba2p, x1_ba2p, x2_ba2p, u_cu, v_cu, profits_cu,  qq_cu, r_cu, ptil_cu, mutil_cu, mc_cu, s_cu, htil_cu, hd_cu, stil_cu, k_cu, h_cu, iv_cu, w_cu, wtil_cu, output_cu, g_cu, mh_cu, m_cu, c_cu, la_cu, pai_cu, f1_cu, f2_cu, x1_cu, x2_cu, u_cup, v_cup, profits_cup,  qq_cup, r_cup, ptil_cup, mutil_cup, mc_cup, s_cup, htil_cup, hd_cup, stil_cup, k_cup, h_cup, iv_cup, w_cup, wtil_cup, output_cup, g_cup, mh_cup, m_cup, c_cup, la_cup, pai_cup, f1_cup, f2_cup, x1_cup, x2_cup, u_fu1, v_fu1, profits_fu1,  qq_fu1, r_fu1, ptil_fu1, mutil_fu1, mc_fu1, s_fu1, htil_fu1, hd_fu1, stil_fu1, k_fu1, h_fu1, iv_fu1, w_fu1, wtil_fu1, output_fu1, g_fu1, mh_fu1, m_fu1, c_fu1, la_fu1, pai_fu1, f1_fu1, f2_fu1, x1_fu1, x2_fu1, u_fu1p, v_fu1p, profits_fu1p,  qq_fu1p, r_fu1p, ptil_fu1p, mutil_fu1p, mc_fu1p, s_fu1p, htil_fu1p, hd_fu1p, stil_fu1p, k_fu1p, h_fu1p, iv_fu1p, w_fu1p, wtil_fu1p, output_fu1p, g_fu1p, mh_fu1p, m_fu1p, c_fu1p, la_fu1p, pai_fu1p, f1_fu1p, f2_fu1p, x1_fu1p, x2_fu1p, muz_cu, muz_cup, muz_ba1, muz_ba1p, muz_fu1, muz_fu1p,muupsilon_cu, muupsilon_cup, muupsilon_ba1, muupsilon_ba1p, muupsilon_fu1, muupsilon_fu1p,MC, RK,RHOMUZ, RHOMUUPSILON, RHOG, STD_EPSMUZ, STD_EPSMUUPSILON, STD_EPSG, G, MUI, MUZ, MUUPSILON, MUZSTAR, MULAMBDA] = ss(R);

R = max(1,R);

[BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, KAPA, CHI, CHITIL, PSSI, PHI1, PHI2, PHI3, PHI4, GAMA1, GAMA2,  RHOMUZ, RHOMUUPSILON, RHOG, STD_EPSMUZ, STD_EPSMUUPSILON, STD_EPSG, G, MUI, MUZ, MUUPSILON] = deep_parameters;

MUZSTAR = MUZ * MUUPSILON^(THETA/(1-THETA));

MULAMBDA = 1/ MUZSTAR;% THIS IS ONLY TRUE IF PHI3=1, ELSE USE  MUZSTAR^((1-PHI3)*(1-PHI4)-1);

PAI = BETTA * R * MULAMBDA; %Inflation rate

QQ = 1; %Tobin's q

V = sqrt(PHI2/PHI1 + 1/PHI1 *(R-1)/R);%Money velocity

ELL = PHI1 * V + PHI2/V - 2 * sqrt(PHI1*PHI2); %Transactions consts

PTIL = ((1-ALFA * PAI^((ETA-1)*(1-CHI))) / (1-ALFA)) ^ (1/(1-ETA)); %relative price of optimizing firms

S = (1-ALFA) * PTIL^(-ETA) / (1-ALFA*PAI^((1-CHI) * ETA)); %distortion between output and total production

WTIL_O_W = ((1-ALFATIL *( MUZSTAR * PAI)^((CHITIL-1)*(1-ETATIL))) / (1-ALFATIL)) ^ (1/(1-ETATIL)); 

%U = 1;

U = sqrt(1+2/GAMA2*(QQ*MUUPSILON/BETTA/MULAMBDA - (1-DELTA)*QQ-GAMA1));

AU = GAMA1*(U-1) + GAMA2/2*(U-1)^2;

RK1 = (QQ*MUUPSILON/BETTA/MULAMBDA-QQ*(1-DELTA) + AU)/U;

RK = GAMA1 + GAMA2 * (U-1);

HTIL_O_HD = WTIL_O_W^(-ETATIL);

STIL = (1-ALFATIL) * WTIL_O_W^(-ETATIL) / (1-ALFATIL * (PAI*MUZSTAR)^((1-CHITIL)*ETATIL));%distortion between hours supplied and hours used as factor input

MC = PTIL * (ETA -1) * (1- ALFA * BETTA * MULAMBDA *  MUZSTAR* PAI^(ETA*(1-CHI))) / ETA / (1-ALFA*BETTA* MULAMBDA * MUZSTAR* PAI^((CHI-1)*(1-ETA))); % marginal cost of producing a unit of final good

H_O_HD = STIL;

K_U_O_HD =  (RK / MC / THETA)^(1/(THETA-1)); %aggregate capital factor input %NOte this is (u*K/(\mu_{I}*h^d))

W = MC * (1-THETA) * (K_U_O_HD)^THETA / (1 + NU*(1-1/R));%average real wage rate

WTIL = WTIL_O_W * W; %Wage charged by optimizing unions


%Prepare coefficients to solve linear equation in HD

b1 = K_U_O_HD^THETA ;

AAA = (1-PHI4) / PHI4  * (ETATIL -1)/ETATIL * (1-B*BETTA * MULAMBDA) / (1-B/MUZSTAR) / (1+2*PHI1*V-2*(PHI1*PHI2)^(1/2)) * WTIL  * (1 - ALFATIL*BETTA*(MUZSTAR*PAI)^((1-CHITIL)*ETATIL)*MULAMBDA*MUZSTAR) /(1 - ALFATIL*BETTA*(MUZSTAR*PAI)^((1-CHITIL)* (ETATIL-1))*MULAMBDA*MUZSTAR);    

b2 = S * ( 1 + ELL) * AAA; 

b3 = S / U *  K_U_O_HD  *	(MUI  -1 + DELTA+AU);

b4 = G * S + PSSI;

HD = (b2+b4) / (b1+STIL*b2-b3); %Composite labor input 

HTIL = HTIL_O_HD * HD; %Labor demanded by reoptimizing unions

H = H_O_HD * HD; %Labor supply

K = K_U_O_HD * HD / U * MUI; %Capital stock

IV = (1-(1-DELTA)/MUI) * K; %investment

%Profits 
PROFITS = ((U*K/MUI)^THETA * HD^(1-THETA) - PSSI)/S - (RK*U*K/MUI + W*HD * (1+NU*(1-1/R)));

OUTPUT = ((U*K/MUI)^THETA * HD^(1-THETA) - PSSI)/S; %Aggregate absorption

C = (OUTPUT-IV-G-AU*K/MUI) / (1+ELL); %Consumption

LA =  (1-B*BETTA*MULAMBDA) * (1-PHI4) / C / (1-B/MUZSTAR) / (1+2*PHI1*V-2*(PHI1*PHI2)^(1/2)); %steady state marginal utility of income

MH = C/V; %Real balances held by households 

M = MH + NU*W*HD; %Total real money holdings

F1 = (ETATIL-1) / ETATIL * WTIL * LA * (W/WTIL)^ETATIL * HD / (1-ALFATIL*BETTA*(MUZSTAR*PAI)^((1-CHITIL)*(ETATIL-1)) *MULAMBDA*MUZSTAR); 

F2 = F1;

X2 = OUTPUT * PTIL^(-ETA) / (1-ALFA* BETTA * PAI^((CHI-1)*(1-ETA)) * MUZSTAR * MULAMBDA);

X1 = (ETA-1) / ETA * X2;

MUTIL = LA * W / PHI4 * (1-H) ; % inverse of wage markup

%MUTIL = (ETATIL / (ETATIL-1)) * ((1-ALFATIL * BETTA * MUZSTAR * MULAMBDA * (MUZSTAR*PAI)^((1-CHITIL)*(ETATIL-1))) / (1 - ALFATIL * BETTA * MUZSTAR * MULAMBDA * (MUZSTAR*PAI)^((1-CHITIL)*ETATIL))) * (1/WTIL_O_W); % inverse of wage markup


muz_cu=MUZ;
muz_cup=MUZ;
muz_ba1=MUZ;
muz_ba1p=MUZ;
muz_fu1=MUZ;
muz_fu1p=MUZ;


muupsilon_cu=MUUPSILON;
u_cu=(U);
v_cu=(V);
profits_cu = PROFITS;
qq_cu = (QQ); 
r_cu = (R);
ptil_cu=(PTIL); 
mutil_cu=(MUTIL); 
mc_cu=(MC); 
s_cu=(S); 
htil_cu=(HTIL);
hd_cu = (HD);
stil_cu=(STIL); 
k_cu=(K); 
h_cu=(H); 
iv_cu=(IV); 
w_cu=(W); 
wtil_cu=(WTIL); 
output_cu=(OUTPUT); 
g_cu=(G); 
mh_cu=(MH);
m_cu = (M);
c_cu=(C); 
la_cu=(LA); 
pai_cu=(PAI); 
f1_cu=(F1);
f2_cu=(F2);
x1_cu=(X1);
x2_cu=(X2);

muupsilon_cup=MUUPSILON;
u_cup=(U);
v_cup=(V);
profits_cup = PROFITS;
qq_cup = (QQ); 
r_cup = (R);
ptil_cup=(PTIL); 
mutil_cup=(MUTIL); 
mc_cup=(MC); 
s_cup=(S); 
htil_cup=(HTIL);
hd_cup = (HD);
stil_cup=(STIL); 
k_cup=(K); 
h_cup=(H); 
iv_cup=(IV); 
w_cup=(W); 
wtil_cup=(WTIL); 
output_cup=(OUTPUT); 
g_cup=(G); 
mh_cup=(MH);
m_cup = (M);
c_cup=(C); 
la_cup=(LA); 
pai_cup=(PAI); 
f1_cup=(F1);
f2_cup=(F2);
x1_cup=(X1);
x2_cup=(X2);

muupsilon_ba1=MUUPSILON;
u_ba1=(U);
v_ba1=(V);
profits_ba1 = PROFITS;
qq_ba1 = (QQ); 
r_ba1 = (R);
ptil_ba1=(PTIL); 
mutil_ba1=(MUTIL); 
mc_ba1=(MC); 
s_ba1=(S); 
htil_ba1=(HTIL);
hd_ba1 = (HD);
stil_ba1=(STIL); 
k_ba1=(K); 
h_ba1=(H); 
iv_ba1=(IV); 
w_ba1=(W); 
wtil_ba1=(WTIL); 
output_ba1=(OUTPUT); 
g_ba1=(G); 
mh_ba1=(MH);
m_ba1 = (M);
c_ba1=(C); 
la_ba1=(LA); 
pai_ba1=(PAI); 
f1_ba1=(F1);
f2_ba1=(F2);
x1_ba1=(X1);
x2_ba1=(X2);

muupsilon_ba1p=MUUPSILON;
u_ba1p=(U);
v_ba1p=(V);
profits_ba1p = PROFITS;
qq_ba1p = (QQ); 
r_ba1p = (R);
ptil_ba1p=(PTIL); 
mutil_ba1p=(MUTIL); 
mc_ba1p=(MC); 
s_ba1p=(S); 
htil_ba1p=(HTIL);
hd_ba1p = (HD);
stil_ba1p=(STIL); 
k_ba1p=(K); 
h_ba1p=(H); 
iv_ba1p=(IV); 
w_ba1p=(W); 
wtil_ba1p=(WTIL); 
output_ba1p=(OUTPUT); 
g_ba1p=(G); 
mh_ba1p=(MH);
m_ba1p = (M);
c_ba1p=(C); 
la_ba1p=(LA); 
pai_ba1p=(PAI); 
f1_ba1p=(F1);
f2_ba1p=(F2);
x1_ba1p=(X1);
x2_ba1p=(X2);

u_ba2=(U);
v_ba2=V;
profits_ba2 = PROFITS;
qq_ba2 = (QQ); 
r_ba2 = (R);
ptil_ba2=(PTIL); 
mutil_ba2=(MUTIL); 
mc_ba2=(MC); 
s_ba2=(S); 
htil_ba2=(HTIL);
hd_ba2 = (HD);
stil_ba2=(STIL); 
k_ba2=(K); 
h_ba2=(H); 
iv_ba2=(IV); 
w_ba2=(W); 
wtil_ba2=(WTIL); 
output_ba2=(OUTPUT); 
g_ba2=(G); 
mh_ba2=(MH);
m_ba2 = (M);
c_ba2=(C); 
la_ba2=(LA); 
pai_ba2=(PAI); 
f1_ba2=(F1);
f2_ba2=(F2);
x1_ba2=(X1);
x2_ba2=(X2);

u_ba2p = (U);
v_ba2p=V;
profits_ba2p = PROFITS;
qq_ba2p = (QQ); 
r_ba2p = (R);
ptil_ba2p=(PTIL); 
mutil_ba2p=(MUTIL); 
mc_ba2p=(MC); 
s_ba2p=(S); 
htil_ba2p=(HTIL);
hd_ba2p = (HD);
stil_ba2p=(STIL); 
k_ba2p=(K); 
h_ba2p=(H); 
iv_ba2p=(IV); 
w_ba2p=(W); 
wtil_ba2p=(WTIL); 
output_ba2p=(OUTPUT); 
g_ba2p=(G); 
mh_ba2p=(MH);
m_ba2p = (M);
c_ba2p=(C); 
la_ba2p=(LA); 
pai_ba2p=(PAI); 
f1_ba2p=(F1);
f2_ba2p=(F2);
x1_ba2p=(X1);
x2_ba2p=(X2);

muupsilon_fu1=MUUPSILON;
v_fu1=(V);
profits_fu1 = PROFITS;
qq_fu1 = (QQ); 
r_fu1 = (R);
ptil_fu1=(PTIL); 
u_fu1=(U); 
mutil_fu1=(MUTIL); 
mc_fu1=(MC); 
s_fu1=(S); 
htil_fu1=(HTIL);
hd_fu1 = (HD);
stil_fu1=(STIL); 
k_fu1=(K); 
h_fu1=(H); 
iv_fu1=(IV); 
w_fu1=(W); 
wtil_fu1=(WTIL); 
output_fu1=(OUTPUT); 
g_fu1=(G); 
mh_fu1=(MH);
m_fu1 = (M);
c_fu1=(C); 
la_fu1=(LA); 
pai_fu1=(PAI); 
f1_fu1=(F1);
f2_fu1=(F2);
x1_fu1=(X1);
x2_fu1=(X2);

muupsilon_fu1p=MUUPSILON;
v_fu1p=(V);
profits_fu1p = PROFITS;
qq_fu1p = (QQ); 
r_fu1p = (R);
ptil_fu1p=(PTIL); 
u_fu1p=(U); 
mutil_fu1p=(MUTIL); 
mc_fu1p=(MC); 
s_fu1p=(S); 
htil_fu1p=(HTIL);
hd_fu1p = (HD);
stil_fu1p=(STIL); 
k_fu1p=(K); 
h_fu1p=(H); 
iv_fu1p=(IV); 
w_fu1p=(W); 
wtil_fu1p=(WTIL); 
output_fu1p=(OUTPUT); 
g_fu1p=(G); 
mh_fu1p=(MH);
m_fu1p = (M);
c_fu1p=(C); 
la_fu1p=(LA); 
pai_fu1p=(PAI); 
f1_fu1p=(F1);
f2_fu1p=(F2);
x1_fu1p=(X1);
x2_fu1p=(X2);


%GET_GAP.M
function [gap,XI] = get_gap(POLGUESS,dLdx,dOdx);

RGUESS = max(1,POLGUESS);


[ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, u_ba1, v_ba1, profits_ba1,  qq_ba1, r_ba1, ptil_ba1, mutil_ba1, mc_ba1, s_ba1, htil_ba1, hd_ba1, stil_ba1, k_ba1, h_ba1, iv_ba1, w_ba1, wtil_ba1, output_ba1, g_ba1, mh_ba1,  m_ba1, c_ba1, la_ba1, pai_ba1, f1_ba1, f2_ba1, x1_ba1, x2_ba1, u_ba1p, v_ba1p, profits_ba1p,  qq_ba1p, r_ba1p, ptil_ba1p, mutil_ba1p, mc_ba1p, s_ba1p, htil_ba1p, hd_ba1p, stil_ba1p, k_ba1p, h_ba1p, iv_ba1p, w_ba1p, wtil_ba1p, output_ba1p, g_ba1p, mh_ba1p,  m_ba1p,  c_ba1p, la_ba1p, pai_ba1p, f1_ba1p, f2_ba1p, x1_ba1p, x2_ba1p, u_ba2, v_ba2, profits_ba2,  qq_ba2, r_ba2, ptil_ba2, mutil_ba2, mc_ba2, s_ba2, htil_ba2, hd_ba2, stil_ba2, k_ba2, h_ba2, iv_ba2, w_ba2, wtil_ba2, output_ba2, g_ba2, mh_ba2,  m_ba2, c_ba2, la_ba2, pai_ba2, f1_ba2, f2_ba2, x1_ba2, x2_ba2, u_ba2p, v_ba2p,  profits_ba2p,  qq_ba2p, r_ba2p, ptil_ba2p, mutil_ba2p, mc_ba2p, s_ba2p, htil_ba2p, hd_ba2p, stil_ba2p, k_ba2p, h_ba2p, iv_ba2p, w_ba2p, wtil_ba2p, output_ba2p, g_ba2p, mh_ba2p, m_ba2p, c_ba2p, la_ba2p, pai_ba2p, f1_ba2p, f2_ba2p, x1_ba2p, x2_ba2p, u_cu, v_cu, profits_cu,  qq_cu, r_cu, ptil_cu, mutil_cu, mc_cu, s_cu, htil_cu, hd_cu, stil_cu, k_cu, h_cu, iv_cu, w_cu, wtil_cu, output_cu, g_cu, mh_cu, m_cu, c_cu, la_cu, pai_cu, f1_cu, f2_cu, x1_cu, x2_cu, u_cup, v_cup, profits_cup,  qq_cup, r_cup, ptil_cup, mutil_cup, mc_cup, s_cup, htil_cup, hd_cup, stil_cup, k_cup, h_cup, iv_cup, w_cup, wtil_cup, output_cup, g_cup, mh_cup, m_cup, c_cup, la_cup, pai_cup, f1_cup, f2_cup, x1_cup, x2_cup, u_fu1, v_fu1, profits_fu1,  qq_fu1, r_fu1, ptil_fu1, mutil_fu1, mc_fu1, s_fu1, htil_fu1, hd_fu1, stil_fu1, k_fu1, h_fu1, iv_fu1, w_fu1, wtil_fu1, output_fu1, g_fu1, mh_fu1, m_fu1, c_fu1, la_fu1, pai_fu1, f1_fu1, f2_fu1, x1_fu1, x2_fu1, u_fu1p, v_fu1p,  profits_fu1p,  qq_fu1p, r_fu1p, ptil_fu1p, mutil_fu1p, mc_fu1p, s_fu1p, htil_fu1p, hd_fu1p, stil_fu1p, k_fu1p, h_fu1p, iv_fu1p, w_fu1p, wtil_fu1p, output_fu1p, g_fu1p, mh_fu1p, m_fu1p, c_fu1p, la_fu1p, pai_fu1p, f1_fu1p, f2_fu1p, x1_fu1p, x2_fu1p, muz_cu, muz_cup, muz_ba1, muz_ba1p, muz_fu1, muz_fu1p,muupsilon_cu, muupsilon_cup, muupsilon_ba1, muupsilon_ba1p, muupsilon_fu1, muupsilon_fu1p,MC, RK,RHOMUZ, RHOMUUPSILON, RHOG, STD_EPSMUZ, STD_EPSMUUPSILON, STD_EPSG, G, MUI, MUZ, MUUPSILON, MUZSTAR, MULAMBDA] =ss(RGUESS);
ndLdx=eval(dLdx);
ndOdx=eval(dOdx);

if r_cu<=1;
     XI = -ndOdx(1:14)*inv(ndLdx(:,1:14));
      gap = (XI*ndLdx+ndOdx);
     sign_right = XI * ndLdx(:,15) + ndOdx(1,15);
     if sign_right<0
        gap = norm(gap(1:14));
     else
        gap = .1;
     end
else
    %Least-Square Method
    XI = -((ndLdx')\ndOdx')';
    residual1 = ndOdx +  XI * ndLdx;
    gap = norm(residual1);
end