%RAMSEY_SS.M
%Ramsey steady state without fiscal policy in ``Optimal,  Simple, and Implementable Monetary and Fiscal Rules,'' by S. Schmitt-Grohe and Martin Uribe (2006)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, April 2006 

function ramsey_ss;

allsyms

%Utility function
objective_cu =  ((c_cu * (1-h_cu)^GAMA )^(1-SIGMA) - 1) / (1-SIGMA);

[constraints_cu, constraints_ba1, constraints_fu1] = constraints;

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

endogvar_cu = [taud_cu k_cup iv_cu c_cu h_cu la_cu w_cu q_cu u_cu pai_cu mc_cu m_cu ptil_cu x1_cu x2_cu output_cu s_cup r_cu ]; 

endogvar_cu=endogvar_cu(:);

endogvar_ba1p = [taud_ba1p k_fu1 iv_ba1p c_ba1p h_ba1p la_ba1p w_ba1p q_ba1p u_ba1p pai_ba1p mc_ba1p m_ba1p ptil_ba1p x1_ba1p x2_ba1p output_ba1p s_fu1 r_ba1p]; 

endogvar_ba1p=endogvar_ba1p(:);

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

O = objective_cu;

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

%Initial guess for gross nominal interest rate and the tax rate
X0=[1.01; 0.17];
%load x0


%Find the optimal policy
OPTIONS=optimset('MaxFunEvals',100000,'MaxIter',100000,'TolFun',1e-11,'TolX',1e-11);

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


%Impose non-negativity of nominal interest rate
oppolicy(1) = max(1,oppolicy(1));

[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);
xi16_cu = XI(16);
xi17_cu = XI(17);

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);
xi16_cup = XI(16);
xi17_cup = XI(17);


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);
xi16_ba1 = XI(16);
xi17_ba1 = XI(17);



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);
xi16_ba1p = XI(16);
xi17_ba1p = XI(17);



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);
xi16_fu1 = XI(16);
xi17_fu1 = XI(17);



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_fu1p = XI(16);
xi17_fu1p = XI(17);


[DELTA, SIGMA, GAMA, NUH, BETTA, THETA, NUF, ALFA, ETA, CHI, GBAR, RHOG, RHOZ, STD_EPSG, STD_EPSZ, PSSI, Q, TAUD, U, PAI, PTIL, S, G, MC, W, H, C, K, IV, LA, M, OUTPUT, X1, X2, r_cu, q_cu, taud_cu, u_cu, pai_cu, ptil_cu, g_cu, mc_cu, w_cu, h_cu, c_cu, k_cu, iv_cu, la_cu, m_cu, output_cu, x1_cu, x2_cu, s_cu, z_cu, r_cup, q_cup, taud_cup, u_cup, pai_cup, ptil_cup, g_cup, mc_cup, w_cup, h_cup, c_cup, k_cup, iv_cup, la_cup, m_cup, output_cup, x1_cup, x2_cup, s_cup, z_cup, r_ba1, q_ba1, taud_ba1, u_ba1, pai_ba1, ptil_ba1, g_ba1, mc_ba1, w_ba1, h_ba1, c_ba1, k_ba1, iv_ba1, la_ba1, m_ba1, output_ba1, x1_ba1, x2_ba1, s_ba1,  z_ba1, r_ba1p, q_ba1p, taud_ba1p, u_ba1p, pai_ba1p, ptil_ba1p, g_ba1p, mc_ba1p, w_ba1p, h_ba1p, c_ba1p, k_ba1p, iv_ba1p, la_ba1p, m_ba1p, output_ba1p, x1_ba1p, x2_ba1p, s_ba1p, z_ba1p, r_fu1, q_fu1, taud_fu1, u_fu1, pai_fu1, ptil_fu1, g_fu1, mc_fu1, w_fu1, h_fu1, c_fu1, k_fu1, iv_fu1, la_fu1, m_fu1, output_fu1, x1_fu1, x2_fu1, s_fu1, z_fu1, r_fu1p, q_fu1p, taud_fu1p, u_fu1p, pai_fu1p, ptil_fu1p, g_fu1p, mc_fu1p, w_fu1p, h_fu1p, c_fu1p, k_fu1p, iv_fu1p, la_fu1p, m_fu1p, output_fu1p, x1_fu1p, x2_fu1p, s_fu1p, z_fu1p, iv_ba2, l_cu, tau_cu, l_cup, tau_cup, l_ba1, tau_ba1, l_ba1p, tau_ba1p, l_fu1, tau_fu1, l_fu1p, tau_fu1p] = ss(oppolicy);


