%% Two Sector Open Economy DSGE model 
%% 2022-08-03

beta = 0.99;
pai = 1;
Rn = 1/beta;
Rn_star = Rn;
sigma =2;
a = 0.5;
ai = 0.5;
eta = 1.5;
etai = 1.5;
v = 1.3;
delta = 0.025;


paiH = 1;
etag = 0.15;
eta_ex = 0.03;

epsilon = 6;
pw = (epsilon-1)/epsilon;

alpha = 0.5;
alphap = 0.3;

s = 1;
pI = 1;
Qk = 1;
Qk_s = Qk;
Qk_p = Qk;


omega = 0.5;
phis = 0.97;
phip = 0.69;

Ap = 1.42;
As = 1;
d = 0.0025;

mu = 3;
psYs_ppYp = 45/55;
L = 1/3;
Ls_Lp = psYs_ppYp;
Lp = L/(1+Ls_Lp);
Ls = L-Lp;
Ks_Kp = psYs_ppYp;
Ys_Yp = (As/Ap)*(Ks_Kp)^alpha*(Ls_Lp)^(1-alpha);
ps_pp = psYs_ppYp/Ys_Yp;
psi = (Ys_Yp^(1/mu)*ps_pp)/(1+Ys_Yp^(1/mu)*ps_pp);
pp = pw/(psi^mu*ps_pp^(1-mu)+(1-psi)^mu)^(1/(1-mu));
ps = ps_pp*pp;


%% solve home price using fsolve
%pH = priceH(1, a, eta, s, ai, etai, pI, phi, omega, delta, alpha, Lp, Ls, Ap, As, psi, mu, etag, eta_ex, pp, ps, pw, psYs_ppYp)
pH = 1;
CH_C = (1-a)*pH^(-eta);
CF_C = a*s^(-eta);
IH_I = (1-ai)*pH^(-etai);
IF_I = ai*(s/pI)^(-etai);
%Rk = pH/(phis*((1-omega)*pH+omega*s));
%Rk = 1/0.99+100/10000;
%phi = pH/(Rk*((1-omega)*pH+omega*s));
Rk_s = pH/(phis*((1-omega)*pH+omega*s));
Rk_p = pH/(phip*((1-omega)*pH+omega*s));
Rr_s = Rk_s;
Rr_p = Rk_p;
Zs = Rk_s-(1-delta);
Zp = Rk_p-(1-delta);
Yp_Kp = (Zp/alphap)/pp;
Ys_Ks = (Zs/alpha)/ps;
Kp = (Yp_Kp)^(1/(alphap-1))*(1/Ap)^(1/(alphap-1))*Lp;
Ks = (Ys_Ks)^(1/(alpha-1))*(1/As)^(1/(alpha-1))*Ls;
Yp = Yp_Kp*Kp;
Ys = Ys_Ks*Ks;
Is = delta*Ks;
Ip = delta*Kp;
I = Is+Ip;
IH = IH_I*I;
IF = IF_I*I;
Y = Ys/(psi^mu*(ps/pw)^(-mu));
g = Y*etag;
ex = Y*eta_ex;
CH = Y-IH-g-ex;
C = CH/CH_C;
CF = CF_C*C;
w = (1-alphap)*pp*Yp/Lp;
xi = (w/C^sigma)/L^v;


up = 0.1;
Np = Kp/(((Rr_p-d)/Rn)^(1/up));
bp = Kp-Np;
us = 0.1;
Ns = Ks/((Rr_s/Rn)^(1/us));
bs = Ks-Ns;
b = bs+bp;

PSI = 1;
Omega_p = 75;
paiH_ss = paiH;
Rn_ss = Rn;
Y_ss = Y;
pai_ss = pai;
s_ss = s;
ex_ss = ex;
Rn_star_ss = Rn_star;
d_ss = d;

Yp_s = 1/Ys_Yp;
bp_b = bp/b;
Kp_s = 1/Ks_Kp;
Lp_s = 1/Ls_Lp;
Yp_Y = 1/(Ys_Yp+1);


SS = [Rn pai C s L w CH CF pH...
Rr_s Qk_s Ks Ns Ls bs Is ps Ys Zs Rk_s...
Rr_p Qk_p Kp Np Lp bp Ip pp Yp Zp Rk_p...
IH I IF pI Y pw paiH b...
Rn_star g ex Yp_s bp_b Kp_s Lp_s CH_C Yp_Y];        
 


param = [beta sigma xi v eta a...
delta PSI us alpha phis omega...
up...
etai ai psi mu...
epsilon Omega_p paiH_ss Rn_ss Y_ss pai_ss s_ss...
etag ex_ss Rn_star_ss d_ss phip alphap];














