%RAMSEY_SS.M
%ramsey_ss(model,X0) computes 
%the exact nonlinar solution to the 
%deterministic equilibrium conditions of the Ramsey problem in the model economy developed in 
%``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by Stephanie Schmitt-Grohe and 
%Martin Uribe (2006).
%Inputs: 
%model = a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%The default value is 'sgu'.
%X0 = initial guess for the Ramsey steady-state value of the gross quarterly nominal interest 
%rate. The default value is 1.0109;
%Ouput is saved in a mat file called ramsey_ss_ehl.mat  or ramsey_ss_sgu.mat depending on the value of the variable model
%dist = a measure of numerical error in the solution.
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.

function  ramsey_ss(model,X0);
if nargin==0
    model='sgu';
    X0 = 1.0109;
elseif nargin==1
    X0 = 1.0109;
end

%define symbolic variables
allsyms

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

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


%Competitive equilibrium conditions
[constraints_cu, constraints_ba1, constraints_fu1] = constraints;

%EHL model
if strcmp(model,'ehl')
constraints_cu=constraints_cu(1:end-2);
constraints_ba1=constraints_ba1(1:end-2);
constraints_fu1=constraints_fu1(1:end-2);
else
%SGU model
constraints_cu=constraints_cu([1:end-4 end-1:end]);
constraints_ba1=constraints_ba1([1:end-4 end-1:end]);
constraints_fu1=constraints_fu1([1:end-4 end-1:end]);
end

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


%Endogenous variables
endogvar_cu = [stil_cup ptil_cu wtil_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  k_cup r_cu ]; 

endogvar_cu=endogvar_cu(:);

endogvar_ba1p = [stil_fu1 ptil_ba1p wtil_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 k_fu1 r_ba1p]; 


endogvar_ba1p=endogvar_ba1p(:);


%First-order conditions of Ramsey problem

%Derivative of the constraints
dLdx =jacobian(L, endogvar_cu)+jacobian(L, endogvar_ba1p);

%Objective function
O = objective_cu  + BETTA * objective_fu1;

%Derivative of the objective function
dOdx =jacobian(O, endogvar_cu)+jacobian(O, endogvar_ba1p);



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

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


%Compute numerical error of solution, dist, and steady-state values of Ramsey Lagrange 
%multipliers, XI
[dist,XI] = get_gap(oppolicy,dLdx,dOdx, model); 

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_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);


    xi16_cu = XI(16);
    xi16_cup = XI(16);
    xi16_ba1 = XI(16);
    xi16_ba1p = XI(16);
    xi16_fu1 = XI(16);
    xi16_fu1p = XI(16);


%Steady-state values
[ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHIL, PHI5, PSSI, THETA, GAMA1, GAMA2, u_ba1, v_ba1, profits_ba1,  qq_ba1, r_ba1, ptil_ba1, mc_ba1, s_ba1, hd_ba1, stil_ba1, k_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, mc_ba1p, s_ba1p, hd_ba1p, stil_ba1p, k_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,  mc_ba2, s_ba2, hd_ba2, stil_ba2, k_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, mc_ba2p, s_ba2p, hd_ba2p, stil_ba2p, k_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,  mc_cu, s_cu, hd_cu, stil_cu, k_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, mc_cup, s_cup, hd_cup, stil_cup, k_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, mc_fu1, s_fu1,  hd_fu1, stil_fu1, k_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, mc_fu1p, s_fu1p, hd_fu1p, stil_fu1p, k_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, model);


%Growth rate of output
muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu;


muzstar_ba1 = muupsilon_ba1^(THETA/(1-THETA))*muz_ba1;

%Utility function
uf_cu = log(c_cu - B * c_ba1 / muzstar_cu)- PHIL * (stil_cup*hd_cu)^(1+PHI5)/(1+PHI5);

vt_cu =   uf_cu/(1- BETTA);
vt_cup = vt_cu;

dist

%save output
eval(['save ramsey_ss_' model])

%SS.M
%Steady-state of the model in %``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by Stephanie Schmitt-Grohe and 
%Martin Uribe (2006).
%This program produces the structural parameters and the associated Ramsey steady-state values for a given value of the gross quarterly nominal interest rate.   
%Inputs:
%R = gross quarterly nominal interest 
%rate.
%model = a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.
  
function [ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHIL, PHI5, PSSI, THETA, GAMA1, GAMA2, u_ba1, v_ba1, profits_ba1,  qq_ba1, r_ba1, ptil_ba1, mc_ba1, s_ba1, hd_ba1, stil_ba1, k_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, mc_ba1p, s_ba1p, hd_ba1p, stil_ba1p, k_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,  mc_ba2, s_ba2, hd_ba2, stil_ba2, k_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, mc_ba2p, s_ba2p, hd_ba2p, stil_ba2p, k_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,  mc_cu, s_cu, hd_cu, stil_cu, k_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, mc_cup, s_cup, hd_cup, stil_cup, k_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, mc_fu1, s_fu1,  hd_fu1, stil_fu1, k_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, mc_fu1p, s_fu1p, hd_fu1p, stil_fu1p, k_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,model);

