function ramsey_ir
%RAMSEY_IR.M
%Compute and plot Ramsey-optimal impulse response functions
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2005.
%
%To choose a shock for plotting impulse responses go to the line labeled `Choose shock' below and set:
%IRR = IRRmuz; for neutral  technology shock
%IRR = IRRmuupsilon; for investment-specific  technology shock
%IRR = IRRg; for gov't spending shock
%Calls: ramsey_ss.m, constraints.m, allsyms.m
%Subfunctions included in this file:
%ramsey_run.m (first-order approximation of Ramsey equilibrium)
%ramsey_f.m (analytical differentiation of Ramsey equilibrium conditions)
%link_ba_cu_fu.m (Convert linearized equilibrim conditions into a first-order vector difference equation) 
%ramsey_foc.m (analytical first-order conditions of the Ramsey problem)

clc

ramsey_run;
load ramsey_run %generated by running RAMSEY_RUN.M

nx=size(hx,1);
ny=size(gx,1);

T=40; %Length of Impulse Responses

var_cu = [controlvar_cu statevar_cu];

%Index of variables

nqq = find(var_cu=='qq_cu');
nhd = find(var_cu=='hd_cu');
noutput = find(var_cu=='output_cu');
nc = find(var_cu=='c_cu');
nla = find(var_cu=='la_cu');
npai = find(var_cu=='pai_cu');
nrk = find(var_cu=='rk_cu');
nf2 = find(var_cu=='f2_cu');
nx2 = find(var_cu=='x2_cu');
nr = find(var_cu=='r_cu');
nrback=find(var_cu=='r_ba1');
niv = find(var_cu=='iv_cu');
nw = find(var_cu=='w_cu');
ns = find(var_cu=='s_cu');
nstil = find(var_cu=='stil_cu');
nk = find(var_cu=='k_cu');
ng = find(var_cu=='g_cu');
nmuz = find(var_cu=='muz_cu');
nmuupsilon = find(var_cu=='muupsilon_cu');
nu = find(var_cu=='u_cu');
nvt = find(var_cu=='vt_cu');

%Initial condition (IC)

%IC to a neutral technology shock (1%)
x0muz = zeros(nx,1);
x0muz(find(statevar_cu=='muz_cu')) = 0.01; 

%IC to an investment specific technology shock (1%) 
x0muupsilon = zeros(nx,1);
x0muupsilon(find(statevar_cu=='muupsilon_cu')) = 0.01; 

%IC to a gov't spending shock (1%) 
x0g = zeros(nx,1);
x0g(find(statevar_cu=='g_cu')) = 0.01;

%Generate Impulse responses

for j=1:T

%IR to neutral technology shock
IRRmuz(j,1:ny) = (gx*x0muz)';
IRRmuz(j,ny+1:ny+nx) = x0muz';
x0muz = hx*x0muz;


%IR to an investment-specific technology shock
IRRmuupsilon(j,1:ny) = (gx*x0muupsilon)';
IRRmuupsilon(j,ny+1:ny+nx) = x0muupsilon';
x0muupsilon = hx*x0muupsilon;


%IR to gov't spending shock
IRRg(j,1:ny) = (gx*x0g)';
IRRg(j,ny+1:ny+nx) = x0g';
x0g = hx*x0g;

end %for j=1:T

%Choose shock
IRR = IRRmuz;

muz_ir = IRR(:,nmuz);

muupsilon_ir=IRR(:,nmuupsilon);

muzstar_ir = (THETA/(1-THETA)) * muupsilon_ir + muz_ir;

output_ir = IRR(:,noutput);
goutput_ir = cumsum(diff([0;output_ir]) + muzstar_ir)*100;

w_ir = IRR(:,nw);
gw_ir = cumsum(diff([0;w_ir]) + muzstar_ir)*100;

iv_ir = IRR(:,niv);
giv_ir = cumsum(diff([0;iv_ir]) + muzstar_ir)*100;

c_ir = IRR(:,nc);
gc_ir = cumsum(diff([0;c_ir]) + muzstar_ir)*100;


hd_ir = IRR(:,nhd)*100;

pai_ir = ((1+IRR(:,npai)+ pai_cu).^4-1)*100;

r_ir = ((1+IRR(:,nr)+ r_cu).^4-1)*100;

realr_ir = [((MUZSTAR/BETTA)^4-1)*100;r_ir(1:end-1)]-pai_ir;

u_ir = IRR(:,nu)*100;

%PLOTS
t=(0:T-1)';

subplot(4,2,1)
plot(t,goutput_ir)
title('Output')

subplot(4,2,2)
plot(t,gc_ir)
title('Consumption')

subplot(4,2,3)
plot(t,giv_ir)
title('Investment')

subplot(4,2,4)
plot(t,hd_ir)
title('Hours')

subplot(4,2,5)
plot(t,gw_ir)
title(' Real wage')

subplot(4,2,6)
plot(t,u_ir)
title('Capacity Utilization')

subplot(4,2,7)
plot(t,r_ir)
title('Nominal Interest Rate')

subplot(4,2,8)
plot(t,pai_ir)
title('Inflation')

if max(abs(IRR(:)-IRRmuz(:)))==0
disp('Impulse Response to a 1-percentage-point increase in the growth rate of the NEUTRAL TECHNOLOGY SHOCK. Figure 5 in the paper')

elseif max(abs(IRR(:)-IRRg(:)))==0
disp('Impulse Response to a 1-percent increase in the GOVERNMENT SPENDING SHOCK. Figure 6 in the paper')


elseif max(abs(IRR(:)-IRRmuupsilon(:)))==0
disp('Impulse Response to a 1-percentage-point increase in the growth rate of the INVESTMENT-SPECIFIC TECHNOLOGY SHOCK. Figure 7 in the paper')

end

shg

function ramsey_run
%RAMSEY_RUN.M
%Compute an approximation approximation of the Ramsey equilibrium conditions of the model developed in ``Optimal Inflation Stabilization in a Medium-Scale Macroeconomic Model,'' by Stephanie Schmitt-Grohe and Martin Uribe (2005)


%Analytical derivatives of Ramsey equilibrium conditions
ramsey_f;
load ramsey_f

%Compute Ramsey Steady State and save it as mat file
ramsey_ss;
load ramsey_ss 

%find the steady state of welfare
h_cu = stil_cup * hd_cu;

muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu;

uf_cu = log(c_cu - B * c_ba1/muzstar_cu)*(1-PHI4)+log(1-h_cu)*PHI4;

vt_cu =   uf_cu/(1- BETTA) ; 
vt_cup = vt_cu;


clear h_cu uf_cu muzstar_cu

%make all variables into logs
u_cu =log(u_cu);
muz_cu = log(muz_cu); 
muupsilon_cu = log(muupsilon_cu); 
qq_cu = log(qq_cu); 
r_cu = log(r_cu);
s_cu = log(s_cu); 
hd_cu = log(hd_cu);
stil_cu = log(stil_cu); 
k_cu = log(k_cu); 
iv_cu = log(iv_cu);
w_cu = log(w_cu);
output_cu = log(output_cu);
g_cu = log(g_cu);
c_cu = log(c_cu);
la_cu = log(la_cu); 
pai_cu = log(pai_cu);
f2_cu = log(f2_cu);
x2_cu = log(x2_cu);

xi1_cu = log(xi1_cu);
xi2_cu = log(xi2_cu);
xi3_cu = log(xi3_cu);
xi4_cu = log(xi4_cu);
xi5_cu = log(xi5_cu);
xi6_cu = log(xi6_cu);
xi7_cu = log(xi7_cu);
xi8_cu = log(xi8_cu);
xi9_cu = log(xi9_cu);
xi10_cu = log(xi10_cu);
xi11_cu = log(xi11_cu);
xi12_cu = log(xi12_cu);
xi13_cu = log(xi13_cu);
xi14_cu = log(xi14_cu);

u_cup =log(u_cup);
muz_cup = log(muz_cup); 
muupsilon_cup = log(muupsilon_cup);
qq_cup = log(qq_cup); 
r_cup = log(r_cup);
s_cup = log(s_cup); 
hd_cup = log(hd_cup);
stil_cup = log(stil_cup); 
k_cup = log(k_cup); 
iv_cup = log(iv_cup);
w_cup = log(w_cup);
output_cup = log(output_cup);
g_cup = log(g_cup);
c_cup = log(c_cup);
la_cup = log(la_cup); 
pai_cup = log(pai_cup);
f2_cup = log(f2_cup);
x2_cup = log(x2_cup);

xi1_cup = log(xi1_cup);
xi2_cup = log(xi2_cup);
xi3_cup = log(xi3_cup);
xi4_cup = log(xi4_cup);
xi5_cup = log(xi5_cup);
xi6_cup = log(xi6_cup);
xi7_cup = log(xi7_cup);
xi8_cup = log(xi8_cup);
xi9_cup = log(xi9_cup);
xi10_cup = log(xi10_cup);
xi11_cup = log(xi11_cup);
xi12_cup = log(xi12_cup);
xi13_cup = log(xi13_cup);
xi14_cup = log(xi14_cup);

u_ba1 =log(u_ba1);
muz_ba1 = log(muz_ba1); 
muupsilon_ba1 = log(muupsilon_ba1);
qq_ba1 = log(qq_ba1); 
r_ba1 = log(r_ba1);
s_ba1 = log(s_ba1); 
hd_ba1 = log(hd_ba1);
stil_ba1 = log(stil_ba1); 
k_ba1 = log(k_ba1); 
iv_ba1 = log(iv_ba1);
w_ba1 = log(w_ba1);
output_ba1 = log(output_ba1);
g_ba1 = log(g_ba1);
c_ba1 = log(c_ba1);
la_ba1 = log(la_ba1); 
pai_ba1 = log(pai_ba1);
f2_ba1 = log(f2_ba1);
x2_ba1 = log(x2_ba1);

xi1_ba1 = log(xi1_ba1);
xi2_ba1 = log(xi2_ba1);
xi3_ba1 = log(xi3_ba1);
xi4_ba1 = log(xi4_ba1);
xi5_ba1 = log(xi5_ba1);
xi6_ba1 = log(xi6_ba1);
xi7_ba1 = log(xi7_ba1);
xi8_ba1 = log(xi8_ba1);
xi9_ba1 = log(xi9_ba1);
xi10_ba1 = log(xi10_ba1);
xi11_ba1 = log(xi11_ba1);
xi12_ba1 = log(xi12_ba1);
xi13_ba1 = log(xi13_ba1);
xi14_ba1 = log(xi14_ba1);

u_ba1p =log(u_ba1p);
muz_ba1p = log(muz_ba1p); 
muupsilon_ba1p = log(muupsilon_ba1p);
qq_ba1p = log(qq_ba1p); 
r_ba1p = log(r_ba1p);
s_ba1p = log(s_ba1p); 
hd_ba1p = log(hd_ba1p);
stil_ba1p = log(stil_ba1p); 
k_ba1p = log(k_ba1p); 
iv_ba1p = log(iv_ba1p);
w_ba1p = log(w_ba1p);
output_ba1p = log(output_ba1p);
g_ba1p = log(g_ba1p);
c_ba1p = log(c_ba1p);
la_ba1p = log(la_ba1p); 
pai_ba1p = log(pai_ba1p);
f2_ba1p = log(f2_ba1p);
x2_ba1p = log(x2_ba1p);

xi1_ba1p = log(xi1_ba1p);
xi2_ba1p = log(xi2_ba1p);
xi3_ba1p = log(xi3_ba1p);
xi4_ba1p = log(xi4_ba1p);
xi5_ba1p = log(xi5_ba1p);
xi6_ba1p = log(xi6_ba1p);
xi7_ba1p = log(xi7_ba1p);
xi8_ba1p = log(xi8_ba1p);
xi9_ba1p = log(xi9_ba1p);
xi10_ba1p = log(xi10_ba1p);
xi11_ba1p = log(xi11_ba1p);
xi12_ba1p = log(xi12_ba1p);
xi13_ba1p = log(xi13_ba1p);
xi14_ba1p = log(xi14_ba1p);

u_ba2 =log(u_ba2);
muz_ba2 = log(muz_ba2); 
muupsilon_ba2 = log(muupsilon_ba2);
qq_ba2 = log(qq_ba2); 
r_ba2 = log(r_ba2);
s_ba2 = log(s_ba2); 
hd_ba2 = log(hd_ba2);
stil_ba2 = log(stil_ba2); 
k_ba2 = log(k_ba2); 
iv_ba2 = log(iv_ba2);
w_ba2 = log(w_ba2);
output_ba2 = log(output_ba2);
g_ba2 = log(g_ba2);
c_ba2 = log(c_ba2);
la_ba2 = log(la_ba2); 
pai_ba2 = log(pai_ba2);
f2_ba2 = log(f2_ba2);
x2_ba2 = log(x2_ba2);

xi1_ba2 = log(xi1_ba2);
xi2_ba2 = log(xi2_ba2);
xi3_ba2 = log(xi3_ba2);
xi4_ba2 = log(xi4_ba2);
xi5_ba2 = log(xi5_ba2);
xi6_ba2 = log(xi6_ba2);
xi7_ba2 = log(xi7_ba2);
xi8_ba2 = log(xi8_ba2);
xi9_ba2 = log(xi9_ba2);
xi10_ba2 = log(xi10_ba2);
xi11_ba2 = log(xi11_ba2);
xi12_ba2 = log(xi12_ba2);
xi13_ba2 = log(xi13_ba2);
xi14_ba2 = log(xi14_ba2);

u_ba2p =log(u_ba2p);
muz_ba2p = log(muz_ba2p); 
muupsilon_ba2p = log(muupsilon_ba2p);
qq_ba2p = log(qq_ba2p); 
r_ba2p = log(r_ba2p);
s_ba2p = log(s_ba2p); 
hd_ba2p = log(hd_ba2p);
stil_ba2p = log(stil_ba2p); 
k_ba2p = log(k_ba2p); 
iv_ba2p = log(iv_ba2p);
w_ba2p = log(w_ba2p);
output_ba2p = log(output_ba2p);
g_ba2p = log(g_ba2p);
c_ba2p = log(c_ba2p);
la_ba2p = log(la_ba2p); 
pai_ba2p = log(pai_ba2p);
f2_ba2p = log(f2_ba2p);
x2_ba2p = log(x2_ba2p);

xi1_ba2p = log(xi1_ba2p);
xi2_ba2p = log(xi2_ba2p);
xi3_ba2p = log(xi3_ba2p);
xi4_ba2p = log(xi4_ba2p);
xi5_ba2p = log(xi5_ba2p);
xi6_ba2p = log(xi6_ba2p);
xi7_ba2p = log(xi7_ba2p);
xi8_ba2p = log(xi8_ba2p);
xi9_ba2p = log(xi9_ba2p);
xi10_ba2p = log(xi10_ba2p);
xi11_ba2p = log(xi11_ba2p);
xi12_ba2p = log(xi12_ba2p);
xi13_ba2p = log(xi13_ba2p);
xi14_ba2p = log(xi14_ba2p);

u_fu1 =log(u_fu1);
muz_fu1 = log(muz_fu1); 
muupsilon_fu1 = log(muupsilon_fu1);
qq_fu1 = log(qq_fu1); 
r_fu1 = log(r_fu1);
s_fu1 = log(s_fu1); 
hd_fu1 = log(hd_fu1);
stil_fu1 = log(stil_fu1); 
k_fu1 = log(k_fu1); 
iv_fu1 = log(iv_fu1);
w_fu1 = log(w_fu1);
output_fu1 = log(output_fu1);
g_fu1 = log(g_fu1);
c_fu1 = log(c_fu1);
la_fu1 = log(la_fu1); 
pai_fu1 = log(pai_fu1);
f2_fu1 = log(f2_fu1);
x2_fu1 = log(x2_fu1);

xi1_fu1 = log(xi1_fu1);
xi2_fu1 = log(xi2_fu1);
xi3_fu1 = log(xi3_fu1);
xi4_fu1 = log(xi4_fu1);
xi5_fu1 = log(xi5_fu1);
xi6_fu1 = log(xi6_fu1);
xi7_fu1 = log(xi7_fu1);
xi8_fu1 = log(xi8_fu1);
xi9_fu1 = log(xi9_fu1);
xi10_fu1 = log(xi10_fu1);
xi11_fu1 = log(xi11_fu1);
xi12_fu1 = log(xi12_fu1);
xi13_fu1 = log(xi13_fu1);
xi14_fu1 = log(xi14_fu1);

u_fu1p =log(u_fu1p);
muz_fu1p = log(muz_fu1p); 
muupsilon_fu1p = log(muupsilon_fu1p);
qq_fu1p = log(qq_fu1p); 
r_fu1p = log(r_fu1p);
s_fu1p = log(s_fu1p); 
hd_fu1p = log(hd_fu1p);
stil_fu1p = log(stil_fu1p); 
k_fu1p = log(k_fu1p); 
iv_fu1p = log(iv_fu1p);
w_fu1p = log(w_fu1p);
output_fu1p = log(output_fu1p);
g_fu1p = log(g_fu1p);
c_fu1p = log(c_fu1p);
la_fu1p = log(la_fu1p); 
pai_fu1p = log(pai_fu1p);
f2_fu1p = log(f2_fu1p);
x2_fu1p = log(x2_fu1p);

xi1_fu1p = log(xi1_fu1p);
xi2_fu1p = log(xi2_fu1p);
xi3_fu1p = log(xi3_fu1p);
xi4_fu1p = log(xi4_fu1p);
xi5_fu1p = log(xi5_fu1p);
xi6_fu1p = log(xi6_fu1p);
xi7_fu1p = log(xi7_fu1p);
xi8_fu1p = log(xi8_fu1p);
xi9_fu1p = log(xi9_fu1p);
xi10_fu1p = log(xi10_fu1p);
xi11_fu1p = log(xi11_fu1p);
xi12_fu1p = log(xi12_fu1p);
xi13_fu1p = log(xi13_fu1p);
xi14_fu1p = log(xi14_fu1p);

approx = 1;

num_eval

[gx,hx] = gx_hx(nfy,nfx,nfyp,nfxp);

save ramsey_run.mat

function ramsey_f
%RAMSEY_f.M computes analytically differentiate the equilibrium conditions of the ramsey problem in `Optimal Inflation Stabilization in a Medium Scale Macro Model,'' by Stephanie Schmitt-Grohe and Martin Uribe.  
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe (October, 2005).
%
%Date October 2005

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