vt_cu = eval(objective_cu) / (1-BETTA);
vt_cup = vt_cu;
output_ba2 = output_ba1;
output_ba2p = output_ba1;
r_ba2 = r_cu;
r_ba2p = r_cu;

PAISTAR = pai_cu;
OUTPUT = output_cu;
RSTAR = r_cu;
LSTAR = l_cu;
TAUSTAR = tau_cu;

format compact

disp(['accuracy criterion, dist = ' num2str(dist)])
disp(['Exit flag = ' num2str(exitflag)])
disp(['annual interest rate = ' num2str(r_cu^4)])
disp(['annual inflation rate = ' num2str(pai_cu^4)])
disp(['tax rate = ' num2str(taud_cu)])


save ramsey_ss.mat 

%SS.M
%Steady-state given R of the model in ``Optimal, Simple, and Implementable Monetary and Fiscal Rules,'' by S. Schmitt-Grohe and Martin Uribe   

function [DELTA, SIGMA, GAMA, NUH, BETTA, THETA, NUF, ALFA, ETA, CHI, GBAR, RHOG, RHOZ, STD_EPSG, STD_EPSZ, PSSI, Q, TAUD, U, PAI, PTIL, S, G, MC, W, H, C, K, IV, LA, M, OUTPUT, X1, X2, r_cu, q_cu, taud_cu, u_cu, pai_cu, ptil_cu, g_cu, mc_cu, w_cu, h_cu, c_cu, k_cu, iv_cu, la_cu, m_cu, output_cu, x1_cu, x2_cu, s_cu, z_cu, r_cup, q_cup, taud_cup, u_cup, pai_cup, ptil_cup, g_cup, mc_cup, w_cup, h_cup, c_cup, k_cup, iv_cup, la_cup, m_cup, output_cup, x1_cup, x2_cup, s_cup, z_cup, r_ba1, q_ba1, taud_ba1, u_ba1, pai_ba1, ptil_ba1, g_ba1, mc_ba1, w_ba1, h_ba1, c_ba1, k_ba1, iv_ba1, la_ba1, m_ba1, output_ba1, x1_ba1, x2_ba1, s_ba1,  z_ba1, r_ba1p, q_ba1p, taud_ba1p, u_ba1p, pai_ba1p, ptil_ba1p, g_ba1p, mc_ba1p, w_ba1p, h_ba1p, c_ba1p, k_ba1p, iv_ba1p, la_ba1p, m_ba1p, output_ba1p, x1_ba1p, x2_ba1p, s_ba1p, z_ba1p, r_fu1, q_fu1, taud_fu1, u_fu1, pai_fu1, ptil_fu1, g_fu1, mc_fu1, w_fu1, h_fu1, c_fu1, k_fu1, iv_fu1, la_fu1, m_fu1, output_fu1, x1_fu1, x2_fu1, s_fu1, z_fu1, r_fu1p, q_fu1p, taud_fu1p, u_fu1p, pai_fu1p, ptil_fu1p, g_fu1p, mc_fu1p, w_fu1p, h_fu1p, c_fu1p, k_fu1p, iv_fu1p, la_fu1p, m_fu1p, output_fu1p, x1_fu1p, x2_fu1p, s_fu1p, z_fu1p, iv_ba2, l_cu, tau_cu, l_cup, tau_cup, l_ba1, tau_ba1, l_ba1p, tau_ba1p, l_fu1, tau_fu1, l_fu1p, tau_fu1p] = ss(pol); 

R = pol(1);
TAUD = pol(2);


%Load structural parameteres
[DELTA, SIGMA, GAMA, NUH, BETTA, THETA, NUF, ALFA, ETA, CHI, GBAR, PSSI, RHOG, STD_EPSG, RHOZ, STD_EPSZ] = structural_parameters; 

Q= 1; %Tobin's Q

U = (1/BETTA + DELTA - 1 - TAUD * DELTA) / (1- TAUD); %Rental rate of capital

PAI = BETTA * R; %Inflation 

PTIL = ((1-ALFA * PAI^(-1+ETA)) / (1 - ALFA))^(1/(1-ETA)); %Relative price of reoptimizing firm's output

S = ((1 - ALFA) * PTIL^(-ETA)) / (1-ALFA* PAI^ETA); %Price distortion

