%CONSTRAINTS.M
function [objective_cu, constraints_cu, objective_ba1, constraints_ba1, objective_fu1, constraints_fu1] = constraints;

allsyms

%CONSTRAINTS_CURRENT

%Production Function
pf_cu = z_cu * (u_cu*k_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 -1)^2);

d1qiv_cu = qiv_cu/iv_cu -KAPA*(iv_cu/iv_ba1-1)*iv_cu/iv_ba1;
d2qiv_cup = (iv_cup/iv_cu)^2*KAPA*(iv_cup / iv_cu -1);

rk_cu = diff(a_cu,'u_cu');
rk_cup = diff(a_cup,'u_cup');


v_cu = sqrt(PHI2/PHI1 + 1/PHI1 * (r_cu-1) / r_cu);

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

mh_cu = c_cu / v_cu;

mc_cu  = rk_cu / (z_cu  * THETA * (u_cu*k_cu /hd_cu )^(THETA-1));

ptil_cu  = ((1 - ALFA * pai_cu ^(ETA-1) * pai_ba1 ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

ptil_cup = ((1 - ALFA * pai_cup^(ETA-1) * pai_ba1p^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

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

wtil_cu  = ((w_cu ^(1-ETATIL) - ALFATIL * w_ba1 ^(1-ETATIL) * (pai_ba1^CHITIL / pai_cu) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

wtil_cup = ((w_cup^(1-ETATIL) - ALFATIL * w_ba1p ^(1-ETATIL) * (pai_ba1p^CHITIL / pai_cup) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

h_cu = stil_cup * hd_cu;
h_cup = stil_fu1p * hd_cup;


f1_cu = f2_cu;
f1_cup = f2_cup;

profits_cu = output_cu - rk_cu  * u_cu * k_cu - w_cu * hd_cu - NU * (1-1/r_cu) * w_cu * hd_cu;

m_cu = mh_cu + NU * w_cu * hd_cu;

tauprofits_cu = tauk_cu;%TAUPROFITS;%

tau_cu = tauk_cu * (u_cu *rk_cu - a_cu - qq_cu * DELTA) * k_cu + tauh_cu * hd_cu * w_cu + tauprofits_cu * profits_cu;

%Utility function
uf_cu = (((c_cu - B * c_ba1)^(1-PHI4) * (1-h_cu)^PHI4)^(1-PHI3)-1) / (1-PHI3);

objective_cu = uf_cu;

e1_cu = -k_cup + (1-DELTA) * k_cu + qiv_cu;

e2_cu = -la_cu * (1+2*PHI1*v_cu-2*(PHI1*PHI2)^(1/2)) + (1-PHI4) * (c_cu - B * c_ba1)^((1-PHI4)*(1-PHI3)-1) * (1-h_cu)^(PHI4*(1-PHI3)) - B * BETTA * (1-PHI4) * (c_cup - B * c_cu)^((1-PHI4)*(1-PHI3)-1) * (1-h_cup)^(PHI4*(1-PHI3));

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

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

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

e6_cu = -f2_cu  + PHI4 * (c_cu-B*c_ba1)^((1-PHI3)*(1-PHI4)) * (1-h_cu)^(PHI4*(1-PHI3)-1) *  (w_cu/wtil_cu)^ETATIL * hd_cu + ALFATIL * BETTA * (pai_cup / pai_cu^CHITIL)^ETATIL * f2_cup* (wtil_cup/wtil_cu)^ETATIL;

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

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

e9_cu = -x2_cu + output_cu * ptil_cu^(-ETA) + ALFA * BETTA * la_cup / la_cu * (pai_cu^CHI / pai_cup )^(1-ETA) * (ptil_cu / ptil_cup)^(-ETA) * 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;

e12_cu = -stil_cup + (1-ALFATIL) * (wtil_cu/w_cu)^(-ETATIL) + ALFATIL * (w_ba1 / w_cu) ^(-ETATIL) * (pai_cu / pai_ba1^CHITIL)^ETATIL * stil_cu;

e13_cu = -assets_cup/r_cu - (1-1/r_cu) * m_cu - tau_cu + assets_cu / pai_cu + g_cu + n_cu;

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

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

 %Create vector of current constraints
constraints_cu = [e1_cu; e2_cu; e3_cu; e4_cu; e5_cu; -e6_cu;  e7_cu; -e8_cu;  e9_cu; -e10_cu; -e11_cu; -e12_cu; -e13_cu; e14_cu; e15_cu];

%CONSTRAINTS_ba1RRENT

%Production Function
pf_ba1 = z_ba1 * (u_ba1*k_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 -1)^2);

d1qiv_ba1 = qiv_ba1/iv_ba1 -KAPA*(iv_ba1/iv_ba2-1)*iv_ba1/iv_ba2;
d2qiv_ba1p = (iv_ba1p/iv_ba1)^2*KAPA*(iv_ba1p / iv_ba1 -1);

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 / (z_ba1  * THETA * (u_ba1*k_ba1 /hd_ba1 )^(THETA-1));

ptil_ba1  = ((1 - ALFA * pai_ba1 ^(ETA-1) * pai_ba2 ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

ptil_ba1p = ((1 - ALFA * pai_ba1p^(ETA-1) * pai_ba2p^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

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

wtil_ba1  = ((w_ba1 ^(1-ETATIL) - ALFATIL * w_ba2 ^(1-ETATIL) * (pai_ba2^CHITIL / pai_ba1) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

wtil_ba1p = ((w_ba1p^(1-ETATIL) - ALFATIL * w_ba2p ^(1-ETATIL) * (pai_ba2p^CHITIL / pai_ba1p) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

h_ba1 = stil_ba1p * hd_ba1;
h_ba1p = stil_cup * hd_ba1p;


f1_ba1 = f2_ba1;
f1_ba1p = f2_ba1p;

profits_ba1 = output_ba1 - rk_ba1  * u_ba1 * k_ba1 - w_ba1 * hd_ba1 - NU * (1-1/r_ba1) * w_ba1 * hd_ba1;

m_ba1 = mh_ba1 + NU * w_ba1 * hd_ba1;

tauprofits_ba1 = tauk_ba1;%TAUPROFITS;%

tau_ba1 = tauk_ba1 * (u_ba1 *rk_ba1 - a_ba1 - qq_ba1 * DELTA) * k_ba1 + tauh_ba1 * hd_ba1 * w_ba1 + tauprofits_ba1 * profits_ba1;

%Utility function
uf_ba1 = (((c_ba1 - B * c_ba2)^(1-PHI4) * (1-h_ba1)^PHI4)^(1-PHI3)-1) / (1-PHI3);

objective_ba1 = uf_ba1;

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

e2_ba1 = -la_ba1 * (1+2*PHI1*v_ba1-2*(PHI1*PHI2)^(1/2)) + (1-PHI4) * (c_ba1 - B * c_ba2)^((1-PHI4)*(1-PHI3)-1) * (1-h_ba1)^(PHI4*(1-PHI3)) - B * BETTA * (1-PHI4) * (c_ba1p - B * c_ba1)^((1-PHI4)*(1-PHI3)-1) * (1-h_ba1p)^(PHI4*(1-PHI3));

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

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

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

e6_ba1 = -f2_ba1  + PHI4 * (c_ba1-B*c_ba2)^((1-PHI3)*(1-PHI4)) * (1-h_ba1)^(PHI4*(1-PHI3)-1) *  (w_ba1/wtil_ba1)^ETATIL * hd_ba1 + ALFATIL * BETTA * (pai_ba1p / pai_ba1^CHITIL)^ETATIL * f2_ba1p* (wtil_ba1p/wtil_ba1)^ETATIL;

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

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

e9_ba1 = -x2_ba1 + output_ba1 * ptil_ba1^(-ETA) + ALFA * BETTA * la_ba1p / la_ba1 * (pai_ba1^CHI / pai_ba1p )^(1-ETA) * (ptil_ba1 / ptil_ba1p)^(-ETA) * 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;

e12_ba1 = -stil_ba1p + (1-ALFATIL) * (wtil_ba1/w_ba1)^(-ETATIL) + ALFATIL * (w_ba2 / w_ba1) ^(-ETATIL) * (pai_ba1 / pai_ba2^CHITIL)^ETATIL * stil_ba1;

e13_ba1 = -assets_ba1p/r_ba1 - (1-1/r_ba1) * m_ba1 - tau_ba1 + assets_ba1 / pai_ba1 + g_ba1 + n_ba1;

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

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

 %Create vector of current constraints
constraints_ba1 = [e1_ba1; e2_ba1; e3_ba1; e4_ba1; e5_ba1; -e6_ba1;  e7_ba1; -e8_ba1;  e9_ba1; -e10_ba1; -e11_ba1; -e12_ba1; -e13_ba1; e14_ba1; e15_ba1];

%CONSTRAINTS_fu1RRENT

%Production Function
pf_fu1 = z_fu1 * (u_fu1*k_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 -1)^2);

d1qiv_fu1 = qiv_fu1/iv_fu1 -KAPA*(iv_fu1/iv_cu-1)*iv_fu1/iv_cu;
d2qiv_fu1p = (iv_fu1p/iv_fu1)^2*KAPA*(iv_fu1p / iv_fu1 -1);

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 / (z_fu1  * THETA * (u_fu1*k_fu1 /hd_fu1 )^(THETA-1));

ptil_fu1  = ((1 - ALFA * pai_fu1 ^(ETA-1) * pai_cu ^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

ptil_fu1p = ((1 - ALFA * pai_fu1p^(ETA-1) * pai_cup^(CHI*(1-ETA))) / (1- ALFA))^(1/(1-ETA));

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

wtil_fu1  = ((w_fu1 ^(1-ETATIL) - ALFATIL * w_cu ^(1-ETATIL) * (pai_cu^CHITIL / pai_fu1) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

wtil_fu1p = ((w_fu1p^(1-ETATIL) - ALFATIL * w_cup ^(1-ETATIL) * (pai_cup^CHITIL / pai_fu1p) ^(1-ETATIL) ) / (1-ALFATIL) )^(1/(1-ETATIL));

h_fu1 = stil_fu1p * hd_fu1;

f1_fu1 = f2_fu1;
f1_fu1p = f2_fu1p;

profits_fu1 = output_fu1 - rk_fu1  * u_fu1 * k_fu1 - w_fu1 * hd_fu1 - NU * (1-1/r_fu1) * w_fu1 * hd_fu1;

m_fu1 = mh_fu1 + NU * w_fu1 * hd_fu1;

tauprofits_fu1 = tauk_fu1;%TAUPROFITS;%

tau_fu1 = tauk_fu1 * (u_fu1 *rk_fu1 - a_fu1 - qq_fu1 * DELTA) * k_fu1 + tauh_fu1 * hd_fu1 * w_fu1 + tauprofits_fu1 * profits_fu1;

%Utility function
uf_fu1 = (((c_fu1 - B * c_cu)^(1-PHI4) * (1-h_fu1)^PHI4)^(1-PHI3)-1) / (1-PHI3);

objective_fu1 = uf_fu1;

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

e2_fu1 = -la_fu1 * (1+2*PHI1*v_fu1-2*(PHI1*PHI2)^(1/2)) + (1-PHI4) * (c_fu1 - B * c_cu)^((1-PHI4)*(1-PHI3)-1) * (1-h_fu1)^(PHI4*(1-PHI3)) - B * BETTA * (1-PHI4) * (c_fu1p - B * c_fu1)^((1-PHI4)*(1-PHI3)-1) * (1-h_fu1p)^(PHI4*(1-PHI3));

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

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

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

e6_fu1 = -f2_fu1  + PHI4 * (c_fu1-B*c_cu)^((1-PHI3)*(1-PHI4)) * (1-h_fu1)^(PHI4*(1-PHI3)-1) *  (w_fu1/wtil_fu1)^ETATIL * hd_fu1 + ALFATIL * BETTA * (pai_fu1p / pai_fu1^CHITIL)^ETATIL * f2_fu1p* (wtil_fu1p/wtil_fu1)^ETATIL;

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

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

e9_fu1 = -x2_fu1 + output_fu1 * ptil_fu1^(-ETA) + ALFA * BETTA * la_fu1p / la_fu1 * (pai_fu1^CHI / pai_fu1p )^(1-ETA) * (ptil_fu1 / ptil_fu1p)^(-ETA) * 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;

e12_fu1 = -stil_fu1p + (1-ALFATIL) * (wtil_fu1/w_fu1)^(-ETATIL) + ALFATIL * (w_cu / w_fu1) ^(-ETATIL) * (pai_fu1 / pai_cu^CHITIL)^ETATIL * stil_fu1;

e13_fu1 = -assets_fu1p/r_fu1 - (1-1/r_fu1) * m_fu1 - tau_fu1 + assets_fu1 / pai_fu1 + g_fu1 + n_fu1;

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

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

 %Create vector of current constraints
constraints_fu1 = [e1_fu1; e2_fu1; e3_fu1; e4_fu1; e5_fu1; -e6_fu1;  e7_fu1; -e8_fu1;  e9_fu1; -e10_fu1; -e11_fu1; -e12_fu1; -e13_fu1; e14_fu1; e15_fu1];