%Define symbols
allsyms

%read Ramsey FOC's
foc = ramsey_foc;

%Get rid of all variables that appear only with _ba1p but not with _ba1
foc=subs(foc,'f2_ba1p', 'f2_cu');
foc=subs(foc,'hd_ba1p', 'hd_cu');
foc=subs(foc,'qq_ba1p', 'qq_cu');
foc=subs(foc,'tauk_ba1p', 'tauk_cu');
foc=subs(foc,'u_ba1p', 'u_cu');
foc=subs(foc,'x2_ba1p', 'x2_cu');

%Get rid of all variables that appear only with _fu1 but not with _fu1p
foc=subs(foc,'c_fu1', 'c_cup');
foc=subs(foc,'hd_fu1', 'hd_cup');
foc=subs(foc,'iv_fu1', 'iv_cup');
foc=subs(foc,'k_fu1', 'k_cup');
foc=subs(foc,'output_fu1', 'output_cup');
foc=subs(foc,'qq_fu1', 'qq_cup');
foc=subs(foc,'r_fu1', 'r_cup');
foc=subs(foc,'s_fu1', 's_cup');
foc=subs(foc,'tauk_fu1', 'tauk_cup');
foc=subs(foc,'u_fu1', 'u_cup');

foc=subs(foc,'xi10_fu1', 'xi10_cup');
foc=subs(foc,'xi11_fu1', 'xi11_cup');
foc=subs(foc,'xi12_fu1', 'xi12_cup');
foc=subs(foc,'xi13_fu1', 'xi13_cup');
foc=subs(foc,'xi1_fu1', 'xi1_cup');
foc=subs(foc,'xi2_fu1', 'xi2_cup');
foc=subs(foc,'xi4_fu1', 'xi4_cup');
foc=subs(foc,'xi5_fu1', 'xi5_cup');
foc=subs(foc,'xi6_fu1', 'xi6_cup');
foc=subs(foc,'xi8_fu1', 'xi8_cup');
foc=subs(foc,'xi9_fu1', 'xi9_cup');
foc=subs(foc,'xi14_fu1', 'xi14_cup');


