%RAMSEY_F.M Computes analytically the equilibrium conditions of the ramsey problem and its derivatives 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_f

%Set order of approximation (1 or 2)
approx = 2;

%Define symbols
allsyms

%read Ramsey FOC's
foc = ramsey_foc;

%Add processes for exogenous variables

%Government Purchases shock
exprocg = log(g_cup/GBAR) - RHOG * log(g_cu/GBAR);

%Evolution of Neutral Technolgy shock
exprocz = log(z_cup/ZBAR) - RHOZ * log(z_cu/ZBAR);

exproc = [exprocg;exprocz];

%Replace all variables that appear in foc as _fu1 but not as _fu1p by variables in _cup
fu1only = [pai_fu1  xi6_fu1 ];

fu1only2cup = [pai_cup  xi6_cup];

foc=subs(foc,fu1only, fu1only2cup);


%Replace all variables that appear in foc as _ba1p but not as _ba1 by variables in _cu
ba1ponly = [pai_ba1p  la_ba1p ];

ba1ponly2cu = [pai_cu la_cu];

foc=subs(foc,ba1ponly, ba1ponly2cu);

%Replace all variables that appear in foc as _cu and as  _ba1  but not as _cup by variables in _ba1p

cuonly = [r_cu xi4_cu xi5_cu ];

cuonly2ba1p = [r_ba1p xi4_ba1p xi5_ba1p ];

foc=subs(foc,cuonly, cuonly2ba1p);

%Period utility
uf_cu = log(c_cu) + DELTA * log(1-h_cu);

%Lifetime utility
vt = -vt_cu + uf_cu + BETTA * vt_cup;

f = [vt; foc;  exproc;-output_cu+z_cu*h_cu;-r_cu+r_ba1p];

statevar_cu = [xi4_ba1; xi5_ba1; a_cu; r_ba1;    g_cu; z_cu]; 

statevar_cu = transpose(statevar_cu);

statevar_cup = [xi4_ba1p; xi5_ba1p; a_cup; r_ba1p;    g_cup; z_cup]; 

statevar_cup = transpose(statevar_cup);

%LAGRANGE_MULTIPLIERS.M
xi_cu= [xi1_cu  xi2_cu  xi3_cu xi6_cu  xi7_cu  xi8_cu];

xi_cup = [xi1_cup  xi2_cup  xi3_cup  xi6_cup  xi7_cup  xi8_cup];

controlvar_cu = [vt_cu c_cu  output_cu h_cu m_cu v_cu mc_cu la_cu pai_cu r_cu tau_cu  xi_cu];

controlvar_cup = [vt_cup c_cup output_cup h_cup m_cup v_cup mc_cup la_cup pai_cup r_cup tau_cup  xi_cup];

%Make f a function of the logarithm of the state and control vector

%variables to substitute from levels to logs
logvar = [g_cu z_cu];

logvarp = [g_cup z_cup];

lv = [logvar logvarp];

f = subs(f, lv, exp(lv));

%Compute analytical derivatives of f
[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,statevar_cu,controlvar_cu,statevar_cup,controlvar_cup,approx);

%Make f and its derivatives a function of the level of its arguments rather than the log
f = subs(f, lv, log(lv));

fx = subs(fx, lv, log(lv));
fy = subs(fy, lv, log(lv));
fxp = subs(fxp, lv, log(lv));
fyp = subs(fyp, lv, log(lv));


if approx==2

fypyp(:) = subs(fypyp(:), lv, log(lv));
fypy(:) = subs(fypy(:), lv, log(lv));
fypxp(:) = subs(fypxp(:), lv, log(lv));
fypx(:) = subs(fypx(:), lv, log(lv));
fyyp = permute(fypy,[1 3 2]);
fyy(:) = subs(fyy(:), lv, log(lv));
fyxp(:) = subs(fyxp(:), lv, log(lv));
fyx(:) = subs(fyx(:), lv, log(lv));
fxpyp = permute(fypxp,[1 3 2]);
fxpy = permute(fyxp,[1 3 2]);
fxpxp(:) = subs(fxpxp(:), lv, log(lv));
fxpx(:) = subs(fxpx(:), lv, log(lv));
fxyp = permute(fypx,[1 3 2]);
fxy = permute(fyx,[1 3 2]);
fxxp = permute(fxpx,[1 3 2]);
fxx(:) = subs(fxx(:), lv, log(lv));

end %If approx ==2
save ramsey_f.mat controlvar_cu statevar_cu f fx fxp fy fyp fypyp fypy fypxp fypx fyyp fyy fyxp fyx fxpyp fxpy fxpxp fxpx fxyp fxy fxxp fxx lv 



%RAMSEY_FOC.M
function foc = ramsey_foc;
%foc = ramsey_foc: first-order conditions of the Ramsey problem 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

%Notes
% x_cu=x_t and x_cup=x_{t+1}
% x_fu1=x_{t+1} 
%x_fu1p=x_{t+2}
% x_ba1=x_{t-1} 
% x_ba1p=x_t
% x_ba2=x_{t-2} and x_ba2p=x_{t-1}

%define symbolic variables
allsyms


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

[constraints_cu, constraints_ba1, constraints_fu1] = constraints;

%The Lagrange Multipliers
xi_cu  = [xi1_cu  xi2_cu  xi3_cu  xi4_cu  xi5_cu  xi6_cu  xi7_cu  xi8_cu];

xi_fu1  = [xi1_fu1  xi2_fu1  xi3_fu1  xi4_fu1  xi5_fu1  xi6_fu1  xi7_fu1  xi8_fu1];

xi_ba1  = [xi1_ba1  xi2_ba1  xi3_ba1  xi4_ba1  xi5_ba1  xi6_ba1  xi7_ba1  xi8_ba1];


%The Lagrangian
lagrangian = BETTA^(-1) * (xi_ba1 * constraints_ba1) + (objective_cu+xi_cu * constraints_cu) + BETTA * (xi_fu1 * constraints_fu1);

endogvar_cu = [c_cu h_cu m_cu v_cu a_cup mc_cu la_cu pai_cu tau_cu r_cu]; 

endogvar_cu=endogvar_cu(:);

endogvar_ba1p = [c_ba1p h_ba1p m_ba1p v_ba1p a_fu1 mc_ba1p la_ba1p pai_ba1p tau_ba1p r_ba1p]; 

endogvar_ba1p=endogvar_ba1p(:);

foc = [jacobian(lagrangian,xi_cu) jacobian(lagrangian,endogvar_cu)+jacobian(lagrangian,endogvar_ba1p)];

foc = foc(:);


%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];