if nargin <2
    model = 'sgu';
end

%Impose zero bound on nominal interest rate
R = max(1,R);

%Deep Structural Parameters
[BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, KAPA, CHI, CHITIL, PSSI, PHI1, PHI2, PHI5, PHIL, GAMA1, GAMA2,  RHOMUZ, RHOMUUPSILON, RHOG, STD_EPSMUZ, STD_EPSMUUPSILON, STD_EPSG, G, MUI, MUZ, MUUPSILON ] = nberma_deep(model);

%Growth rate of output
MUZSTAR = MUZ * MUUPSILON^(THETA/(1-THETA));

%Growth rate of marginal utility of consumption
MULAMBDA = 1/ MUZSTAR;

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)); 

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

%Cost of capacity utilization
AU = GAMA1*(U-1) + GAMA2/2*(U-1)^2;

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

%Rental rate of capital
RK = GAMA1 + GAMA2 * (U-1);

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

%EHL model
if strcmp(model,'ehl')
    STIL = ((1-ALFATIL) * WTIL_O_W^(-ETATIL*(1+PHI5)) / (1-ALFATIL * (PAI*MUZSTAR)^((1-CHITIL)*ETATIL*(1+PHI5))))^(1/(1+PHI5));%distortion between hours supplied and hours used as fac
end
   
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

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 quadratic equation in HD

constant_L = (1-BETTA * B * MULAMBDA) / (1- B / MUZSTAR) / (1 + ELL + V * (PHI1 - PHI2 / V^2));

constant_f1 = (ETATIL -1 ) / ETATIL * WTIL * WTIL_O_W ^(-ETATIL ) / (1 - ALFATIL * BETTA *(MUZSTAR * PAI)^((1-CHITIL)*(ETATIL -1))*MULAMBDA * MUZSTAR);

constant_f2 = PHIL * STIL ^PHI5 * WTIL_O_W^(-ETATIL) /(1 - ALFATIL * BETTA * (MUZSTAR * PAI)^((1-CHITIL)* ETATIL));

if strcmp(model, 'ehl')
    constant_f2 = PHIL * WTIL_O_W^(-ETATIL*(1+PHI5)) /(1 - ALFATIL * BETTA * (MUZSTAR * PAI)^((1-CHITIL)* ETATIL*(1+PHI5)));
end

%this is only valid if PHI5 =1%
aa = K_U_O_HD^THETA / S - K_U_O_HD/U *(MUI -(1-DELTA) + AU);
bb = - PSSI/S - G; 
cc = - constant_L * constant_f1 / constant_f2 *( 1+ ELL) ;

HD = - bb/2/aa + sqrt(bb^2/4/aa^2 - cc/aa); %Composite labor input 

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 =  constant_L / C; %steady state marginal utility of income

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

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

F1 = constant_f1 * LA * HD; 

F2 = F1;

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

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

%Assign values to endogenous variables dated t-1, t, and t+1 (or _ba1, _cu, and _fu1, respectively)
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); 
mc_cu=(MC); 
s_cu=(S); 
hd_cu = (HD);
stil_cu=(STIL); 
k_cu=(K); 
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); 
mc_cup=(MC); 
s_cup=(S); 
hd_cup = (HD);
stil_cup=(STIL); 
k_cup=(K); 
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); 
mc_ba1=(MC); 
s_ba1=(S); 
hd_ba1 = (HD);
stil_ba1=(STIL); 
k_ba1=(K); 
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); 
mc_ba1p=(MC); 
s_ba1p=(S); 
hd_ba1p = (HD);
stil_ba1p=(STIL); 
k_ba1p=(K); 
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); 
mc_ba2=(MC); 
s_ba2=(S); 
hd_ba2 = (HD);
stil_ba2=(STIL); 
k_ba2=(K); 
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); 
mc_ba2p=(MC); 
s_ba2p=(S); 
hd_ba2p = (HD);
stil_ba2p=(STIL); 
k_ba2p=(K); 
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); 
mc_fu1=(MC); 
s_fu1=(S); 
hd_fu1 = (HD);
stil_fu1=(STIL); 
k_fu1=(K); 
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); 
mc_fu1p=(MC); 
s_fu1p=(S); 
hd_fu1p = (HD);
stil_fu1p=(STIL); 
k_fu1p=(K); 
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
%[gap] = get_gap(POLGUESS,dLdx,dOdx, model) computes the numerical error associated with a 
%particular guess for the Ramsey-optimal steady-state nominal interest rate. 
%[gap,XI] = get_gap(POLGUESS,dLdx,dOdx, model) computes the value of the Lagrange multipliers of the constraints of the Ramsey problem in the deterministic steady state.
%Inputs:
%POLGUESS = guess for the gross quarterly nominal interest rate
%dLdx = derivative of competitive equilibrium conditions
%dOdx = derivative of the Ramsey objective function
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.

