%RAMSEY_SS.M
%Ramsey steady state under income taxation in ``Optimal Fiscal and Monetary Policy in a Medium-Scale Macroeconomic Model'' by S. Schmitt-Grohe and Martin Uribe (2005)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, April 2005 
function [ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, TAUPROFITS, s_cu, stil_cu, assets_cu, k_cu, n_cu, g_cu, z_cu, qq_cu, hd_cu, output_cu, c_cu, la_cu, pai_cu, u_cu, f2_cu, x2_cu, r_cu, tauh_cu, tauk_cu, iv_cu, w_cu, xi1_cu, xi2_cu, xi3_cu, xi4_cu, xi5_cu, xi6_cu, xi7_cu, xi8_cu, xi9_cu, xi10_cu, xi11_cu, xi12_cu, xi13_cu, xi14_cu, xi15_cu, s_cup, stil_cup, assets_cup, k_cup, n_cup, g_cup, z_cup, qq_cup, hd_cup, output_cup, c_cup, la_cup, pai_cup, u_cup, f2_cup, x2_cup, r_cup, tauh_cup, tauk_cup, iv_cup, w_cup, xi1_cup, xi2_cup, xi3_cup, xi4_cup, xi5_cup, xi6_cup, xi7_cup, xi8_cup, xi9_cup, xi10_cup, xi11_cup, xi12_cup, xi13_cup, xi14_cup, xi15_cup, s_ba1, stil_ba1, assets_ba1, k_ba1, n_ba1, g_ba1, z_ba1, qq_ba1, hd_ba1, output_ba1, c_ba1, la_ba1, pai_ba1, u_ba1, f2_ba1, x2_ba1, r_ba1, tauh_ba1, tauk_ba1, iv_ba1, w_ba1, xi1_ba1, xi2_ba1, xi3_ba1, xi4_ba1, xi5_ba1, xi6_ba1, xi7_ba1, xi8_ba1, xi9_ba1, xi10_ba1, xi11_ba1, xi12_ba1, xi13_ba1, xi14_ba1, xi15_ba1, s_ba1p, stil_ba1p, assets_ba1p, k_ba1p, n_ba1p, g_ba1p, z_ba1p, qq_ba1p, hd_ba1p, output_ba1p, c_ba1p, la_ba1p, pai_ba1p, u_ba1p, f2_ba1p, x2_ba1p, r_ba1p, tauh_ba1p, tauk_ba1p, iv_ba1p, w_ba1p, xi1_ba1p, xi2_ba1p, xi3_ba1p, xi4_ba1p, xi5_ba1p, xi6_ba1p, xi7_ba1p, xi8_ba1p, xi9_ba1p, xi10_ba1p, xi11_ba1p, xi12_ba1p, xi13_ba1p, xi14_ba1p, xi15_ba1p, s_fu1, stil_fu1, assets_fu1, k_fu1, n_fu1, g_fu1, z_fu1, qq_fu1, hd_fu1, output_fu1, c_fu1, la_fu1, pai_fu1, u_fu1, f2_fu1, x2_fu1, r_fu1, tauh_fu1, tauk_fu1, iv_fu1, w_fu1, xi1_fu1, xi2_fu1, xi3_fu1, xi4_fu1, xi5_fu1, xi6_fu1, xi7_fu1, xi8_fu1, xi9_fu1, xi10_fu1, xi11_fu1, xi12_fu1, xi13_fu1, xi14_fu1, xi15_fu1, s_fu1p, stil_fu1p, assets_fu1p, k_fu1p, n_fu1p, g_fu1p, z_fu1p, qq_fu1p, hd_fu1p, output_fu1p, c_fu1p, la_fu1p, pai_fu1p, u_fu1p, f2_fu1p, x2_fu1p, r_fu1p, tauh_fu1p, tauk_fu1p, iv_fu1p, w_fu1p, xi1_fu1p, xi2_fu1p, xi3_fu1p, xi4_fu1p, xi5_fu1p, xi6_fu1p, xi7_fu1p, xi8_fu1p, xi9_fu1p, xi10_fu1p, xi11_fu1p, xi12_fu1p, xi13_fu1p, xi14_fu1p, xi15_fu1p, s_ba2, stil_ba2, assets_ba2, k_ba2, n_ba2, g_ba2, z_ba2, qq_ba2, hd_ba2, output_ba2, c_ba2, la_ba2, pai_ba2, u_ba2, f2_ba2, x2_ba2, r_ba2, tauh_ba2, tauk_ba2, iv_ba2, w_ba2, xi1_ba2, xi2_ba2, xi3_ba2, xi4_ba2, xi5_ba2, xi6_ba2, xi7_ba2, xi8_ba2, xi9_ba2, xi10_ba2, xi11_ba2, xi12_ba2, xi13_ba2, xi14_ba2, xi15_ba2, s_ba2p, stil_ba2p, assets_ba2p, k_ba2p, n_ba2p, g_ba2p, z_ba2p, qq_ba2p, hd_ba2p, output_ba2p, c_ba2p, la_ba2p, pai_ba2p, u_ba2p, f2_ba2p, x2_ba2p, r_ba2p, tauh_ba2p, tauk_ba2p, iv_ba2p, w_ba2p, xi1_ba2p, xi2_ba2p, xi3_ba2p, xi4_ba2p, xi5_ba2p, xi6_ba2p, xi7_ba2p, xi8_ba2p, xi9_ba2p, xi10_ba2p, xi11_ba2p, xi12_ba2p, xi13_ba2p, xi14_ba2p, xi15_ba2p, MC, RK, N, G, RHOG, RHOZ, RHON, STD_EPSG, STD_EPSZ, STD_EPSN, dist, exitflag] = ramsey_ss;

