

%Labor is calibrated
% keyboard
L1_u=Lss1;
L1_c = Lss1;
L2_u=Lss2;
L2_c=Lss2;

%CBDC parameters
midc1=Theta1*mi2;
midc2=Theta1*mi2;
sigmadc1 = (sigmam1 + sigmam1*(1-Theta1));
sigmadc2 = (sigmam2 + sigmam2*(1-Theta2));
rdc1_u=1;
rdc1_c = 1;
%MC derives from optimality (i.e. in the equilibrium the mark-up is
%proportional to mrk power)
MC1=((ni1-1)/(ni1))*(1+tau1);
MC2=((ni2-1)/(ni2))*(1+tau2);

%Inflation is 1
pi1=1;
pi2=1;


%prices on domestic and foreign goods equal CPI
p11=1;
p22=1;
p21=1;
p12=1;
PO11=1;
PO22=1;
PO21=1;
PO12=1;
pi11=1;
pi22=1;
pi21=1;
pi12=1;

%Exchange rate
NEER12=1;
REER12=1;

%Risk free rate
% beta1=1;
r1=1/beta1;
r2=1/beta2;

%Rate of return on risky assets (this is from eq 17 and 18) net of loss
rk1=((1-beta1*(1-deltak1))/(beta1));
rk2=((1-beta2*(1-deltak2))/(beta2));

%Capital (from production function)
K1=(rk1/(alpha1*MC1))^(1/(alpha1-1));
K2=(rk2/(alpha2*MC2))^(1/(alpha2-1));
%Investments
I1_u=deltak1*K1;
I1_c=deltak1*K1;
I2_u=deltak2*K2;
I2_c=deltak2*K2;
ratio1=1;
ratio2=1;
%Investment consumtion of domestic goods
I1_1 = omega11*(I1_c+I1_u);
I2_2 = omega22*(I2_c+I2_u);

%Investment consumption of foreign goods
I1_2 = (omega12)*(I1_c+I1_u);
I2_1 = (omega21)*(I2_c+I2_u);

%Total output
Ytot1=K1^(alpha1);
Ytot2=K2^(alpha2);

