%RAMSEY_RUN.M
%[gx, hx]=ramsey_run(model) computes 
%a first-order approximation to the Ramsey equilibrium conditions of the model developed in ``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by Stephanie Schmitt-Grohe and Martin Uribe (2006).
%[gx, hx, controlvar_cu, statevar_cu, ETAMATRIX] = ramsey_run(model) computes in addition symbolic vectors of control and state variables and a numeric matrix  (ETAMATRIX) that premultiplies the vector of innovations to the state vector.  
%[gx, hx, controlvar_cu, statevar_cu, ETAMATRIX, gxx, hxx, gss, hss, conditional_welfare, unconditional_welfare] = ramsey_run(model,2) computes 
%the second-order approximation to the Ramsey policy functions.
%[gx, hx, gxx, hxx, gss, hss, conditional_welfare, unconditional_welfare]=ramsey_run(model,2) computes a second-order approximation to conditional and unconditional welfare in the Ramsey equilibrium.
%inputs: 
%model is a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%The default value is 'sgu'.
%approx takes the values 1 for a linear approximation or 2 for a second-order approximation. The default value is 1. 
%Ouput is saved in a mat file called ramsey_run_ehl.mat  or ramsey_run_sgu.mat depending on the value of the variable model
%Subfunctions: ramsey_f.m, ramsey_foc.m, link_ba_cu_fu.m 
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.
function [gx, hx, controlvar_cu, statevar_cu, ETAMATRIX, gxx, hxx, gss, hss, conditional_welfare, unconditional_welfare] = ramsey_run(model,approx);

if nargin==0
 model='sgu';
 approx = 1;
elseif nargin == 1
 approx=1;
end

ramsey_f(model);

eval(['load ramsey_f_' model])

ramsey_ss(model)
eval(['load ramsey_ss_' model]) %This is saved automatically, when running ramsey_ss.m

%make all variables into logs (but not vt_cu) 
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);
ptil_cu = log(ptil_cu);
wtil_cu = log(wtil_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);
xi15_cu = log(xi15_cu);
xi16_cu = log(xi16_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);
ptil_cup = log(ptil_cup);
wtil_cup = log(wtil_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);
xi15_cup = log(xi15_cup);
xi16_cup = log(xi16_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);
ptil_ba1 = log(ptil_ba1);
wtil_ba1 = log(wtil_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);
xi15_ba1 = log(xi15_ba1);
xi16_ba1 = log(xi16_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);
ptil_ba1p = log(ptil_ba1p);
wtil_ba1p = log(wtil_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);
xi15_ba1p = log(xi15_ba1p);
xi16_ba1p= log(xi16_ba1p);

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);
ptil_fu1 = log(ptil_fu1);
wtil_fu1 = log(wtil_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);
xi15_fu1 = log(xi15_fu1);
xi16_fu1 = log(xi16_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);
ptil_fu1p = log(ptil_fu1p);
wtil_fu1p = log(wtil_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);
xi15_fu1p = log(xi15_fu1p);
xi16_fu1p = log(xi16_fu1p);


%numerial evaluation of function f and its derivatives
num_eval 

%Compute firt-order approximation to policy functions
[gx,hx,exitflag] = gx_hx(nfy,nfx,nfyp,nfxp);

%Matrix Scaling Standard Deviations of Innovation to State Vector

ne=3; 
ETAMATRIX = zeros(size(fx,2),ne); 
ETAMATRIX(end-2:end,:) = [STD_EPSG 0 0; 0 STD_EPSMUZ 0; 0 0 STD_EPSMUUPSILON];


%Second-order approximation to policy function
if approx==2

[gxx,hxx] = gxx_hxx_sparse(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx);

[gss,hss] = gss_hss(nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx,hx,gx,gxx,ETAMATRIX);

%Scalar scaling the the standard deviations of exogenous shocks
sig =1;

%Compute welfare conditional on the initial state being the ramsey deterministic steady state
conditional_welfare = vt_cu  +  gss(1)  * sig^2 / 2

%Unconditional expectations
[Ey,Ex] = unconditional_mean(gx, hx, gxx, hxx, gss, hss, ETAMATRIX, 1);

%Unconditional expectation of welfare 
unconditional_welfare = Ey(1) + vt_cu; 

else %if approx==2
 gxx = [];
 hxx = [];
 gss = [];
 hss = [];
 conditional_welfare = [];
 unconditional_welfare = [];
