%COMPETITIVE_F.M

%This program computes analytically the equilibrium conditions and its derivatives in ``Optimal,  Simple, and Implementable Monetary and Fiscal Rules,'' by S. Schmitt-Grohe and Martin Uribe (2006)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, May 2006


function competitive_f

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

%Define symbols
allsyms

%read Equilibrium conditions
eqmc = constraints;


%Interest-rate rule
rule =  crule * (- ALFA_R_CU * log(r_cu/RSTAR) + ALFA_R_BA1 * log(r_ba1/RSTAR) + ALFA_PAI * log(pai_cu/PAISTAR) + ALFA_Y * log(output_cu/OUTPUT)) + brule * (- ALFA_R_CU * log(r_cu/RSTAR) + ALFA_R_BA1 * log(r_ba1/RSTAR) + ALFA_PAI * log(pai_ba1/PAISTAR) + ALFA_Y * log(output_ba1/OUTPUT)) + frule * (- ALFA_R_CU * log(r_cu/RSTAR) + ALFA_R_BA1 * log(r_ba1/RSTAR) + ALFA_PAI * log(pai_cup/PAISTAR) + ALFA_Y * log(output_cup/OUTPUT)) + drule * (- ALFA_R_CU * log(r_cu/r_ba1) + ALFA_R_BA1 * log(r_ba1/r_ba2) +ALFA_PAI * log(pai_ba1/PAISTAR) + ALFA_Y * log(output_ba1/output_ba2));  


%welfare at time t:

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

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


%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 equation linking past-,and  current--dated investment, ouput, inflation, and the interest rate

link_bwd_fwd = [-iv_ba1p+iv_cu; -r_ba1p+r_cu; -pai_ba1p+pai_cu; -output_ba1p+output_cu; -output_ba2p+output_ba1 ; -r_ba2p+r_ba1];

f = [vt; eqmc; exproc; link_bwd_fwd; rule];

statevar_cu = [r_ba1; r_ba2; pai_ba1; output_ba1; output_ba2; iv_ba1; s_cu;  k_cu;    g_cu; z_cu]; 

statevar_cu = transpose(statevar_cu);

statevar_cup = [r_ba1p; r_ba2p; pai_ba1p; output_ba1p; output_ba2p; iv_ba1p; s_cup;  k_cup;    g_cup; z_cup]; 

statevar_cup = transpose(statevar_cup);

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

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

%variables to substitute from levels to logs (all variables but money and welfare)
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 r_cu iv_ba1 r_ba1 r_ba2 output_ba1 output_ba2 pai_ba1];

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 r_cup iv_ba1p r_ba1p r_ba2p output_ba1p output_ba2p pai_ba1p];

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