allsyms

[objective_cu, constraints_cu, objective_ba1, constraints_ba1, objective_fu1, 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);


%Note: We aren't taking foc w.r.t. assets because it is redundant

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 tauk_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 tauk_ba1p]; 

endogvar_ba1p=endogvar_ba1p(:);

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

O = BETTA^(-1) * (objective_ba1) + (objective_cu) + BETTA * (objective_fu1);

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


%Initial guess for tax rate and gross nominal interest rate
X0 = [0.3; 1.01];

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);
xi15_cu = XI(15);

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);
xi15_cup = XI(15);

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);
xi15_ba1 = XI(15);

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);
xi15_ba1p = XI(15);

xi1_ba2 = XI(1);
xi2_ba2 = XI(2);
xi3_ba2 = XI(3);
xi4_ba2 = XI(4);
xi5_ba2 = XI(5);
xi6_ba2 = XI(6);
xi7_ba2 = XI(7);
xi8_ba2 = XI(8);
xi9_ba2 = XI(9);
xi10_ba2 = XI(10);
xi11_ba2 = XI(11);
xi12_ba2 = XI(12);
xi13_ba2 = XI(13);
xi14_ba2 = XI(14);
xi15_ba2 = XI(15);

xi1_ba2p = XI(1);
xi2_ba2p = XI(2);
xi3_ba2p = XI(3);
xi4_ba2p = XI(4);
xi5_ba2p = XI(5);
xi6_ba2p = XI(6);
xi7_ba2p = XI(7);
xi8_ba2p = XI(8);
xi9_ba2p = XI(9);
xi10_ba2p = XI(10);
xi11_ba2p = XI(11);
xi12_ba2p = XI(12);
xi13_ba2p = XI(13);
xi14_ba2p = XI(14);
xi15_ba2p = XI(15);

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);
xi15_fu1 = XI(15);

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);
xi15_fu1p = XI(15);

[ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, TAUPROFITS, u_ba1, v_ba1, tauh_ba1, tauk_ba1, profits_ba1, z_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, n_ba1, m_ba1, assets_ba1, tau_ba1, c_ba1, la_ba1, pai_ba1, f1_ba1, f2_ba1, x1_ba1, x2_ba1, u_ba1p, v_ba1p, tauh_ba1p, tauk_ba1p, profits_ba1p, z_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, n_ba1p, m_ba1p, assets_ba1p, tau_ba1p, c_ba1p, la_ba1p, pai_ba1p, f1_ba1p, f2_ba1p, x1_ba1p, x2_ba1p, u_ba2, v_ba2, tauh_ba2, tauk_ba2, profits_ba2, z_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, n_ba2, m_ba2, assets_ba2, tau_ba2, c_ba2, la_ba2, pai_ba2, f1_ba2, f2_ba2, x1_ba2, x2_ba2, u_ba2p, v_ba2p, tauh_ba2p, tauk_ba2p, profits_ba2p, z_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, n_ba2p, m_ba2p, assets_ba2p, tau_ba2p, c_ba2p, la_ba2p, pai_ba2p, f1_ba2p, f2_ba2p, x1_ba2p, x2_ba2p, u_cu, v_cu, tauh_cu, tauk_cu, profits_cu, z_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, n_cu, m_cu, assets_cu, tau_cu, c_cu, la_cu, pai_cu, f1_cu, f2_cu, x1_cu, x2_cu, u_cup, v_cup, tauh_cup, tauk_cup, profits_cup, z_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, n_cup, m_cup, assets_cup, tau_cup, c_cup, la_cup, pai_cup, f1_cup, f2_cup, x1_cup, x2_cup, u_fu1, v_fu1, tauh_fu1, tauk_fu1, profits_fu1, z_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, n_fu1, m_fu1, assets_fu1, tau_fu1, c_fu1, la_fu1, pai_fu1, f1_fu1, f2_fu1, x1_fu1, x2_fu1, u_fu1p, v_fu1p, tauh_fu1p, tauk_fu1p, profits_fu1p, z_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, n_fu1p, m_fu1p, assets_fu1p, tau_fu1p, c_fu1p, la_fu1p, pai_fu1p, f1_fu1p, f2_fu1p, x1_fu1p, x2_fu1p, MC, RK, RHOZ, RHOG, RHON, N, G, STD_EPSG, STD_EPSZ, STD_EPSN]=ss(oppolicy(1), oppolicy(2));

%SS.M
%Steady-state of the model in ``Stabilization Policy for the United States,'' by S. Schmitt-Grohe and Martin Uribe   
function [ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, TAUPROFITS, u_ba1, v_ba1, tauh_ba1, tauk_ba1, profits_ba1, z_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, n_ba1, m_ba1, assets_ba1, tau_ba1, c_ba1, la_ba1, pai_ba1, f1_ba1, f2_ba1, x1_ba1, x2_ba1, u_ba1p, v_ba1p, tauh_ba1p, tauk_ba1p, profits_ba1p, z_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, n_ba1p, m_ba1p, assets_ba1p, tau_ba1p, c_ba1p, la_ba1p, pai_ba1p, f1_ba1p, f2_ba1p, x1_ba1p, x2_ba1p, u_ba2, v_ba2, tauh_ba2, tauk_ba2, profits_ba2, z_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, n_ba2, m_ba2, assets_ba2, tau_ba2, c_ba2, la_ba2, pai_ba2, f1_ba2, f2_ba2, x1_ba2, x2_ba2, u_ba2p, v_ba2p, tauh_ba2p, tauk_ba2p, profits_ba2p, z_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, n_ba2p, m_ba2p, assets_ba2p, tau_ba2p, c_ba2p, la_ba2p, pai_ba2p, f1_ba2p, f2_ba2p, x1_ba2p, x2_ba2p, u_cu, v_cu, tauh_cu, tauk_cu, profits_cu, z_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, n_cu, m_cu, assets_cu, tau_cu, c_cu, la_cu, pai_cu, f1_cu, f2_cu, x1_cu, x2_cu, u_cup, v_cup, tauh_cup, tauk_cup, profits_cup, z_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, n_cup, m_cup, assets_cup, tau_cup, c_cup, la_cup, pai_cup, f1_cup, f2_cup, x1_cup, x2_cup, u_fu1, v_fu1, tauh_fu1, tauk_fu1, profits_fu1, z_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, n_fu1, m_fu1, assets_fu1, tau_fu1, c_fu1, la_fu1, pai_fu1, f1_fu1, f2_fu1, x1_fu1, x2_fu1, u_fu1p, v_fu1p, tauh_fu1p, tauk_fu1p, profits_fu1p, z_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, n_fu1p, m_fu1p, assets_fu1p, tau_fu1p, c_fu1p, la_fu1p, pai_fu1p, f1_fu1p, f2_fu1p, x1_fu1p, x2_fu1p, MC, RK, RHOZ, RHOG, RHON, N, G, STD_EPSG, STD_EPSZ, STD_EPSN, DEBT_OVER_GDP]=ss(TAUK, R);

%This program produces the structural parameters and the Ramsey steady-state values of endogenous variables of a standard rbc model augmented to allow for price and wage  stickiness and fiscal and monetary policy  
%(c) Stephanie Schmitt-Grohe and Martin Uribe, December 17, 2004. 

[BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, KAPA, CHI, CHITIL, PSSI, PHI1, PHI2, PHI3, PHI4, GAMA1, GAMA2,  RHOZ, RHOG, RHON, STD_EPSZ, STD_EPSG, STD_EPSN, Z, G, N, DEBT_OVER_GDP] = deep_parameters;

if nargin==0
TAUK= (0.407 + 0.285)/2;
R = 1.042^(1/4) * 1.04^(1/4)
end

TAUH = TAUK;

%Profit taxes
TAUPROFITS = TAUK; %-0.29%1;%1;%

PAI = BETTA * R; %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 * PAI^((CHITIL-1)*(1-ETATIL))) / (1-ALFATIL)) ^ (1/(1-ETATIL)); 

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

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

RK = (QQ * (1-BETTA * (1-DELTA + DELTA * TAUK)) / U/ BETTA / (1-TAUK))+AU/U;

%GAMA1 = RK;

%GAMA2 = GAMA2_O_GAMA1 * GAMA1;

X1_O_X2 = (ETA-1) / ETA;

HTIL_O_HD = WTIL_O_W^(-ETATIL);

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

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

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

K_U_O_HD = (RK / MC / THETA)^(1/(THETA-1)); %aggregate capital factor input

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

a1 = (K_U_O_HD^THETA /S);
a2 = W * (1-TAUH) * (1-PHI4) / (PHI4 * MUTIL) * (1-B*BETTA) * (1+ELL) / (1-B) / (1+2*PHI1*V-2*(PHI1*PHI2)^(1/2));    

a3 = K_U_O_HD / U *	(DELTA + AU);

a4 = G + PSSI/S;

HD = (a2+a4) / (a1+STIL*a2-a3); %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; %Capital stock

IV = DELTA * K; %investment

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

TAU = TAUK * (RK*U - AU -QQ*DELTA)*K + TAUH * W *HD + TAUPROFITS * PROFITS; %Total taxes

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

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

LA =  (1-B*BETTA) * (1-PHI4) / C / (1-B) / (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

ASSETS = (G+N-TAU-M*(1-1/R)) / (1/R-1/PAI);

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

F2 = F1;

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

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

u_cu=(U);
v_cu=(V);
tauh_cu = TAUH;
tauk_cu = TAUK;
profits_cu = PROFITS;
z_cu=(Z); 
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);
n_cu = (N);
m_cu = (M);
assets_cu = ASSETS; 
tau_cu = (TAU);
c_cu=(C); 
la_cu=(LA); 
pai_cu=(PAI); 
f1_cu=(F1);
f2_cu=(F2);
x1_cu=(X1);
x2_cu=(X2);

u_cup=(U);
v_cup=(V);
tauh_cup = TAUH;
tauk_cup = TAUK;
profits_cup = PROFITS;
z_cup=(Z); 
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);
n_cup = (N);
m_cup = (M);
assets_cup = ASSETS; 
tau_cup = (TAU);
c_cup=(C); 
la_cup=(LA); 
pai_cup=(PAI); 
f1_cup=(F1);
f2_cup=(F2);
x1_cup=(X1);
x2_cup=(X2);

u_ba1=(U);
v_ba1=(V);
tauh_ba1 = TAUH;
tauk_ba1 = TAUK;
profits_ba1 = PROFITS;
z_ba1=(Z); 
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);
n_ba1 = (N);
m_ba1 = (M);
assets_ba1 = ASSETS; 
tau_ba1 = (TAU);
c_ba1=(C); 
la_ba1=(LA); 
pai_ba1=(PAI); 
f1_ba1=(F1);
f2_ba1=(F2);
x1_ba1=(X1);
x2_ba1=(X2);

u_ba1p=(U);
v_ba1p=(V);
tauh_ba1p = TAUH;
tauk_ba1p = TAUK;
profits_ba1p = PROFITS;
z_ba1p=(Z); 
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);
n_ba1p = (N);
m_ba1p = (M);
assets_ba1p = ASSETS; 
tau_ba1p = (TAU);
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;
tauh_ba2 = TAUH;
tauk_ba2 = TAUK;
profits_ba2 = PROFITS;
z_ba2=(Z); 
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);
n_ba2 = (N);
m_ba2 = (M);
assets_ba2 = ASSETS; 
tau_ba2 = (TAU);
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;
tauh_ba2p = TAUH;
tauk_ba2p = TAUK;
profits_ba2p = PROFITS;
z_ba2p=(Z); 
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);
n_ba2p = (N);
m_ba2p = (M);
assets_ba2p = ASSETS; 
tau_ba2p = (TAU);
c_ba2p=(C); 
la_ba2p=(LA); 
pai_ba2p=(PAI); 
f1_ba2p=(F1);
f2_ba2p=(F2);
x1_ba2p=(X1);
x2_ba2p=(X2);

