%CONSTRAINTS.M
%[constraints_cu, constraints_ba1, constraints_fu1] = constraints produces symbolic expressions for the equilibrium conditions (constraints of the Ramsey problem) of the model Comparing Two 
%Variants of Calvo-Type Wage Stickiness,'' by Stephanie Schmitt-Grohe and 
%Martin Uribe (2006). 
%Suffix notation: _cu=current (or period t)
%                         _ba1=back 1 period (or period t-1)
%                         _fu1=future 1 (or period t+1)
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.

function [constraints_cu, constraints_ba1, constraints_fu1] = constraints;

allsyms

%CONSTRAINTS_CURRENT

%Define muzstar_cu, muI_cu (growth rates)
muzstar_cu = muupsilon_cu^(THETA/(1-THETA))*muz_cu;
muzstar_cup = muupsilon_cup^(THETA/(1-THETA))*muz_cup;

muI_cu = muupsilon_cu*muzstar_cu;
muI_cup = muupsilon_cup*muzstar_cup;

mula_cu = 1/muzstar_cu;%holds only if phi3=1
mula_cup = 1/muzstar_cup;

%Production Function
pf_cu =  (u_cu*k_cu/muI_cu)^THETA * hd_cu^(1-THETA) - PSSI;

%Capital utilization cost
a_cu = GAMA1 * (u_cu-1) + GAMA2/2 * (u_cu-1)^2;
a_cup = GAMA1 * (u_cup-1) + GAMA2/2 * (u_cup-1)^2;

%Capital adjustment costs
qiv_cu = iv_cu * (1 - KAPA /2 * (iv_cu/iv_ba1* muI_cu - MUI)^2);
d1qiv_cu = qiv_cu/iv_cu - KAPA*( muI_cu * iv_cu/iv_ba1-MUI )*(muI_cu * iv_cu / iv_ba1);
d2qiv_cup = (muI_cup * iv_cup/iv_cu)^2 * KAPA * (muI_cup*iv_cup / iv_cu -MUI);

%Rental rate of capital
rk_cu = diff(a_cu,'u_cu');
rk_cup = diff(a_cup,'u_cup');

%Money velocity (households)
v_cu = sqrt(PHI2/PHI1 + 1/PHI1 * (r_cu-1) / r_cu);

%Monetary transaction costs (households)
ell_cu = PHI1 * v_cu + PHI2 / v_cu - 2 * sqrt(PHI1*PHI2); %transaction cost

%Real balances (households)
mh_cu = c_cu / v_cu;

%marginal cost
mc_cu  = rk_cu /  THETA /  (u_cu*k_cu /hd_cu / muI_cu )^(THETA-1);

x1_cu  = (ETA -1 ) / ETA * x2_cu;
x1_cup = (ETA -1 ) / ETA * x2_cup;

f1_cu = f2_cu;
f1_cup = f2_cup;

%Equilibrium conditions
e1_cu = -k_cup + (1-DELTA) * k_cu/muI_cu + qiv_cu;

e2_cu = -la_cu * (1+2*PHI1*v_cu-2*(PHI1*PHI2)^(1/2)) +  (c_cu - B * c_ba1/muzstar_cu)^(-1)  - B * BETTA *  (muzstar_cup*c_cup - B * c_cu)^(-1) ;

e3_cu = -la_cu * qq_cu + BETTA * mula_cup/muupsilon_cup * la_cup * ((u_cup *rk_cup - a_cup) + qq_cup * (1-DELTA ));

e4_cu = -la_cu + la_cu * qq_cu * d1qiv_cu  + BETTA * mula_cup / muupsilon_cup * la_cup * qq_cup  * d2qiv_cup;  

e5_cu = -f1_cu + (ETATIL-1) / ETATIL * wtil_cu * la_cu * (w_cu/wtil_cu)^ETATIL * hd_cu + ALFATIL * BETTA * (pai_cup / (MUZSTAR*pai_cu)^CHITIL)^(ETATIL-1) * (muzstar_cup*wtil_cup/wtil_cu)^(ETATIL-1) * mula_cup * muzstar_cup * f1_cup; 