%get rid of variables that appear with _ba1 but with nothing higher than
%_cu
foc=subs(foc,'xi3_cu', 'xi3_ba1p');
foc=subs(foc,'xi7_cu', 'xi7_ba1p');%Note this means that we need to take xi3_cu and x7_cu out of xi_cu

%welfare at time t:

%Total hours worked
h_cu = stil_cup * hd_cu;

%Growth rate of output
muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu;


%Utility function
uf_cu = log(c_cu - B * c_ba1 / muzstar_cu)*(1-PHI4)+log(1-h_cu)*PHI4;

syms vt_cu vt_cup

%Add processes for exogenous variables

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

%Evolution of Neutral Technolgy shock
exprocmuz = log(muz_cup/MUZ) - RHOMUZ * log(muz_cu/MUZ);


%Evolution of Investment-Specific Technolgy shock
exprocmuupsilon = log(muupsilon_cup/MUUPSILON) - RHOMUUPSILON * log(muupsilon_cu/MUUPSILON);

exproc = [exprocg;exprocmuz;exprocmuupsilon];

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

[link_bwd_fwd, variables_ba1, variables_ba1p, variables_ba2, variables_ba2p, variables_fu1, variables_fu1p] = link_ba_cu_fu;

link_muz = [-muz_ba1p+muz_cu;-muz_fu1+muz_cup];