v_fu1=(V);
tauh_fu1 = TAUH;
tauk_fu1 = TAUK;
profits_fu1 = PROFITS;
z_fu1=(Z); 
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);
n_fu1 = (N);
m_fu1 = (M);
assets_fu1 = ASSETS; 
tau_fu1 = (TAU);
c_fu1=(C); 
la_fu1=(LA); 
pai_fu1=(PAI); 
f1_fu1=(F1);
f2_fu1=(F2);
x1_fu1=(X1);
x2_fu1=(X2);


v_fu1p=(V);
tauh_fu1p = TAUH;
tauk_fu1p = TAUK;
profits_fu1p = PROFITS;
z_fu1p=(Z); 
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);
n_fu1p = (N);
m_fu1p = (M);
assets_fu1p = ASSETS; 
tau_fu1p = (TAU);
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);

TAUKGUESS= POLGUESS(1);
RGUESS = POLGUESS(2);

TAUHGUESS = TAUKGUESS;

[ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHI3, PHI4, PSSI, THETA, GAMA1, GAMA2, TAUPROFITS, u_ba1, v_ba1, tauh_ba1, tauk_ba1, profits_ba1, z_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, n_ba1, m_ba1, assets_ba1, tau_ba1, c_ba1, la_ba1, pai_ba1, f1_ba1, f2_ba1, x1_ba1, x2_ba1, u_ba1p, v_ba1p, tauh_ba1p, tauk_ba1p, profits_ba1p, z_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, n_ba1p, m_ba1p, assets_ba1p, tau_ba1p, c_ba1p, la_ba1p, pai_ba1p, f1_ba1p, f2_ba1p, x1_ba1p, x2_ba1p, u_ba2, v_ba2, tauh_ba2, tauk_ba2, profits_ba2, z_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, n_ba2, m_ba2, assets_ba2, tau_ba2, c_ba2, la_ba2, pai_ba2, f1_ba2, f2_ba2, x1_ba2, x2_ba2, u_ba2p, v_ba2p, tauh_ba2p, tauk_ba2p, profits_ba2p, z_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, n_ba2p, m_ba2p, assets_ba2p, tau_ba2p, c_ba2p, la_ba2p, pai_ba2p, f1_ba2p, f2_ba2p, x1_ba2p, x2_ba2p, u_cu, v_cu, tauh_cu, tauk_cu, profits_cu, z_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, n_cu, m_cu, assets_cu, tau_cu, c_cu, la_cu, pai_cu, f1_cu, f2_cu, x1_cu, x2_cu, u_cup, v_cup, tauh_cup, tauk_cup, profits_cup, z_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, n_cup, m_cup, assets_cup, tau_cup, c_cup, la_cup, pai_cup, f1_cup, f2_cup, x1_cup, x2_cup, u_fu1, v_fu1, tauh_fu1, tauk_fu1, profits_fu1, z_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, n_fu1, m_fu1, assets_fu1, tau_fu1, c_fu1, la_fu1, pai_fu1, f1_fu1, f2_fu1, x1_fu1, x2_fu1, u_fu1p, v_fu1p, tauh_fu1p, tauk_fu1p, profits_fu1p, z_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, n_fu1p, m_fu1p, assets_fu1p, tau_fu1p, c_fu1p, la_fu1p, pai_fu1p, f1_fu1p, f2_fu1p, x1_fu1p, x2_fu1p, MC, RK, RHOZ, RHOG, RHON, N, G, STD_EPSG, STD_EPSZ, STD_EPSN,DEBT_OVER_GDP]=ss(TAUKGUESS, RGUESS);

ndLdx=eval(dLdx);
ndOdx=eval(dOdx);

%Least-Square Method

XI = -((ndLdx')\ndOdx')';

residual1 = ndOdx +  XI * ndLdx;

residual2 = (assets_cu-m_cu)/output_cu/r_cu-DEBT_OVER_GDP; 

residual3 = 100 * (1-r_cu) * (r_cu<1);

gap = norm(residual1) + abs(residual2) + residual3;