%real wages
w1_u=(1-alpha1)*MC1*K1^(alpha1);
w1_c=(1-alpha1)*MC1*K1^(alpha1);
w2_u=(1-alpha2)*MC2*K2^(alpha2);
w2_c=(1-alpha2)*MC2*K2^(alpha2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Define domestic output, in the equilibrium
Y1=omega11*Ytot1;
Y2=omega22*Ytot2;

%Define foreign output in the equilibrium
% //foreign demand (country 1)
% C2_1 + G2_1 + I2_1 = Y21;
% C3_1 + G3_1 + I3_1 = Y31;
Y21 = omega21*Ytot2;
% //foreign demand (country 2)
% C1_2 + G1_2 + I1_2 = Y12;
% C3_2 + G3_2 + I3_2 = Y32;
Y12 = omega12*Ytot1;
% //foreign demand (country 3)
% C1_3 + G1_3 + I1_3 = Y13;
% C2_3 + G2_3 + I2_3 = Y23;


%Degine G from the share of total output (share is calibrated 0.2)
G1=GY1ss*Ytot1;
G2=GY2ss*Ytot2;
G1_1=omega11*G1;
G2_2=omega22*G2;
G1_2 = (omega12)*G1;
G2_1 = (omega21)*G2;

%define consumption. First notice that
%domestic consumption of domestic goods is
%  Yh = Ch_h +Gh_h +Ih_h;
%domestic consumption for foreign goods is
%  Ch_f + Gh_f + Ih_f = Yf_for;
C1_1  = Y1 - (G1_1 + I1_1);
C2_2  = Y2 - (G2_2 + I2_2);

% foreign demand (country 1)
C2_1  = Y21 - (G2_1 + I2_1);

% foreign demand (country 2)
C1_2  = Y12 - (G1_2 + I1_2);


% demand for domestic goods in consumption
C1_u = C1_1 / omega11;
C1_c  = C1_1 / omega11;
C2_u = C2_2 / omega22;
C2_c = C2_2 / omega22;
C1 = C1_u+C1_c;
C2 = C2_u+C2_c;

% keyboard
%derive optimal shares

if phib1==0
    B12=0;
    B13=0;
else
B12  = (beta1 * r2-1)/phib1 ;
end

if phib2==0
    B21=0;
    B23=0;
else
B21  = (beta2* r1-1)/phib2 ;
end


B11=-(B21);
B22=-(B12);


% keyboard
%Definition of Lambda
% Lambda1=((C1 - h1*C1)^(-sigma1)) - beta1*h1*((C1 - h1*C1)^(-sigma1));
% Lambda2=((C2 - h2*C2)^(-sigma2)) - beta2*h2*((C2 - h2*C2)^(-sigma2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CASH, CBDC AND LAGRANGIAN MULTIPLIERS BLOCK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%*****************************HOME COUNTRY*********************************
F = @(x) [((C1_u)^(rho_dc))*midc1*x(5)*x(2)^(-rho_dc) - chi_dc1*((2*x(1)*(Gamma1 - x(3)))/(x(1)+x(2))^2)-x(4)+beta1*x(4)*rdc1_u/pi1;
          ((C1_u)^(rho_dc))*mi1*x(5)*x(1)^(-rho_dc) - chi_dc1*(2*(Gamma1 - x(3))*(x(3)/(x(1)+x(2)) - 1/(x(1)+x(2)))) - x(4) + beta1*x(4)*storage1/pi1(+1);
          x(3) - x(1)/(x(1)+x(2));
          ((C1_u - h1*C1_u)^(-sigma1)) - beta1*h1*((C1_u - h1*C1_u)^(-sigma1))-x(4)-x(5);
          C1_u-(mi1*x(1)^(1-rho_dc)+midc1*x(2)^(1-rho_dc))^(1/(1-rho_dc))];  
%initial value for x(2) guessed at 0.5
options = optimoptions('fsolve','Display','none','FunctionTolerance',1e-8,'OptimalityTolerance',1e-8); 
[x,fval,exitflag,output,jacobian] = fsolve(F,[1 1 0.5 1 1],options); 
trs1_u=x(3);
M1=x(1);
DC1_u=x(2);
Lambda1_u=x(4);
gam1=x(5);

F = @(x) [((C1_c)^(rho_dc))*midc1*x(5)*x(2)^(-rho_dc) - chi_dc1*((2*x(1)*(Gamma1 - x(3)))/(x(1)+x(2))^2)-x(4)+beta1*x(4)*rdc1_c/pi1;
          ((C1_c)^(rho_dc))*mu_d*x(5)*x(1)^(-rho_dc) - chi_dc1*(2*(Gamma1 - x(3))*(x(3)/(x(1)+x(2)) - 1/(x(1)+x(2)))) - x(4) + beta1*x(4)*storage1/pi1(+1);
          x(3) - x(1)/(x(1)+x(2));
          ((C1_c - h1*C1_c)^(-sigma1)) - beta1*h1*((C1_c - h1*C1_c)^(-sigma1))-x(4)-x(5);
          C1_c-(mu_d*x(1)^(1-rho_dc)+midc1*x(2)^(1-rho_dc))^(1/(1-rho_dc))];  
%initial value for x(2) guessed at 0.5
options = optimoptions('fsolve','Display','none','FunctionTolerance',1e-8,'OptimalityTolerance',1e-8); 
[x,fval,exitflag,output,jacobian] = fsolve(F,[1 1 0.5 1 1],options); 
trs1_c=x(3);
I1_c=x(1);
DC1_c=x(2);
Lambda1_c=x(4);
gam1=x(5);
%**************************FOREIGN COUNTRY*********************************
F = @(x) [((C2_u)^(rho_dc))*midc2*x(5)*x(2)^(-rho_dc)-chi_dc2*((2*x(1)*(Gamma2-x(3)))/((x(1)+x(2))^2)) - x(4) - x(4)*phidc2*(x(2))+beta2*x(4)*rdc1_u*(1/pi2);
         ((C2_u)^(rho_dc))*mi2*x(5)*x(1)^(-rho_dc)-chi_dc2*(2*(x(3)/(x(1)+x(2)) - 1/(x(1)+x(2)))*(Gamma2 - x(3))) - x(4) + beta2*x(4)*storage2/pi2;
         x(3) - x(1)/(x(1)+x(2));
         ((C2_u - h2*C2_u)^(-sigma2))-beta2*h2*((C2_u - h2*C2_u)^(-sigma2))-x(4)-x(5);
         C2_u-(mi1*x(1)^(1-rho_dc)+midc1*x(2)^(1-rho_dc))^(1/(1-rho_dc))];
options = optimoptions('fsolve','Display','none','FunctionTolerance',1e-8,'OptimalityTolerance',1e-8); 
[x,fval,exitflag,output,jacobian] = fsolve(F,[1 1 0.5 1 1],options); 
trs2_u=x(3);
M2=x(1);
DC2_u=x(2);
Lambda2_u=x(4);
gam2=x(5);


F = @(x) [((C2_c)^(rho_dc))*midc2*x(5)*x(2)^(-rho_dc)-chi_dc2*((2*x(1)*(Gamma2-x(3)))/((x(1)+x(2))^2)) - x(4) - x(4)*phidc2*(x(2))+beta2*x(4)*rdc1_c*(1/pi2);
         ((C2_c)^(rho_dc))*mu_d*x(5)*x(1)^(-rho_dc)-chi_dc2*(2*(x(3)/(x(1)+x(2)) - 1/(x(1)+x(2)))*(Gamma2 - x(3))) - x(4) + beta2*x(4)*storage2/pi2;
         x(3) - x(1)/(x(1)+x(2));
         ((C2_c - h2*C2_c)^(-sigma2))-beta2*h2*((C2_c - h2*C2_c)^(-sigma2))-x(4)-x(5);
         C2_c-(mu_d*x(1)^(1-rho_dc)+midc1*x(2)^(1-rho_dc))^(1/(1-rho_dc))];
options = optimoptions('fsolve','Display','none','FunctionTolerance',1e-8,'OptimalityTolerance',1e-8); 
[x,fval,exitflag,output,jacobian] = fsolve(F,[1 1 0.5 1 1],options); 
trs2_c=x(3);
I2_c=x(1);
DC2_c=x(2);
Lambda2_c=x(4);
gam2=x(5);

%define Q
Q1 = beta1*(Lambda1_u+Lambda1_c)* rk1/(1-beta1*(1-deltak1));
Q2 = beta2*(Lambda2_u+Lambda2_c)* rk2/(1-beta2*(1-deltak2));

% Definition of ratio
ratio1=1;
ratio2=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Price setting, home good sold at home
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition of K & F for the domestic goods 
FF11 = Y1*(Lambda1_c+Lambda1_u)/(1- beta1 * csi1);
FF22 = Y2*(Lambda2_c+Lambda2_u)/(1- beta2 * csi2);
KK11 = FF11;
KK22 = FF22;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Price setting, exports of country 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FF21 = (Lambda1_c+Lambda1_u)*Y21/(1-beta1*csi1);
KK21 = (Lambda1_c+Lambda1_u)*Y21/(1-beta1*csi1);
PO21=1;
pi21 = 1;
p21 = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Price setting, exports of country 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PO12=1;                                       
FF12 = (Lambda2_c+Lambda2_u)*Y12/(1-beta2*csi2);
KK12 = (Lambda2_c+Lambda2_u)*Y12/(1-beta2*csi2);
pi12 = 1;
p12 = 1;


% now back out chi
chi1 = (w1_c+w1_u)*(Lambda1_c+Lambda1_u) / ((L1_u+L1_c)^(phi1));
chi2 = (w2_c+w2_u)*(Lambda2_c+Lambda2_u) / ((L2_u+L2_c)^(phi2));


% Aggregators of output are 1 in ss
d11 = 1;
d22 = 1;
d21  = 1;
d12  = 1;


% define the term in output needed for computationally solve output
% (sometimes there is a discrpancy of 10^-8 and dynare does not solve)
% keyboard
agg1 = d11*Y1+d21*Y21 - Ytot1;
agg2 = d22*Y2+d12*Y12 - Ytot2;


%compute the tax
%Notice that profits from invoicing in foreign goods are included in the
%budget constraint, assume they go to reduce taxes
% (B11+REER12*B12+REER13*B13+C1+I1+M1) = 
% (r2(-1)-1)*BG12(-1)*REER12/pi1+r1(-1)*B11(-1)/pi1 + REER12*r2(-1)*B12(-1)/pi1 + REER13*r3(-1)*B13(-1)/pi1 -(phib1/2)*((REER12*B12)^2)- (phib1/2)*((REER13*B13)^2) +rk1*K1 + ((chi1 * L1^(phi1))/Lambda1)*L1 +T1ss + M1(-1)/pi1 + (size1/size1)*(p11-MC1*d11)*Y1 + (size2/size1)*(((1-s1b)*1+s1b*REER12)*p21-MC1*d21)*Y21 + (size3/size1)*(((1-s1-s1b)*REER13+s1*1+s1b*REER12)*p31-MC1*d31)*Y31;
part1_u = (C1_u+REER12*B12+(B11+I1_u+M1+DC1_u)*zeta1);
part2_u =M1*storage1/pi1 + DC1_u*rdc1_u/pi1+r1*B11/pi1 + REER12*r2*B12/pi1 -(phib1/2)*((REER12*B12)^2) +rk1*K1 + (w1_c+w1_u)*(L1_u+L1_c)+ (p11-MC1*d11)*Y1 + (REER12*p21-MC1*d21)*Y21 ; 
T1ss_u = part1_u-part2_u;

part1_c = (C1_c+(I1_c+DC1_c)*zeta2);
part2_c = DC1_c*rdc1_c/pi1 +rk1*K1 + (w1_c+w1_u)*(L1_c+L1_c)+ (p11-MC1*d11)*Y1 + (REER12*p21-MC1*d21)*Y21 ; 
T1ss_c = part1_c - part2_c;

part1_u = ((B21/REER12)+C2_u+(B22+I2_u+M2+DC2_u)*zeta2);
part2_u = M2*storage2/pi2 + DC2_u*rdc1_u/pi2 + r2*B22/pi2 + (1/REER12)*r1*B21/pi2  -(phib2/2)*((B21/REER12)^2)+rk2*K2 + (w2_c+w2_u)*(L2_u+L2_c) + (p22-MC2*d22)*Y2 + (p12*(1/REER12)-MC2*d12)*Y12 ;
T2ss_u = part1_u-part2_u;

part1_c = (C2_u+(I2_c+DC2_c)*zeta2);
part2_c = DC2_c*rdc1_c/pi2+rk2*K2 + (w2_c+w2_u)*(L2_u+L2_c) + (p22-MC2*d22)*Y2 + (p12*(1/REER12)-MC2*d12)*Y12 ;
T2ss_c = part1_c-part2_c;

%utility function
U1_u = (log(C1_u - h1*C1_u)) - ((chi1/(1+phi1))*L1_u^(1+phi1)) - chi_dc1*(Gamma1-trs1_u)^2 ;
U1_c = (log(C1_c - h1*C1_c)) - ((chi1/(1+phi1))*L1_c^(1+phi1)) - chi_dc1*(Gamma1_c-trs1_c)^2 ;
U2_u = (log(C2_u - h2*C2_u)) - ((chi2/(1+phi2))*L2_u^(1+phi2)) - chi_dc2*(Gamma2-trs2_u)^2 ;
U2_c = (log(C2_c - h2*C2_c)) - ((chi2/(1+phi2))*L2_c^(1+phi2)) - chi_dc2*(Gamma2_c-trs2_c)^2 ;


%welfare
WEL1=(U1_u+U1_c)/(1-beta1);

WEL2=(U2_u+U2_c)/(1-beta2);


% Shocks are zero in the steady state
RR1=0;
RR2=0;


A1=0;
A2=0;


G1ss = G1;
G2ss = G2;


EC1 = 0;
EC2 = 0;


PSI1 = 0;
PSI2 = 0;


P1 = 0;
P2 = 0;

UIP = 100;
ly1=log(Ytot1);
ly2=log(Ytot2);