link_muupsilon = [-muupsilon_ba1p+muupsilon_cu;-muupsilon_fu1+muupsilon_cup]; 

vt = -vt_cu + uf_cu + BETTA * vt_cup;

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

statevar_cu = [variables_ba2; variables_ba1; s_cu; stil_cu;  k_cu;   muupsilon_ba1; muz_ba1; g_cu; muz_cu; muupsilon_cu]; 

statevar_cu = transpose(statevar_cu);

statevar_cup = [variables_ba2p; variables_ba1p; s_cup; stil_cup;  k_cup;  muupsilon_ba1p; muz_ba1p; g_cup; muz_cup; muupsilon_cup]; 

statevar_cup = transpose(statevar_cup);

%LAGRANGE_MULTIPLIERS.M
xi_cu= [xi1_cu  xi2_cu  xi4_cu  xi5_cu  xi6_cu  xi8_cu  xi9_cu  xi10_cu  xi11_cu  xi12_cu  xi13_cu xi14_cu ];

xi_cup=[xi1_cup xi2_cup xi4_cup xi5_cup xi6_cup xi8_cup xi9_cup xi10_cup xi11_cup xi12_cup xi13_cup xi14_cup];

controlvar_cu = [vt_cu qq_cu hd_cu output_cu c_cu la_cu pai_cu u_cu f2_cu x2_cu r_cu  iv_cu w_cu muz_fu1 muupsilon_fu1   transpose(variables_fu1) xi_cu];
 