G = GBAR; %Government spending

MC = (ETA -1) / ETA * (1-ALFA * BETTA * PAI^ETA) / (1 - ALFA * BETTA * PAI^(ETA-1)) * PTIL; %Marginal cost

KAPA = (U / MC / THETA)^(1/(THETA-1)); %Capital-labor ratio

W = MC * (1-THETA) * KAPA^THETA / (1+NUF * (1-1/R)); %Wage rate

ZETAC = W * (1-TAUD) / GAMA / (1+NUH * (1- 1/R)); %Consumption-leisure ratio

H = (S * ZETAC + S*G + CHI) / (S* ZETAC - S* DELTA * KAPA + KAPA^THETA); %Hours

C = ZETAC * (1-H); %Consumption

K = KAPA * H; %Capital

IV = DELTA * K; %Investment

LA = C^(-SIGMA) * (1-H)^(GAMA*(1-SIGMA)) / (1 + NUH * (1- 1/R)); %Marginal utility of wealth

M = NUH * C + NUF * W * H; %Money

OUTPUT = C + IV + G; 

X1 = (PTIL^(-1-ETA) * OUTPUT * MC ) / (1- ALFA * BETTA * PAI^ETA); %PDV of marginal revenues

X2 = ETA/ (ETA-1) * X1; %PDV of marginal costs

TAU = TAUD * OUTPUT;

L = (R*(G-TAU) - M * (R-1)) / (1-R/PAI);

r_cu = (R);
q_cu = (Q);
taud_cu = TAUD;
u_cu = (U);
pai_cu = (PAI);
ptil_cu = (PTIL);
s_cu = (S);
g_cu = (G); 
mc_cu = (MC);
w_cu = (W); 
h_cu = (H);
c_cu = (C);
k_cu = (K);
iv_cu = (IV);
la_cu = (LA);
m_cu = (M);
output_cu = (OUTPUT);
x1_cu = (X1);
x2_cu = (X2);
z_cu = (1);
taud_cu = TAUD;
tau_cu = TAU;
l_cu = L;


r_cup = (R);
q_cup = (Q);
taud_cup = TAUD;
u_cup = (U);
pai_cup = (PAI);
ptil_cup = (PTIL);
s_cup = (S);
g_cup = (G); 
mc_cup = (MC);
w_cup = (W); 
h_cup = (H);
c_cup = (C);
k_cup = (K);
iv_cup = (IV);
la_cup = (LA);
m_cup = (M);
output_cup = (OUTPUT);
x1_cup = (X1);
x2_cup = (X2);
z_cup = (1);
taud_cup = TAUD;
tau_cup = TAU;
l_cup = L;
r_ba1 = (R);
q_ba1 = (Q);
taud_ba1 = TAUD;
u_ba1 = (U);
pai_ba1 = (PAI);
ptil_ba1 = (PTIL);
s_ba1 = (S);
g_ba1 = (G); 
mc_ba1 = (MC);
w_ba1 = (W); 
h_ba1 = (H);
c_ba1 = (C);
k_ba1 = (K);
iv_ba1 = (IV);
la_ba1 = (LA);
m_ba1 = (M);
output_ba1 = (OUTPUT);
x1_ba1 = (X1);
x2_ba1 = (X2);
z_ba1 = (1);
taud_ba1 = TAUD;
tau_ba1 = TAU;
l_ba1 = L;

r_ba1p = (R);
q_ba1p = (Q);
taud_ba1p = TAUD;
u_ba1p = (U);
pai_ba1p = (PAI);
ptil_ba1p = (PTIL);
s_ba1p = (S);
g_ba1p = (G); 
mc_ba1p = (MC);
w_ba1p = (W); 
h_ba1p = (H);
c_ba1p = (C);
k_ba1p = (K);
iv_ba1p = (IV);
la_ba1p = (LA);
m_ba1p = (M);
output_ba1p = (OUTPUT);
x1_ba1p = (X1);
x2_ba1p = (X2);
z_ba1p = (1);
taud_ba1p = TAUD;
tau_ba1p = TAU;
l_ba1p = L;