function [gap,XI] = get_gap(POLGUESS,dLdx,dOdx, model);

RGUESS = max(1,POLGUESS);

%Steady-state values
[ALFA, ALFATIL, B, BETTA, CHI, CHITIL, DELTA, ETA, ETATIL, KAPA, NU, PHI1, PHI2,  PHIL, PHI5, PSSI, THETA, GAMA1, GAMA2, u_ba1, v_ba1, profits_ba1,  qq_ba1, r_ba1, ptil_ba1, mc_ba1, s_ba1, hd_ba1, stil_ba1, k_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, mc_ba1p, s_ba1p, hd_ba1p, stil_ba1p, k_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,  mc_ba2, s_ba2, hd_ba2, stil_ba2, k_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, mc_ba2p, s_ba2p, hd_ba2p, stil_ba2p, k_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,  mc_cu, s_cu, hd_cu, stil_cu, k_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, mc_cup, s_cup, hd_cup, stil_cup, k_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, mc_fu1, s_fu1,  hd_fu1, stil_fu1, k_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, mc_fu1p, s_fu1p, hd_fu1p, stil_fu1p, k_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, model);

%Evaluate derivatives
ndLdx=eval(dLdx);
ndOdx=eval(dOdx);

%Zero bound is binding
if r_cu<=1;
     XI = -ndOdx(1:end-1)*inv(ndLdx(:,1:end-1));
      gap = (XI*ndLdx+ndOdx);
     sign_right = XI * ndLdx(:,end) + ndOdx(1,end);
     if sign_right<0
        gap = norm(gap(1:end-1));
     else
        gap = 1;
     end