controlvar_cup = [vt_cup qq_cup hd_cup output_cup c_cup la_cup pai_cup u_cup f2_cup x2_cup r_cup iv_cup w_cup muz_fu1p muupsilon_fu1p  transpose(variables_fu1p) xi_cup];  
  

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

%controls in logs
lc =  find(controlvar_cu~='vt_cu');

%variables to substitute from levels to logs
vs=[statevar_cu, controlvar_cu(lc), statevar_cup, controlvar_cup(lc)];

fold = f;

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

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

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 fold


%LINK_BA_CU_FU.M
function [link_bwd_fwd, variables_ba1, variables_ba1p, variables_ba2, variables_ba2p, variables_fu1, variables_fu1p] = link_ba_cu_fu;
%LINK_BA_CU_FU.M
%(c) Stephanie Schmitt-Grohe and Martin Uribe  (October 2005)


allsyms    


variables_cu = transpose([c_cu, iv_cu, la_cu, pai_cu, r_cu, w_cu, xi2_cu,   xi4_cu, xi5_cu, xi6_cu,  xi8_cu, xi9_cu]);

variables_ba1p = transpose([c_ba1p, iv_ba1p, la_ba1p, pai_ba1p, r_ba1p, w_ba1p,  xi2_ba1p,   xi4_ba1p, xi5_ba1p, xi6_ba1p, xi8_ba1p, xi9_ba1p]);