end %if approx==2


eval(['save ramsey_run_' model '.mat'])



%RAMSEY_f.M
%ramsey_f(model) computes symbolically the equilibrium conditions of the ramsey problem in ``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by Stephanie Schmitt-Grohe and Martin Uribe (2006)  
%Input:
%model is a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%The default value is 'sgu'.
%(c) Stephanie Schmitt-Grohe and Martin Uribe (October, 2006).


function ramsey_f(model)

if nargin < 1
    model = 'sgu'
end

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

%Define symbols
allsyms

%read Ramsey FOC's
foc = ramsey_foc(model);

%The following algebraic operations are performed to reduced the size of the system

%Get rid of all variables that appear only with _ba1p but not with _ba1
foc=subs(foc,'f2_ba1p', 'f2_cu', 0);
foc=subs(foc,'muupsilon_ba1p', 'muupsilon_cu', 0);
foc=subs(foc,'muz_ba1p', 'muz_cu', 0);
foc=subs(foc,'qq_ba1p', 'qq_cu', 0);
foc=subs(foc,'u_ba1p', 'u_cu', 0);
foc=subs(foc,'x2_ba1p', 'x2_cu', 0);

%Get rid of all variables that appear only with _fu1 but not with _fu1p
foc=subs(foc,'c_fu1', 'c_cup', 0);
foc=subs(foc,'hd_fu1', 'hd_cup', 0);
foc=subs(foc,'iv_fu1', 'iv_cup', 0);
foc=subs(foc,'k_fu1', 'k_cup', 0);
foc=subs(foc,'la_fu1', 'la_cup', 0);
foc=subs(foc,'muupsilon_fu1', 'muupsilon_cup', 0);
foc=subs(foc,'muz_fu1', 'muz_cup', 0);
foc=subs(foc,'output_fu1', 'output_cup', 0);
foc=subs(foc,'pai_fu1', 'pai_cup', 0);
foc=subs(foc,'ptil_fu1', 'ptil_cup', 0);
foc=subs(foc,'qq_fu1', 'qq_cup', 0);
foc=subs(foc,'r_fu1', 'r_cup', 0);
foc=subs(foc,'s_fu1', 's_cup', 0);
foc=subs(foc,'stil_fu1', 'stil_cup', 0);
foc=subs(foc,'u_fu1', 'u_cup', 0);
foc=subs(foc,'w_fu1', 'w_cup', 0);
foc=subs(foc,'wtil_fu1', 'wtil_cup', 0);
foc=subs(foc,'xi10_fu1', 'xi10_cup', 0);
foc=subs(foc,'xi11_fu1', 'xi11_cup', 0);
foc=subs(foc,'xi12_fu1', 'xi12_cup', 0);
foc=subs(foc,'xi13_fu1', 'xi13_cup', 0);
foc=subs(foc,'xi14_fu1', 'xi14_cup', 0);
foc=subs(foc,'xi16_fu1', 'xi16_cup', 0);
foc=subs(foc,'xi1_fu1', 'xi1_cup', 0);
foc=subs(foc,'xi2_fu1', 'xi2_cup', 0);
foc=subs(foc,'xi4_fu1', 'xi4_cup', 0);
foc=subs(foc,'xi7_fu1', 'xi7_cup', 0);
foc=subs(foc,'xi9_fu1', 'xi9_cup', 0);

%get rid of variables that appear with _ba1 but with nothing higher than
%_cu
foc=subs(foc,'xi15_cu', 'xi15_ba1p', 0);
foc=subs(foc,'xi3_cu', 'xi3_ba1p', 0);
foc=subs(foc,'xi5_cu', 'xi5_ba1p', 0);
foc=subs(foc,'xi6_cu', 'xi6_ba1p', 0);
foc=subs(foc,'xi8_cu', 'xi8_ba1p', 0);%Note this means that we need to take xi15_cu, xi3_cu, xi5_cu, xi6_cu, and xi8_cu out of xi_cu

%welfare at time t:

%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)- PHIL * (stil_cup*hd_cu)^(1+PHI5)/(1+PHI5);


%Declare new symbolic varibles
syms vt_cu vt_cup

%Expceted value of welfare conditional on information available at time t
vt = -vt_cu + uf_cu + BETTA * 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] = link_ba_cu_fu;