r_fu1 = (R);
q_fu1 = (Q);
taud_fu1 = TAUD;
u_fu1 = (U);
pai_fu1 = (PAI);
ptil_fu1 = (PTIL);
s_fu1 = (S);
g_fu1 = (G); 
mc_fu1 = (MC);
w_fu1 = (W); 
h_fu1 = (H);
c_fu1 = (C);
k_fu1 = (K);
iv_fu1 = (IV);
la_fu1 = (LA);
m_fu1 = (M);
output_fu1 = (OUTPUT);
x1_fu1 = (X1);
x2_fu1 = (X2);
z_fu1 = (1);
taud_fu1 = TAUD;
tau_fu1 = TAU;
l_fu1 = L;

r_fu1p = (R);
q_fu1p = (Q);
taud_fu1p = TAUD;
u_fu1p = (U);
pai_fu1p = (PAI);
ptil_fu1p = (PTIL);
s_fu1p = (S);
g_fu1p = (G); 
mc_fu1p = (MC);
w_fu1p = (W); 
h_fu1p = (H);
c_fu1p = (C);
k_fu1p = (K);
iv_fu1p = (IV);
la_fu1p = (LA);
m_fu1p = (M);
output_fu1p = (OUTPUT);
x1_fu1p = (X1);
x2_fu1p = (X2);
z_fu1p = (1);
taud_fu1p = TAUD;
tau_fu1p = TAU;
l_fu1p = L;

iv_ba2 = (IV);


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

RGUESS = max(1,POLGUESS(1));
TAUDGUESS = POLGUESS(2);



[DELTA, SIGMA, GAMA, NUH, BETTA, THETA, NUF, ALFA, ETA, CHI, GBAR, RHOG, RHOZ, STD_EPSG, STD_EPSZ, PSSI, Q, TAUD, U, PAI, PTIL, S, G, MC, W, H, C, K, IV, LA, M, OUTPUT, X1, X2, r_cu, q_cu, taud_cu, u_cu, pai_cu, ptil_cu, g_cu, mc_cu, w_cu, h_cu, c_cu, k_cu, iv_cu, la_cu, m_cu, output_cu, x1_cu, x2_cu, s_cu, z_cu, r_cup, q_cup, taud_cup, u_cup, pai_cup, ptil_cup, g_cup, mc_cup, w_cup, h_cup, c_cup, k_cup, iv_cup, la_cup, m_cup, output_cup, x1_cup, x2_cup, s_cup, z_cup, r_ba1, q_ba1, taud_ba1, u_ba1, pai_ba1, ptil_ba1, g_ba1, mc_ba1, w_ba1, h_ba1, c_ba1, k_ba1, iv_ba1, la_ba1, m_ba1, output_ba1, x1_ba1, x2_ba1, s_ba1,  z_ba1, r_ba1p, q_ba1p, taud_ba1p, u_ba1p, pai_ba1p, ptil_ba1p, g_ba1p, mc_ba1p, w_ba1p, h_ba1p, c_ba1p, k_ba1p, iv_ba1p, la_ba1p, m_ba1p, output_ba1p, x1_ba1p, x2_ba1p, s_ba1p, z_ba1p, r_fu1, q_fu1, taud_fu1, u_fu1, pai_fu1, ptil_fu1, g_fu1, mc_fu1, w_fu1, h_fu1, c_fu1, k_fu1, iv_fu1, la_fu1, m_fu1, output_fu1, x1_fu1, x2_fu1, s_fu1, z_fu1, r_fu1p, q_fu1p, taud_fu1p, u_fu1p, pai_fu1p, ptil_fu1p, g_fu1p, mc_fu1p, w_fu1p, h_fu1p, c_fu1p, k_fu1p, iv_fu1p, la_fu1p, m_fu1p, output_fu1p, x1_fu1p, x2_fu1p, s_fu1p, z_fu1p, iv_ba2, l_cu, tau_cu, l_cup, tau_cup, l_ba1, tau_ba1, l_ba1p, tau_ba1p, l_fu1, tau_fu1, l_fu1p, tau_fu1p] = ss([RGUESS; TAUDGUESS]);

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

if r_cu==1;
     XI = -(ndLdx(:,1:end-1)'\ndOdx(1:end-1)')';
      gap1 = (XI*ndLdx+ndOdx);
     sign_right = XI * ndLdx(:,end) + ndOdx(1,end);
     if sign_right<=0
        gap1 = norm(gap1(1:end-1));
     else
        gap1 = 10;
     end
else
    %Least-Square Method
    XI = -((ndLdx')\ndOdx')';
    gap1 = norm(ndOdx +  XI * ndLdx);
end
    gap2 = abs((l_cu + m_cu)/r_cu / output_cu / 4 - 0.44);
    gap = gap1+gap2