%RAMSEY_F.M

%This program computes analytically the equilibrium conditions of the ramsey problem and its derivatives in ``Optimal,  Simple, and Implementable Monetary and Fiscal Rules,'' by S. Schmitt-Grohe and Martin Uribe (2006)
%Modified on May 12, 2011 to avoid use the function permute, which had a bug in matlab releases R2009 and R2010. This bug was 
%supposedly fixed in R2011. 
%(c) Stephanie Schmitt-Grohe and Martin Uribe, June 2006, 2011


function ramsey_f

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

%Define symbols
allsyms %run allsyms.m

%read Ramsey FOC's
foc = ramsey_foc;

%Replace all variables that appear in foc as _fu1 but not as _fu1p by variables in _cup
fu1only = [r_fu1 h_fu1 iv_fu1 k_fu1 la_fu1 mc_fu1 pai_fu1 q_fu1 xi14_fu1 xi16_fu1 xi1_fu1 xi4_fu1 xi7_fu1 xi8_fu1 xi17_fu1 z_fu1];

fu1only2cup = [r_cup h_cup iv_cup k_cup la_cup mc_cup pai_cup q_cup xi14_cup xi16_cup xi1_cup xi4_cup xi7_cup xi8_cup xi17_cup z_cup ];

foc=subs(foc,fu1only, fu1only2cup);

%Replace all variables that appear in foc as _ba1p but not as _ba1 by variables in _cu
ba1ponly = [taud_ba1p pai_ba1p q_ba1p u_ba1p x1_ba1p x2_ba1p ];

ba1ponly2cu = [taud_cu pai_cu q_cu u_cu x1_cu x2_cu ];

foc=subs(foc,ba1ponly, ba1ponly2cu);


%The variables r, xi11, xi12, xi5, xi6,  appears only as _cu and _ba1

cuonly = [xi11_cu xi12_cu xi5_cu xi6_cu];

cuonly2ba1p = [xi11_ba1p xi12_ba1p xi5_ba1p xi6_ba1p];

foc =subs(foc,cuonly, cuonly2ba1p);

%welfare at time t:

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

%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) - RHOZ * log(z_cu);

exproc = [exprocg;exprocz];

%Add equations linking past-, current-, and future-dated variables

[link_bwd_fwd, variables_ba1, variables_ba1p] = link_ba_cu_fu;

vt = -vt_cu + uf_cu + BETTA * vt_cup;

f = [vt; foc;  exproc; link_bwd_fwd];

cuonly2ba1 = [xi11_ba1 xi12_ba1 xi5_ba1 xi6_ba1];

statevar_cu = [variables_ba1; cuonly2ba1(:); s_cu;  k_cu;  l_cu;  g_cu; z_cu]; 

statevar_cu = transpose(statevar_cu);

statevar_cup = [variables_ba1p; cuonly2ba1p(:); s_cup;  k_cup;  l_cup;  g_cup; z_cup]; 

statevar_cup = transpose(statevar_cup);

%LAGRANGE_MULTIPLIERS.M
xi_cu= [xi1_cu  xi2_cu  xi3_cu xi4_cu  xi7_cu  xi8_cu  xi9_cu  xi10_cu   xi13_cu xi14_cu xi15_cu  xi16_cu xi17_cu];

xi_cup = [xi1_cup  xi2_cup  xi3_cup xi4_cup  xi7_cup  xi8_cup  xi9_cup  xi10_cup   xi13_cup xi14_cup xi15_cup  xi16_cup xi17_cup];

controlvar_cu = [vt_cu taud_cu c_cu h_cu iv_cu la_cu m_cu mc_cu output_cu pai_cu ptil_cu q_cu r_cu  u_cu w_cu x1_cu x2_cu xi_cu];


controlvar_cup = [vt_cup taud_cup c_cup h_cup iv_cup la_cup m_cup mc_cup output_cup pai_cup ptil_cup q_cup  r_cup u_cup w_cup x1_cup x2_cup xi_cup];

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

%variables to substitute from levels to logs
logvar = [c_cu h_cu iv_cu la_cu mc_cu output_cu pai_cu ptil_cu q_cu u_cu w_cu x1_cu x2_cu s_cu k_cu g_cu z_cu iv_ba1 la_ba1 ptil_ba1 r_ba1 r_cu];

logvarp = [c_cup h_cup iv_cup la_cup mc_cup output_cup pai_cup ptil_cup q_cup u_cup w_cup x1_cup x2_cup s_cup k_cup g_cup z_cup iv_ba1p la_ba1p ptil_ba1p r_ba1p r_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);


disp('done with anal_deriv')

%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]);
fyyp(:) = subs(fyyp(:), lv, log(lv)); %%%
fyy(:) = subs(fyy(:), lv, log(lv));
fyxp(:) = subs(fyxp(:), lv, log(lv));
fyx(:) = subs(fyx(:), lv, log(lv));
%fxpyp = permute(fypxp,[1 3 2]);
fxpyp(:) = subs(fxpyp(:), lv, log(lv)); %%%
%fxpy = permute(fyxp,[1 3 2]);
fxpy(:) = subs(fxpy(:), lv, log(lv)); %%%
fxpxp(:) = subs(fxpxp(:), lv, log(lv));
fxpx(:) = subs(fxpx(:), lv, log(lv));
%fxyp = permute(fypx,[1 3 2]);
fxyp(:) = subs(fxyp(:), lv, log(lv)); %%%
%fxy = permute(fyx,[1 3 2]);
fxy(:) = subs(fxy(:), lv, log(lv)); %%%
%fxxp = permute(fxpx,[1 3 2]);
fxxp(:) = subs(fxxp(:), lv, log(lv)); %%%
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