%RAMSEY_SS.M
%Ramsey steady state under income taxation in ``Optimal Fiscal and Monetary Policy Under Sticky Prices'' by S. Schmitt-Grohe and Martin Uribe (JET, February 2004)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, May 2006 
function ramsey_ss;

allsyms

%Utility function
objective_cu =  log(c_cu) + DELTA * log(1-h_cu);

[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 = [c_cu h_cu m_cu v_cu mc_cu la_cu pai_cu r_cu tau_cu]; 

endogvar_cu=endogvar_cu(:);

endogvar_ba1p = [c_ba1p h_ba1p m_ba1p v_ba1p mc_ba1p la_ba1p pai_ba1p r_ba1p tau_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, X0(1), and tax rate X0(2)
X0 = [1  0.18];

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

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

dist

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

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

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

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

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

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

[BETTA, DELTA, GBAR, ZBAR, ETA, THETA, B, A, SIG, RHOG, RHOZ,  STD_EPSG, STD_EPSZ, DEBT_GDP_RATIO, c_cu,   la_cu, v_cu, h_cu, tau_cu, z_cu, mc_cu, r_cu, m_cu, a_cu, pai_cu, g_cu, vt_cu, c_cup,   la_cup, v_cup, h_cup, tau_cup, z_cup, mc_cup, r_cup, m_cup, a_cup, pai_cup, g_cup, vt_cup, c_ba1, la_ba1, v_ba1, h_ba1, tau_ba1, z_ba1, mc_ba1, r_ba1, m_ba1, a_ba1, pai_ba1, g_ba1,vt_ba1, c_ba1p, la_ba1p, v_ba1p, h_ba1p, tau_ba1p, z_ba1p, mc_ba1p, r_ba1p, m_ba1p, a_ba1p, pai_ba1p, g_ba1p, vt_ba1p, c_fu1, la_fu1, v_fu1, h_fu1, tau_fu1, z_fu1, mc_fu1, r_fu1, m_fu1, a_fu1, pai_fu1, g_fu1,vt_fu1, c_fu1p, la_fu1p, v_fu1p, h_fu1p, tau_fu1p, z_fu1p, mc_fu1p, r_fu1p, m_fu1p, a_fu1p, pai_fu1p, g_fu1p, vt_fu1p] = steady_state(oppolicy);

output_cu = z_cu * h_cu;
output_cup = output_cu;

save ramsey_ss



%STEADY_STATE.M
%Steady-state given R and TAU of the model in ``Optimal, Fiscal and Monetary Policy Under Sticky Prices,'' by S. Schmitt-Grohe and Martin Uribe (JET, February 2004).

    
function [BETTA, DELTA, GBAR, ZBAR, ETA, THETA, B, A, SIG, RHOG, RHOZ,  STD_EPSG, STD_EPSZ, DEBT_GDP_RATIO, c_cu, la_cu, v_cu, h_cu, tau_cu, z_cu, mc_cu, r_cu, m_cu, a_cu, pai_cu, g_cu, vt_cu, c_cup, la_cup, v_cup, h_cup, tau_cup, z_cup, mc_cup, r_cup, m_cup, a_cup, pai_cup, g_cup, vt_cup, c_ba1, la_ba1, v_ba1, h_ba1, tau_ba1, z_ba1, mc_ba1, r_ba1, m_ba1, a_ba1, pai_ba1, g_ba1,vt_ba1, c_ba1p, la_ba1p, v_ba1p, h_ba1p, tau_ba1p, z_ba1p, mc_ba1p, r_ba1p, m_ba1p, a_ba1p, pai_ba1p, g_ba1p, vt_ba1p, c_fu1, la_fu1, v_fu1, h_fu1, tau_fu1, z_fu1, mc_fu1, r_fu1, m_fu1, a_fu1, pai_fu1, g_fu1,vt_fu1, c_fu1p, la_fu1p, v_fu1p, h_fu1p, tau_fu1p, z_fu1p, mc_fu1p, r_fu1p, m_fu1p, a_fu1p, pai_fu1p, g_fu1p, vt_fu1p] = steady_state(x);

    
%Load structural parameteres
[BETTA,DELTA,GBAR,ZBAR, ETA,THETA,B,A,SIG,RHOG,RHOZ, STD_EPSG, STD_EPSZ,DEBT_GDP_RATIO] = structural_parameters;

r_cu = x(1);
tau_cu = x(2);

z_cu = ZBAR;

g_cu = GBAR;

pai_cu = BETTA * r_cu;

v_cu = sqrt(( (r_cu - 1) / r_cu + B) / A );

% introducing variables to solve the quadratic equation for h_cu
h1 = (DELTA * z_cu * ( 1 + 2 * A * v_cu - 2 * sqrt(A *B))) / ( 1 + A * v_cu + B / v_cu - 2 * sqrt (A*B));

h2 = (1 - tau_cu) * z_cu * (1 + ETA) / ETA;

h3 = (1 - tau_cu) * z_cu * pai_cu * (pai_cu -1) * THETA * (1 - BETTA) / ETA;

h4 = DELTA * (g_cu + THETA/2  * ( pai_cu - 1)^2 ) * ( 1 + 2 * A * v_cu - 2 * sqrt(A *B)) / ( 1 + A * v_cu + B / v_cu - 2 * sqrt (A*B));

h11 = h1 + h2 ;

h22 = - h4 - h3 - h2;

h33= h3;


% here we choose the root that works (+ or -)
h_cu = (-h22 + sqrt( (h22)^2 - 4 * h11 * h33))/ (2 * h11);

mc_cu = ( 1 + ETA )/ ETA - (pai_cu * (pai_cu - 1) * THETA * ( 1- BETTA)) / ( ETA * z_cu * h_cu );

c_cu = (z_cu  * h_cu - g_cu - THETA / 2 * ( pai_cu -1 )^2) / ( 1 + A * v_cu + B / v_cu - 2 * sqrt (A*B));

la_cu = 1 / ( c_cu * ( 1 + 2 * A * v_cu - 2 * sqrt(A *B)) );

m_cu = c_cu / v_cu;

a_cu = (m_cu * ( 1 - 1/ r_cu) + tau_cu * z_cu * mc_cu * h_cu -g_cu) / ( 1 / pai_cu - 1 / r_cu);


vt_cu = (log(c_cu) + DELTA * log(1-h_cu)) / (1-BETTA); %Lifetime utility

% CURRENT 
c_cup = c_cu;

la_cup = la_cu;

v_cup = v_cu;

h_cup = h_cu;

tau_cup = tau_cu;

z_cup = z_cu;

mc_cup = mc_cu;

r_cup = r_cu;

m_cup = m_cu;

a_cup = a_cu;

pai_cup = pai_cu;

g_cup = g_cu;

vt_cup = vt_cu;



% PAST
c_ba1 = c_cu;

la_ba1 = la_cu;

v_ba1 = v_cu;

h_ba1 = h_cu;

tau_ba1 = tau_cu;

z_ba1 = z_cu;

mc_ba1 = mc_cu;

r_ba1 = r_cu;

m_ba1 = m_cu;

a_ba1 = a_cu;

pai_ba1 = pai_cu;

g_ba1 = g_cu;

vt_ba1 = vt_cu;


c_ba1p = c_cu;

la_ba1p = la_cu;

v_ba1p = v_cu;

h_ba1p = h_cu;

tau_ba1p = tau_cu;

z_ba1p = z_cu;

mc_ba1p = mc_cu;

r_ba1p = r_cu;

m_ba1p = m_cu;

a_ba1p = a_cu;

pai_ba1p = pai_cu;

g_ba1p = g_cu;

vt_ba1p = vt_cu;


% FUTURE
c_fu1 = c_cu;

la_fu1 = la_cu;

v_fu1  = v_cu;

h_fu1  = h_cu;

tau_fu1 = tau_cu; 

z_fu1  = z_cu;

mc_fu1 = mc_cu;

r_fu1  = r_cu;

m_fu1 = m_cu;

a_fu1  = a_cu;

pai_fu1 = pai_cu;

g_fu1  = g_cu;

vt_fu1 = vt_cu;


c_fu1p = c_cu;

la_fu1p = la_cu;

v_fu1p  = v_cu;

h_fu1p  = h_cu;

tau_fu1p = tau_cu; 

z_fu1p  = z_cu;

mc_fu1p = mc_cu;

r_fu1p  = r_cu;

m_fu1p = m_cu;

a_fu1p  = a_cu;

pai_fu1p = pai_cu;

g_fu1p  = g_cu;

vt_fu1p = vt_cu;

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

[BETTA, DELTA, GBAR, ZBAR, ETA, THETA, B, A, SIG, RHOG, RHOZ,  STD_EPSG, STD_EPSZ, DEBT_GDP_RATIO, c_cu, la_cu, v_cu, h_cu, tau_cu, z_cu, mc_cu, r_cu, m_cu, a_cu, pai_cu, g_cu, vt_cu, c_cup, la_cup, v_cup, h_cup, tau_cup, z_cup, mc_cup, r_cup, m_cup, a_cup, pai_cup, g_cup, vt_cup, c_ba1, la_ba1, v_ba1, h_ba1, tau_ba1, z_ba1, mc_ba1, r_ba1, m_ba1, a_ba1, pai_ba1, g_ba1,vt_ba1, c_ba1p, la_ba1p, v_ba1p, h_ba1p, tau_ba1p, z_ba1p, mc_ba1p, r_ba1p, m_ba1p, a_ba1p, pai_ba1p, g_ba1p, vt_ba1p, c_fu1, la_fu1, v_fu1, h_fu1, tau_fu1, z_fu1, mc_fu1, r_fu1, m_fu1, a_fu1, pai_fu1, g_fu1,vt_fu1, c_fu1p, la_fu1p, v_fu1p, h_fu1p, tau_fu1p, z_fu1p, mc_fu1p, r_fu1p, m_fu1p, a_fu1p, pai_fu1p, g_fu1p, vt_fu1p] = steady_state(POLGUESS);

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

%Least-Square Method

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

residual1 = ndOdx +  XI * ndLdx;

residual2 = (a_cu-m_cu)/(z_cu*h_cu) / r_cu -DEBT_GDP_RATIO; 

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


%CONSTRAINTS.M
%Equilibrium conditions (constraints of the Ramsey problem) of the model in ``Optimal Fiscal and Monetary Policy Under Sticky Prices'' by S. Schmitt-Grohe and Martin Uribe (JET, February 2004)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, May 2006 


function [constraints_cu, constraints_ba1, constraints_fu1] = constraints;

allsyms

%CONSTRAINTS_CURRENT

e1_cu = - 1 / c_cu + la_cu * ( 1 + 2 * A * v_cu - 2 * sqrt(A*B) );

e2_cu = - DELTA * c_cu / (1-h_cu) + (1- tau_cu) * z_cu * mc_cu   / ( 1 +2 * A * v_cu - 2 * sqrt( A*B));

e3_cu = - A * (v_cu)^2 + B + (r_cu - 1)/r_cu;

e4_cu = - la_cu + BETTA * r_cu * la_cup / pai_cup;

e5_cu = - la_cu * pai_cu * (pai_cu - 1) + BETTA * la_cup * pai_cup * (pai_cup - 1) + la_cu * ETA * z_cu * h_cu / THETA * ( (1+ETA)/ETA - mc_cu);

e6_cu = - m_cu * ( 1- 1/r_cu) - a_cup / r_cu - tau_cu * z_cu * mc_cu * h_cu + a_cu / pai_cu + g_cu; 

e7_cu = - ( 1 + A * v_cu + B / v_cu - 2 * sqrt (A*B)) * c_cu - g_cu - THETA/2 *( pai_cu -1)^2 + z_cu * h_cu;

e8_cu = - v_cu + c_cu / m_cu;

%Create vector of current constraints
constraints_cu = [e1_cu; -e2_cu; e3_cu; e4_cu; e5_cu; e6_cu;  -e7_cu; e8_cu];


%CONSTRAINTS_PAST
e1_ba1 = - 1 / c_ba1 + la_ba1 * ( 1 + 2 * A * v_ba1 - 2 * sqrt(A*B) );

e2_ba1 = - DELTA * c_ba1 / (1-h_ba1) + (1- tau_ba1) * z_ba1 * mc_ba1  / ( 1 +2 * A * v_ba1 - 2 * sqrt( A*B));

e3_ba1 = - A * (v_ba1)^2 + B + (r_ba1 - 1)/r_ba1;

e4_ba1 = - la_ba1 + BETTA * r_ba1 * la_ba1p / pai_ba1p;

e5_ba1 = - la_ba1 * pai_ba1 * (pai_ba1 - 1) + BETTA * la_ba1p * pai_ba1p * (pai_ba1p - 1) + la_ba1 * ETA * z_ba1 * h_ba1 / THETA * ( (1+ETA)/ETA - mc_ba1);

e6_ba1 = - m_ba1 * ( 1- 1/r_ba1) - a_ba1p / r_ba1 - tau_ba1 * z_ba1 * mc_ba1 * h_ba1 + a_ba1 / pai_ba1 + g_ba1; 

e7_ba1 = - ( 1 + A * v_ba1 + B / v_ba1 - 2 * sqrt (A*B)) * c_ba1 - g_ba1 - THETA/2 *( pai_ba1 -1)^2 + z_ba1 * h_ba1;

e8_ba1 = - v_ba1 + c_ba1 / m_ba1;

%Create vector of past constraints
constraints_ba1 = [e1_ba1; -e2_ba1; e3_ba1; e4_ba1; e5_ba1; e6_ba1;  -e7_ba1; e8_ba1];


%CONSTRAINTS_FUTURE

e1_fu1 = - 1 / c_fu1 + la_fu1  * ( 1 + 2 * A * v_fu1 - 2 * sqrt(A*B) );

e2_fu1 = - DELTA * c_fu1 / (1-h_fu1) + (1- tau_fu1) * z_fu1 * mc_fu1  / ( 1 +2 * A * v_fu1 - 2 * sqrt( A*B));

e3_fu1 = - A * (v_fu1)^2 + B + (r_fu1 - 1)/r_fu1;

e4_fu1 = - la_fu1 + BETTA * r_fu1 * la_fu1p / pai_fu1p;

e5_fu1 = - la_fu1 * pai_fu1 * (pai_fu1 - 1) + BETTA * la_fu1p * pai_fu1p * (pai_fu1p - 1) + la_fu1 * ETA * z_fu1 * h_fu1 / THETA * ( (1+ETA)/ETA - mc_fu1);

e6_fu1 = - m_fu1 * ( 1- 1/r_fu1) - a_fu1p / r_fu1 - tau_fu1 * z_fu1 * mc_fu1 * h_fu1 + a_fu1 / pai_fu1 + g_fu1; 

e7_fu1 = - ( 1 + A * v_fu1 + B / v_fu1 - 2 * sqrt (A*B)) * c_fu1 - g_fu1 - THETA/2 *( pai_fu1 -1)^2 + z_fu1 * h_fu1;

e8_fu1 = - v_fu1 + c_fu1 / m_fu1;


%Create vector of future constraints
constraints_fu1 = [e1_fu1; -e2_fu1; e3_fu1; e4_fu1; e5_fu1; e6_fu1;  -e7_fu1; e8_fu1];



%STRUCTURAL_PARAMETERS.M
function [BETTA,DELTA,GBAR,ZBAR, ETA,THETA,B,A,SIG,RHOG,RHOZ, STD_EPSG, STD_EPSZ,DEBT_GDP_RATIO] = structural_parameters;
%
%Structural Parameters

H=0.2; %Used only to identify DELTA

BETTA=1/1.04; %Discount factor

SIG = 1; %scalar scaling standard deviations of shocks

%U(c,h)=log c + DELTA log (1-h) 
DELTA =2.9;

GBAR = .2*.2; %20% of GDP

ZBAR = 1; %Steady-state productivity shock

DEBT_GDP_RATIO = .44; %44% of GDP

ETA = 1.2/(1-1.2); %1+MARKUP=ETA/(1+ETA);
THETA=(-17.5/4*(1+ETA)*H); %ADJUSTMENT COST OF INFLATION PARAMETER

%Transaction cost technology s(v)
A=0.0111;
B=0.07524; %this assumes s(v) =A v + B/v -2sqrt(AB); see RATS program data/money.prg

RHOG = 0.9;
STD_EPSG = 0.03018950178319;
RHOZ = 0.82;
STD_EPSZ = 0.02289454083401;