e7_cu = -la_cu + BETTA * r_cu * mula_cup *  la_cup / pai_cup;

e8_cu = -x1_cu + output_cu * mc_cu * ptil_cu^(-ETA-1) + ALFA * BETTA * mula_cup * la_cup / la_cu * (ptil_cu / ptil_cup)^(-ETA-1) * (pai_cu^CHI / pai_cup)^(-ETA) * muzstar_cup *  x1_cup;

e9_cu = -x2_cu + output_cu * ptil_cu^(-ETA) + ALFA * BETTA * mula_cup * la_cup / la_cu * (pai_cu^CHI / pai_cup )^(1-ETA) * (ptil_cu / ptil_cup)^(-ETA) * muzstar_cup * x2_cup;

e10_cu = -pf_cu + output_cu * s_cup;

e11_cu = - s_cup + (1 - ALFA) * ptil_cu^(-ETA) + ALFA * (pai_cu/ pai_ba1^CHI)^ETA * s_cu;

%Investment
e14_cu = -iv_cu + output_cu  - c_cu  * (1+ell_cu) - g_cu - a_cu * k_cu / muI_cu;

%real wage
e15_cu = -w_cu + mc_cu * (1-THETA) * (u_cu * k_cu / muI_cu  / hd_cu )^(THETA)   / (1 + NU * (r_cu -1)/r_cu );