else %Zero bound is not binding
    %Least-Square Method
    XI = -((ndLdx')\ndOdx')';
    residual1 = ndOdx +  XI * ndLdx;
    gap = norm(residual1);
end



%NBERMA_DEEP.M
%nberma_deep(model) computes 
%the deep Structural parameters of the competitive economy developed in 
%``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by Stephanie Schmitt-Grohe and 
%Martin Uribe (2006).
%Input: 
%model = a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%The default value is 'sgu'.
%Output: the deep structural parameters 
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.

function [BETTA, THETA, DELTA, NU,  ETATIL, ETA, ALFA, ALFATIL, B, KAPA, CHI, CHITIL, PSSI, PHI1, PHI2, PHI5, PHIL, GAMA1, GAMA2,  RHOMUZ, RHOMUUPSILON, RHOG, STD_EPSMUZ, STD_EPSMUUPSILON, STD_EPSG, G, MUI, MUZ, MUUPSILON,ALFATILSGU, ALFATILEHL] = nberma_deep(model);

if nargin <1
    model = 'sgu';
end


U = 1; %SS capital utilization

GAMA2_O_GAMA1 =  1.46; %cost of capital utilization

PAISTAR = 1.042^(1/4); %Inflation target

%Exogenous processes, from ACEL
RHOMUZ = 0.89; %Serial correlation of neutral tech. shock growth, 

STD_EPSMUZ = 0.0007; %Std. dev.  of innovation to neutral tech. shock

RHOMUUPSILON = 0.20; %Serial correlation of neutral tech. shock growth, 

STD_EPSMUUPSILON = 0.0031; %Std. dev.  of innovation to neutral tech. shock

% GOV'T PROCESS FROM Ravn (2005)
RHOG = 0.9; %Serial correlation of gov't consumption 
STD_EPSG = 0.008; %Std. Dev. of  innovation to gov't consumption 


SG = 0.17; %Share of gov't expenditures in gdp (own estimate)


%Calibration of Money demand block
SMH = 0.44; %Share of household money in total money (from CEE; ratio M1/M2). 

SM = 0.1695*4; %M1/GDP. Sample: 1959:1-2004:3. Source: GDP NIPA and M1 FRB. Produced with m1_gdp.m  

EPS_MH_R = -0.81; %annualized inerest rate semielasticity of money demand (ACEL)

BETTA =  1.03^(-1/4); %Discount factor;

PHIL = 1; %Preference Parameter

PHI5 = 1; %Preference Parameter

THETA = 0.36;%Capital share; from ACEL 

DELTA = 0.025; %from ACEL 

ETATIL = 21; %Labor elasticity of subst (ACEL)

ETA = 6; %goods elasticity of substitution (ACEL)

ALFA = 0.8; %Degree of price stickiness (ACEL, rental market case)

%Degree of wage stickiness in the EHL model
ALFATILEHL = 0.69;

%Coefficient of the Wage Phillips curve
GAMA = (1-ALFATILEHL) .* (1-BETTA*ALFATILEHL) ./ ALFATILEHL/(1+ETATIL*PHI5); 

%Degree of wage stickiness in the SGU model
ALFATILSGU = 1./2./BETTA.*(GAMA+BETTA+1-(GAMA.^2+2.*GAMA.*BETTA+2.*GAMA+BETTA.^2-2.*BETTA+1).^(1./2));

if strcmp(model, 'ehl')
    ALFATIL = ALFATILEHL; %degree of wage stickiness (ACEL)
else

ALFATIL = ALFATILSGU;

 end

    
B = 0.69; %degree of habit formation (ACEL)

KAPA = 2.79; %Capital adjustment cost (ACEL)

CHI = 0; %Degree of price ndexation (Cogley and Sbordone, LOWW)

CHITIL = 1; %Degree of wage indexation (ACEL and LOWW)

MUZ = 1.00213; %steady state growth rate of neutral technological growth, ACEL%The report 1.00013, but to square with their text it should be 1.0045/1.0042^(THeTA/(1-THETA))

MUUPSILON = 1.0042; %ACEL

MUZSTAR = MUZ * MUUPSILON^(THETA/(1-THETA)); %Steady-state growth rate of output

MUI = MUUPSILON * MUZSTAR;  %steady state growth rate of investment measured in tems of investment goods

MULAMBDA = 1/ MUZSTAR;%growth rate of martinal utility of consumption

PAI = PAISTAR; %Inflation

QQ = 1; %Tobin's q

SK = THETA; %capital share

R = PAI / BETTA / MULAMBDA; %Nominal interest rate

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

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

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

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

SI = (MUI - (1-DELTA) ) * SK / (MUUPSILON/BETTA/MULAMBDA -1 + DELTA); %share of investment

RK = MUUPSILON/BETTA/MULAMBDA -1 + DELTA; %rental rate of capital

K_O_HD = (RK/MC/THETA)^(1/(THETA-1)) * MUI/U; %aggregate capital factor input

NU = SM * (1-SMH) / (1-SK - SM * (1-SMH) *(1- 1/R)); %Working capital constraint

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

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

GAMA1 = RK; %Parameter of capacity utilization cost function

GAMA2 = GAMA2_O_GAMA1 * GAMA1; %Parameter of capacity utilization cost function


PHI2=-1/8*(8*EPS_MH_R*R^2-8*EPS_MH_R*R+1)/EPS_MH_R/R^2; %Paramter of transactions cost technology

VTILDE = sqrt(PHI2+1-1/R); %auxiliary variable, V=VTILDE/ sqrt(PHI1)

SC = 1 - SI - SG; %Consumption share in gdp (this includes the transactions cost)

PHI1 = VTILDE^2 / (SC/SMH/SM +2* VTILDE* sqrt(PHI2) - VTILDE^2-PHI2)^2; %Paramter of transactions cost technology

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

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

%Construct some constants to solve for HD

C2_EHL = PHIL * WTIL_O_W^(-ETATIL*(1+PHI5))/(1-ALFATIL*BETTA*(PAI*MUZSTAR)^((1-CHITIL)*ETATIL*(1+PHI5)));

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

C2_SGU = PHIL * STIL_SGU^PHI5 * WTIL_O_W^(-ETATIL) / (1-ALFATIL*BETTA*(PAI*MUZSTAR)^((1-CHITIL)*ETATIL));

%NOTE THAT C2_SGU=C2_EHL IN OUR SETUPT (because chitil =1) THUS IT DOES NOT MATTER WHICH ONE WE
%USE IN WHAT FOLLOWS 

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

CL = (1 - B * BETTA * MULAMBDA) / (1-B/MUZSTAR) /( 1 + ELL + V * PHI1 - PHI2/V); 

HD = (C1 * CL / NU / W/C2_SGU * (1-SMH) * SM / SC *(1+ELL))^(1/(1+PHI5)); %Hours employed by firms

K = K_O_HD * HD;  %capital

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

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

OUTPUT = ((K/MUI)^THETA * HD^(1-THETA) - PSSI) / S; %aggregate demand  

C = (OUTPUT*(1-SG)-IV) / (1+ELL);%Steady state consumption, one should subtract a(u)k, but this is zero

G = SG * OUTPUT; %Government spending