variables_ba1 = transpose([pai_ba1,  w_ba1 ]);

variables_ba2p = transpose([pai_ba2p,  w_ba2p]);


variables_fu1 = transpose([f2_fu1, la_fu1, pai_fu1, stil_fu1, w_fu1, x2_fu1]); 

variables_cup = transpose([f2_cup, la_cup, pai_cup, stil_cup, w_cup, x2_cup]); 


link_ba1p_cu = -variables_ba1p  + variables_cu;

link_ba2p_ba1 = -variables_ba2p  + variables_ba1;

link_cup_fu1 = -variables_fu1 + variables_cup;

link_bwd_fwd = [link_ba2p_ba1; link_ba1p_cu; link_cup_fu1]; 


variables_ba2 = transpose([pai_ba2,  w_ba2]);

variables_ba1 = transpose([c_ba1, iv_ba1, la_ba1, pai_ba1, r_ba1, w_ba1, xi2_ba1, xi3_ba1, xi4_ba1, xi5_ba1, xi6_ba1, xi7_ba1, xi8_ba1, xi9_ba1]);

variables_ba1p = transpose([c_ba1p, iv_ba1p, la_ba1p, pai_ba1p, r_ba1p, w_ba1p,  xi2_ba1p, xi3_ba1p,  xi4_ba1p, xi5_ba1p, xi6_ba1p, xi7_ba1p, xi8_ba1p, xi9_ba1p]);