%ptil
e16_cu = - ptil_cu  + ((1 - ALFA * pai_cu ^(ETA-1) * pai_ba1 ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

%wtil
e17_cu = - wtil_cu  + ((w_cu ^(1-ETATIL) - ALFATIL * (w_ba1/muzstar_cu)^(1-ETATIL) * ((MUZSTAR*pai_ba1)^CHITIL / pai_cu) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

%f2 for the EHL model!
e6_ehl_cu = [-f2_cu + PHIL *(wtil_cu/w_cu)^(-ETATIL*(1+PHI5)) * hd_cu^(1+PHI5) + ALFATIL * BETTA * (pai_cup/(MUZSTAR * pai_cu)^(CHITIL)*muzstar_cup * wtil_cup/wtil_cu)^(ETATIL*(1+PHI5)) * f2_cup;
             -stil_cup^(1+PHI5) + (1-ALFATIL) *(wtil_cu / w_cu)^(-ETATIL*(1+PHI5)) + ALFATIL * (w_ba1 *(MUZSTAR*pai_ba1)^CHITIL / pai_cu /muzstar_cu / w_cu)^(-ETATIL*(1+PHI5))*stil_cu^(1+PHI5)];


%f2 for the SGU model: 
e6_sgu_cu = [-f2_cu + PHIL *stil_cup^PHI5*(wtil_cu/w_cu)^(-ETATIL) * hd_cu^(1+PHI5) + ALFATIL * BETTA * (pai_cup/(MUZSTAR * pai_cu)^(CHITIL)*muzstar_cup * wtil_cup/wtil_cu)^(ETATIL) * f2_cup;
             -stil_cup + (1-ALFATIL) *(wtil_cu / w_cu)^(-ETATIL) + ALFATIL * (w_ba1 *(MUZSTAR*pai_ba1)^CHITIL / pai_cu /muzstar_cu / w_cu)^(-ETATIL)*stil_cu];

%Create vector of current constraints
constraints_cu = [e1_cu; e2_cu; e3_cu; e4_cu; e5_cu;  e7_cu; -e8_cu;  e9_cu; -e10_cu; -e11_cu;  e14_cu; e15_cu; e16_cu; e17_cu;-e6_ehl_cu;-e6_sgu_cu];

%CONSTRAINTS_ba1RRENT

%Define muzstar_ba1, muI_ba1
muzstar_ba1 = muupsilon_ba1^(THETA/(1-THETA))*muz_ba1;
muzstar_ba1p = muupsilon_ba1p^(THETA/(1-THETA))*muz_ba1p;

muI_ba1 = muupsilon_ba1*muzstar_ba1;
muI_ba1p = muupsilon_ba1p*muzstar_ba1p;

mula_ba1 = 1/muzstar_ba1;%holds only if phi3=1
mula_ba1p = 1/muzstar_ba1p;

%Production Function
pf_ba1 =  (u_ba1*k_ba1/muI_ba1)^THETA * hd_ba1^(1-THETA) - PSSI;

%Capital utilization cost
a_ba1 = GAMA1 * (u_ba1-1) + GAMA2/2 * (u_ba1-1)^2;
a_ba1p = GAMA1 * (u_ba1p-1) + GAMA2/2 * (u_ba1p-1)^2;

%Capital adjustment costs
qiv_ba1 = iv_ba1 * (1 - KAPA /2 * (iv_ba1/iv_ba2* muI_ba1 - MUI)^2);
d1qiv_ba1 = qiv_ba1/iv_ba1 - KAPA*( muI_ba1 * iv_ba1/iv_ba2-MUI )*(muI_ba1 * iv_ba1 / iv_ba2);
d2qiv_ba1p = (muI_ba1p * iv_ba1p/iv_ba1)^2 * KAPA * (muI_ba1p*iv_ba1p / iv_ba1 -MUI);

rk_ba1 = diff(a_ba1,'u_ba1');
rk_ba1p = diff(a_ba1p,'u_ba1p');

v_ba1 = sqrt(PHI2/PHI1 + 1/PHI1 * (r_ba1-1) / r_ba1);

ell_ba1 = PHI1 * v_ba1 + PHI2 / v_ba1 - 2 * sqrt(PHI1*PHI2); %transaction cost

mh_ba1 = c_ba1 / v_ba1;

mc_ba1  = rk_ba1 /  THETA /  (u_ba1*k_ba1 /hd_ba1 / muI_ba1 )^(THETA-1);

x1_ba1  = (ETA -1 ) / ETA * x2_ba1;
x1_ba1p = (ETA -1 ) / ETA * x2_ba1p;

f1_ba1 = f2_ba1;
f1_ba1p = f2_ba1p;

e1_ba1 = -k_ba1p + (1-DELTA) * k_ba1/muI_ba1 + qiv_ba1;

e2_ba1 = -la_ba1 * (1+2*PHI1*v_ba1-2*(PHI1*PHI2)^(1/2)) +  (c_ba1 - B * c_ba2/muzstar_ba1)^(-1)  - B * BETTA *  (muzstar_ba1p*c_ba1p - B * c_ba1)^(-1) ;

e3_ba1 = -la_ba1 * qq_ba1 + BETTA * mula_ba1p/muupsilon_ba1p * la_ba1p * ((u_ba1p *rk_ba1p - a_ba1p) + qq_ba1p * (1-DELTA ));

e4_ba1 = -la_ba1 + la_ba1 * qq_ba1 * d1qiv_ba1  + BETTA * mula_ba1p / muupsilon_ba1p * la_ba1p * qq_ba1p  * d2qiv_ba1p;  

e5_ba1 = -f1_ba1 + (ETATIL-1) / ETATIL * wtil_ba1 * la_ba1 * (w_ba1/wtil_ba1)^ETATIL * hd_ba1 + ALFATIL * BETTA * (pai_ba1p / (MUZSTAR*pai_ba1)^CHITIL)^(ETATIL-1) * (muzstar_ba1p*wtil_ba1p/wtil_ba1)^(ETATIL-1) * mula_ba1p * muzstar_ba1p * f1_ba1p; 

e7_ba1 = -la_ba1 + BETTA * r_ba1 * mula_ba1p *  la_ba1p / pai_ba1p;

e8_ba1 = -x1_ba1 + output_ba1 * mc_ba1 * ptil_ba1^(-ETA-1) + ALFA * BETTA * mula_ba1p * la_ba1p / la_ba1 * (ptil_ba1 / ptil_ba1p)^(-ETA-1) * (pai_ba1^CHI / pai_ba1p)^(-ETA) * muzstar_ba1p *  x1_ba1p;

e9_ba1 = -x2_ba1 + output_ba1 * ptil_ba1^(-ETA) + ALFA * BETTA * mula_ba1p * la_ba1p / la_ba1 * (pai_ba1^CHI / pai_ba1p )^(1-ETA) * (ptil_ba1 / ptil_ba1p)^(-ETA) * muzstar_ba1p * x2_ba1p;

e10_ba1 = -pf_ba1 + output_ba1 * s_ba1p;

e11_ba1 = - s_ba1p + (1 - ALFA) * ptil_ba1^(-ETA) + ALFA * (pai_ba1/ pai_ba2^CHI)^ETA * s_ba1;

%Investment
e14_ba1 = -iv_ba1 + output_ba1  - c_ba1  * (1+ell_ba1) - g_ba1 - a_ba1 * k_ba1 / muI_ba1;

%real wage
e15_ba1 = -w_ba1 + mc_ba1 * (1-THETA) * (u_ba1 * k_ba1 / muI_ba1  / hd_ba1 )^(THETA)   / (1 + NU * (r_ba1 -1)/r_ba1 );

%ptil
e16_ba1 = - ptil_ba1  + ((1 - ALFA * pai_ba1 ^(ETA-1) * pai_ba2 ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

%wtil
e17_ba1 = - wtil_ba1  + ((w_ba1 ^(1-ETATIL) - ALFATIL * (w_ba2/muzstar_ba1)^(1-ETATIL) * ((MUZSTAR*pai_ba2)^CHITIL / pai_ba1) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

%f2 for the EHL model!
e6_ehl_ba1 = [-f2_ba1 + PHIL *(wtil_ba1/w_ba1)^(-ETATIL*(1+PHI5)) * hd_ba1^(1+PHI5) + ALFATIL * BETTA * (pai_ba1p/(MUZSTAR * pai_ba1)^(CHITIL)*muzstar_ba1p * wtil_ba1p/wtil_ba1)^(ETATIL*(1+PHI5)) * f2_ba1p;
             -stil_ba1p^(1+PHI5) + (1-ALFATIL) *(wtil_ba1 / w_ba1)^(-ETATIL*(1+PHI5)) + ALFATIL * (w_ba2 *(MUZSTAR*pai_ba2)^CHITIL / pai_ba1 /muzstar_ba1 / w_ba1)^(-ETATIL*(1+PHI5))*stil_ba1^(1+PHI5)];
                        

%f2 for the SGU model: 
e6_sgu_ba1 = [-f2_ba1 + PHIL *stil_ba1p^PHI5*(wtil_ba1/w_ba1)^(-ETATIL) * hd_ba1^(1+PHI5) + ALFATIL * BETTA * (pai_ba1p/(MUZSTAR * pai_ba1)^(CHITIL)*muzstar_ba1p * wtil_ba1p/wtil_ba1)^(ETATIL) * f2_ba1p;
             -stil_ba1p + (1-ALFATIL) *(wtil_ba1 / w_ba1)^(-ETATIL) + ALFATIL * (w_ba2 *(MUZSTAR*pai_ba2)^CHITIL / pai_ba1 /muzstar_ba1 / w_ba1)^(-ETATIL)*stil_ba1];

%Create vector of current constraints
constraints_ba1 = [e1_ba1; e2_ba1; e3_ba1; e4_ba1; e5_ba1;  e7_ba1; -e8_ba1;  e9_ba1; -e10_ba1; -e11_ba1;  e14_ba1; e15_ba1; e16_ba1; e17_ba1;-e6_ehl_ba1;-e6_sgu_ba1];


%CONSTRAINTS_fu1RRENT

%Define muzstar_fu1, muI_fu1
muzstar_fu1 = muupsilon_fu1^(THETA/(1-THETA))*muz_fu1;
muzstar_fu1p = muupsilon_fu1p^(THETA/(1-THETA))*muz_fu1p;

muI_fu1 = muupsilon_fu1*muzstar_fu1;
muI_fu1p = muupsilon_fu1p*muzstar_fu1p;

mula_fu1 = 1/muzstar_fu1;%holds only if phi3=1
mula_fu1p = 1/muzstar_fu1p;

%Production Function
pf_fu1 =  (u_fu1*k_fu1/muI_fu1)^THETA * hd_fu1^(1-THETA) - PSSI;

%Capital utilization cost
a_fu1 = GAMA1 * (u_fu1-1) + GAMA2/2 * (u_fu1-1)^2;
a_fu1p = GAMA1 * (u_fu1p-1) + GAMA2/2 * (u_fu1p-1)^2;

%Capital adjustment costs
qiv_fu1 = iv_fu1 * (1 - KAPA /2 * (iv_fu1/iv_cu* muI_fu1 - MUI)^2);
d1qiv_fu1 = qiv_fu1/iv_fu1 - KAPA*( muI_fu1 * iv_fu1/iv_cu-MUI )*(muI_fu1 * iv_fu1 / iv_cu);
d2qiv_fu1p = (muI_fu1p * iv_fu1p/iv_fu1)^2 * KAPA * (muI_fu1p*iv_fu1p / iv_fu1 -MUI);

rk_fu1 = diff(a_fu1,'u_fu1');
rk_fu1p = diff(a_fu1p,'u_fu1p');

v_fu1 = sqrt(PHI2/PHI1 + 1/PHI1 * (r_fu1-1) / r_fu1);

ell_fu1 = PHI1 * v_fu1 + PHI2 / v_fu1 - 2 * sqrt(PHI1*PHI2); %transaction cost

mh_fu1 = c_fu1 / v_fu1;

mc_fu1  = rk_fu1 /  THETA /  (u_fu1*k_fu1 /hd_fu1 / muI_fu1 )^(THETA-1);

x1_fu1  = (ETA -1 ) / ETA * x2_fu1;
x1_fu1p = (ETA -1 ) / ETA * x2_fu1p;

f1_fu1 = f2_fu1;
f1_fu1p = f2_fu1p;

e1_fu1 = -k_fu1p + (1-DELTA) * k_fu1/muI_fu1 + qiv_fu1;

e2_fu1 = -la_fu1 * (1+2*PHI1*v_fu1-2*(PHI1*PHI2)^(1/2)) +  (c_fu1 - B * c_cu/muzstar_fu1)^(-1)  - B * BETTA *  (muzstar_fu1p*c_fu1p - B * c_fu1)^(-1) ;

e3_fu1 = -la_fu1 * qq_fu1 + BETTA * mula_fu1p/muupsilon_fu1p * la_fu1p * ((u_fu1p *rk_fu1p - a_fu1p) + qq_fu1p * (1-DELTA ));

e4_fu1 = -la_fu1 + la_fu1 * qq_fu1 * d1qiv_fu1  + BETTA * mula_fu1p / muupsilon_fu1p * la_fu1p * qq_fu1p  * d2qiv_fu1p;  

e5_fu1 = -f1_fu1 + (ETATIL-1) / ETATIL * wtil_fu1 * la_fu1 * (w_fu1/wtil_fu1)^ETATIL * hd_fu1 + ALFATIL * BETTA * (pai_fu1p / (MUZSTAR*pai_fu1)^CHITIL)^(ETATIL-1) * (muzstar_fu1p*wtil_fu1p/wtil_fu1)^(ETATIL-1) * mula_fu1p * muzstar_fu1p * f1_fu1p; 

e7_fu1 = -la_fu1 + BETTA * r_fu1 * mula_fu1p *  la_fu1p / pai_fu1p;

e8_fu1 = -x1_fu1 + output_fu1 * mc_fu1 * ptil_fu1^(-ETA-1) + ALFA * BETTA * mula_fu1p * la_fu1p / la_fu1 * (ptil_fu1 / ptil_fu1p)^(-ETA-1) * (pai_fu1^CHI / pai_fu1p)^(-ETA) * muzstar_fu1p *  x1_fu1p;

e9_fu1 = -x2_fu1 + output_fu1 * ptil_fu1^(-ETA) + ALFA * BETTA * mula_fu1p * la_fu1p / la_fu1 * (pai_fu1^CHI / pai_fu1p )^(1-ETA) * (ptil_fu1 / ptil_fu1p)^(-ETA) * muzstar_fu1p * x2_fu1p;

e10_fu1 = -pf_fu1 + output_fu1 * s_fu1p;

e11_fu1 = - s_fu1p + (1 - ALFA) * ptil_fu1^(-ETA) + ALFA * (pai_fu1/ pai_cu^CHI)^ETA * s_fu1;

%Investment
e14_fu1 = -iv_fu1 + output_fu1  - c_fu1  * (1+ell_fu1) - g_fu1 - a_fu1 * k_fu1 / muI_fu1;

%real wage
e15_fu1 = -w_fu1 + mc_fu1 * (1-THETA) * (u_fu1 * k_fu1 / muI_fu1  / hd_fu1 )^(THETA)   / (1 + NU * (r_fu1 -1)/r_fu1 );

%ptil
e16_fu1 = - ptil_fu1  + ((1 - ALFA * pai_fu1 ^(ETA-1) * pai_cu ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

%wtil
e17_fu1 = - wtil_fu1  + ((w_fu1 ^(1-ETATIL) - ALFATIL * (w_cu/muzstar_fu1)^(1-ETATIL) * ((MUZSTAR*pai_cu)^CHITIL / pai_fu1) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

%f2 for the EHL model!
e6_ehl_fu1 = [-f2_fu1 + PHIL *(wtil_fu1/w_fu1)^(-ETATIL*(1+PHI5)) * hd_fu1^(1+PHI5) + ALFATIL * BETTA * (pai_fu1p/(MUZSTAR * pai_fu1)^(CHITIL)*muzstar_fu1p * wtil_fu1p/wtil_fu1)^(ETATIL*(1+PHI5)) * f2_fu1p;
             -stil_fu1p^(1+PHI5) + (1-ALFATIL) *(wtil_fu1 / w_fu1)^(-ETATIL*(1+PHI5)) + ALFATIL * (w_cu *(MUZSTAR*pai_cu)^CHITIL / pai_fu1 /muzstar_fu1 / w_fu1)^(-ETATIL*(1+PHI5))*stil_fu1^(1+PHI5)];

%f2 for the SGU model: 
e6_sgu_fu1 = [-f2_fu1 + PHIL *stil_fu1p^PHI5*(wtil_fu1/w_fu1)^(-ETATIL) * hd_fu1^(1+PHI5) + ALFATIL * BETTA * (pai_fu1p/(MUZSTAR * pai_fu1)^(CHITIL)*muzstar_fu1p * wtil_fu1p/wtil_fu1)^(ETATIL) * f2_fu1p;
             -stil_fu1p + (1-ALFATIL) *(wtil_fu1 / w_fu1)^(-ETATIL) + ALFATIL * (w_cu *(MUZSTAR*pai_cu)^CHITIL / pai_fu1 /muzstar_fu1 / w_fu1)^(-ETATIL)*stil_fu1];

%Create vector of current constraints
constraints_fu1 = [e1_fu1; e2_fu1; e3_fu1; e4_fu1; e5_fu1;  e7_fu1; -e8_fu1;  e9_fu1; -e10_fu1; -e11_fu1;  e14_fu1; e15_fu1; e16_fu1; e17_fu1;-e6_ehl_fu1;-e6_sgu_fu1];