%Full set of Rasmey equilibrium conditions inclusing conditional welfare
f = [vt; foc;  exproc; link_bwd_fwd];

%Vector of states
statevar_cu =  [variables_ba1;  k_cu;  s_cu;  stil_cu;  g_cu;  muz_cu;  muupsilon_cu]; 

statevar_cu = transpose(statevar_cu);

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

statevar_cup = transpose(statevar_cup);

%LAGRANGE_MULTIPLIERS.M
xi_cu= [xi1_cu  xi2_cu  xi4_cu  xi7_cu  xi9_cu xi10_cu  xi11_cu  xi12_cu  xi13_cu xi14_cu  xi16_cu];%xi3_cu, xi5_cu, xi6_cu, xi8_cu, and xi15_cu are excluded. They appear only in _ba1 and are included as such in variables_ba1

xi_cup=[xi1_cup xi2_cup xi4_cup xi7_cup xi9_cup xi10_cup xi11_cup xi12_cup xi13_cup xi14_cup xi16_cup];


%Vector of control variables 
controlvar_cu = [vt_cu  c_cu  iv_cu  la_cu pai_cu ptil_cu r_cu w_cu  wtil_cu f2_cu hd_cu output_cu qq_cu u_cu  x2_cu xi_cu];
 
controlvar_cup = [vt_cup c_cup  iv_cup  la_cup pai_cup ptil_cup r_cup w_cup  wtil_cup f2_cup hd_cup output_cup qq_cup u_cup x2_cup 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
eval(['save ramsey_f_' model '.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'])


%RAMSEY_FOC.M
%foc = ramsey_foc(model) computes analyticallythe  first-order conditions of the 
%Ramsey problem in ``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by Stephanie 
%Schmitt-Grohe and Martin Uribe (2006).
%Input:
%model = a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%The default value is 'sgu'.
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 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}

function foc = ramsey_foc(model);

if nargin<1
    model = 'sgu'
end

%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 =  log(c_cu  - B * c_ba1/muzstar_cu ) -  PHIL * (stil_fu1*hd_cu  )^(1+PHI5)/(1+PHI5);
objective_fu1 = log(c_fu1 - B * c_cu /muzstar_fu1) -  PHIL * (stil_fu1p*hd_fu1)^(1+PHI5)/(1+PHI5);
objective_ba1 = log(c_ba1 - B * c_ba2/muzstar_ba1) -  PHIL * (stil_cu*hd_ba1  )^(1+PHI5)/(1+PHI5);

[constraints_cu, constraints_ba1, constraints_fu1] = constraints;

%if you want to use the EHL model
if strcmp(model,'ehl')
constraints_cu=constraints_cu(1:end-2);
constraints_ba1=constraints_ba1(1:end-2);
constraints_fu1=constraints_fu1(1:end-2);


else
%if you want the SGU model
constraints_cu=constraints_cu([1:end-4 end-1:end]);
constraints_ba1=constraints_ba1([1:end-4 end-1:end]);
constraints_fu1=constraints_fu1([1:end-4 end-1:end]);
end

%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  xi15_cu  xi16_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 xi15_fu1 xi16_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 xi15_ba1 xi16_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 = [ptil_cu wtil_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 = [ptil_ba1p wtil_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(:);



%LINK_BA_CU_FU.M
%[link_bwd_fwd, variables_ba1, variables_ba1p ] = link_ba_cu_fu produces symbolic auxiliary equations and variables useful in trasforming a system of difference equations of order higher than first into a a first-order system
function [link_bwd_fwd, variables_ba1, variables_ba1p ] = link_ba_cu_fu;

allsyms    


variables_cu = transpose([c_cu, iv_cu, la_cu, pai_cu, ptil_cu, r_cu, w_cu,  wtil_cu, xi2_cu, xi4_cu, xi7_cu ]);

variables_ba1p = transpose([c_ba1p, iv_ba1p, la_ba1p, pai_ba1p, ptil_ba1p, r_ba1p, w_ba1p,  wtil_ba1p, xi2_ba1p, xi4_ba1p, xi7_ba1p ]);


link_ba1p_cu = - variables_ba1p  + variables_cu;


link_bwd_fwd = [link_ba1p_cu ]; 



%Note: in variables_ba1 we also include those variables that appear in ba_1
%only: xi3_ba1, xi5_ba1, xi6_ba1, and xi8_ba1; and xi15_ba1

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

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