variables_fu1p = transpose([f2_fu1p, la_fu1p, pai_fu1p, stil_fu1p, w_fu1p, x2_fu1p]); 


%RAMSEY_FOC.M
function foc = ramsey_foc;
%foc = ramsey_foc: first-order conditions of the Ramsey problem in `Optimal  Inflation Stabilization in a Medium Scale Macro Model'' by Stephanie Schmitt-Grohe and Martin Uribe.  
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%
%Date October 2005

%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
muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu;
muzstar_fu1 = muupsilon_fu1^(THETA/(1-THETA))*muz_fu1;
muzstar_ba1 = muupsilon_ba1^(THETA/(1-THETA))*muz_ba1;

objective_cu =  (1-PHI4)*log((c_cu - B * c_ba1/muzstar_cu )) +  PHI4 * log((1-stil_fu1*hd_cu));

objective_fu1 = (1-PHI4)*log((c_fu1 - B * c_cu/muzstar_fu1 )) +  PHI4 * log( (1-stil_fu1p*hd_fu1));

objective_ba1 = (1-PHI4)*log((c_ba1 - B * c_ba2/muzstar_ba1 )) +  PHI4 * log( (1-stil_cu*hd_ba1));

[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  xi9_cu  xi10_cu  xi11_cu  xi12_cu  xi13_cu  xi14_cu ];

xi_fu1 = [xi1_fu1 xi2_fu1 xi3_fu1 xi4_fu1 xi5_fu1 xi6_fu1 xi7_fu1 xi8_fu1 xi9_fu1 xi10_fu1 xi11_fu1 xi12_fu1 xi13_fu1 xi14_fu1]; 

xi_ba1 = [xi1_ba1 xi2_ba1 xi3_ba1 xi4_ba1 xi5_ba1 xi6_ba1 xi7_ba1 xi8_ba1 xi9_ba1 xi10_ba1 xi11_ba1 xi12_ba1 xi13_ba1 xi14_ba1];;

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

endogvar_cu = [u_cu qq_cu iv_cu w_cu hd_cu output_cu c_cu la_cu pai_cu f2_cu x2_cu s_cup stil_cup k_cup r_cu]; 

endogvar_cu=endogvar_cu(:);

endogvar_ba1p = [u_ba1p qq_ba1p iv_ba1p w_ba1p hd_ba1p output_ba1p c_ba1p la_ba1p pai_ba1p f2_ba1p x2_ba1p s_fu1 stil_fu1  k_fu1 r_ba1p]; 

endogvar_ba1p=endogvar_ba1p(:);

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

foc = foc(:);