clear
syms y yW c i G gdp k lambda q H m m2 HS HP yS yP pW pS pP kS kP nS nP bS bP w wSe wPe RSL RPL Rk Rd RS RP RSI RPI pi tao taoT miu  v vS vP U omgS omgP fomgS fomgP gomgS gomgP fdomgS fdomgP gdomgS gdomgP FomgS FomgP;% 算上技术冲击AS共59 variables ，10个参数方程，1个冲击方程，48个内生方程
syms ys yWs cs is Gs gdps ks lambdas qs Hs ms m2s HPs HSs ySs yPs pWs pSs pPs kSs kPs nSs nPs bSs bPs ws wSes wPes RSLs RPLs Rks Rds RSs RPs RSIs RPIs pis taos taoTs mius vs vSs vPs Us omgSs omgPs fomgSs fomgPs gomgSs gomgPs fdomgSs fdomgPs gdomgS gdomgPs FomgSs FomgPs
syms AP AS APm ASm Abar rhoa beta eta delta g faiK psiH psiM  alphaP alphaS  theta ksiS ksiP ebslon pis taoTs vs Us taos mius faiTgamma faiTpi faitrv faitru faitaogamma faitaopi faitaov faitaou faimiugamma faimiupi faimiuv faimiuu gc mS mP  fai sigma omgm kappa lS lP;                                    
 %在论文中，RSL写成了RSB
% Parameter Calibration
beta=0.995; eta=1.5; psiH=16; delta=0.05; APms=1; ASms=1; 
ebslon = 10; g = 1.0125; alphaS = 0.5; alphaP=0.85; theta = 0.97; sigma = 3; 
mS = 0.15; mP = 0.15; ksiS = 0.97; ksiP = 0.92;%mP设定的不能太小，否则会出现gdomgS/Ps为负，无解的情形，而mS设定也要适度，否则要么HSs太大，要么为0;ksiP不能太大，也会出现无解的情形。
gc=0.143; fai=0.4; lS = 0.379;  lP = 0.001;
psiM =0.06; omgm = 0.37; kappa =1.587 ; 
faiP =4.6; rhoaS = 0.95; rhoaP = 0.95; faiK = 0.142; 

x = sym('x', [1,5])%1行两列
HSs = x(1)
omgSs = x(2)
omgPs = x(3)
cs = x(4)
Hs = x(5)

%需要对fomgs等进行重新设定调整
fomgSs = 1/(kappa-1) * omgm^kappa * omgSs^(1-kappa)
fomgPs = 1/(kappa-1) * omgm^kappa * omgPs^(1-kappa)

fdomgSs = - omgm^kappa * omgSs^(-kappa)
fdomgPs = - omgm^kappa * omgPs^(-kappa)

gomgSs = (1 - lS) * (1 - mS) * kappa * omgm / (kappa-1) + lS * omgSs + (1-lS) * (1- (1-mS) * kappa / (kappa-1)) * omgSs^ (1-kappa) * omgm^kappa
gomgPs = (1 - lP) * (1 - mP) * kappa * omgm / (kappa-1) + lP * omgPs + (1-lP) * (1- (1-mP) * kappa / (kappa-1)) * omgPs^ (1-kappa) * omgm^kappa

gdomgSs = lS + (1-lS) * ( 1-mS * kappa ) * omgSs^(-kappa) * omgm^kappa
gdomgPs = lP + (1-lP) * ( 1-mP * kappa ) * omgPs^(-kappa) * omgm^kappa
 
FomgSs = kappa * omgSs^(1-kappa) * omgm^kappa / (1-kappa) - kappa * omgm / (1-kappa)
FomgPs = kappa * omgPs^(1-kappa) * omgm^kappa / (1-kappa) - kappa * omgm / (1-kappa)

%求稳态,求解对于初值非常敏感，设置有些初始值，会得到实数解，有些则得到复数解

ASs = 1; APs = 1; taos = 0.15; taoTs = 0.0075; pis = 1+0.035/4;
%Us = 0.0406;  vs = 0.04; %三个方程无效，则校准三个变量

mius = 1;
qs = 1; 
Rks = qs * g / beta - qs * ( 1 - delta )
Rds = pis * g / beta
RSs = (Rds - 1) / (1-taos) + 1
RPs = (Rds - 1) / (1-taos + taoTs) + 1
HPs = Hs -HSs
ws = cs * psiH * (Hs) ^eta
kSs = (1-alphaS) * ws * HSs / (alphaS * theta)/ Rks
kPs = (1-alphaP) * ws * HPs / (alphaP * theta)/ Rks

ks = g * (kSs + kPs)

is = (1 - (1-delta) / g) * ks
%A8
ySs = ASs * kSs ^(1-alphaS) * HSs^(theta * alphaS)
% A16
yPs = APs * kPs ^(1-alphaP) * HPs^(theta * alphaP)

%A26

pWs =(ebslon - 1) / ebslon
% A35
yWs = (fai * ySs^((sigma - 1)/sigma) + (1-fai) * yPs^((sigma - 1)/sigma))^(sigma / (sigma - 1))

%A36
pSs = pWs * fai * yWs^(1/sigma) / ySs^(1 / sigma)
%A37
pPs = pWs * (1 - fai) * yWs^(1/sigma) / yPs^(1 / sigma)
%A9
% A17
%计算稳态时用WH，验证稳态时换bn
%A10
wSes = alphaS * (1-theta) * ws * HSs / (alphaS * theta)
%A18
wPes = alphaP * (1-theta) * ws * HPs / (alphaP * theta)

%A12
RSIs = pSs * ySs / (ws * HSs / (alphaS * theta))
%A20
RPIs = pPs * yPs / (ws * HPs / (alphaP * theta))

%A15
nSs = wSes + ksiS * RSIs * ws * HSs / (alphaS * theta) * fomgSs
% A23
nPs = wPes + ksiP * RPIs * ws * HPs / (alphaP * theta) * fomgPs

%A13
bSs = gomgSs * ws * HSs / (alphaS * theta) * (RSIs / RSs)
%A21
bPs = gomgPs * ws * HPs / (alphaP * theta) * (RPIs / RPs)

% //A14///////////////////////////////////////////////////////
fun1 = (bSs + nSs/(pis * g)) - ws *  HSs / (alphaS * theta) 
fun11 = nSs / (pis * g) / (ws * HSs / (alphaS * theta)) - (- gdomgSs * RSIs * fomgSs / (fdomgSs * RSs))
% //A22
fun2 = (bPs + nPs/(pis * g)) - ws *  HPs / (alphaP * theta) 
fun22 = nPs / (pis * g) / (ws * HPs / (alphaP * theta)) - (- gdomgPs * RPIs * fomgPs / (fdomgPs * RPs))

%A45
RSLs = omgSs * RSIs * ws * HSs / (alphaS * theta) / bSs
%A46
RPLs = omgPs * RPIs * (nPs/(pis * g) + bPs) / bPs
%A43
vSs = 1 - RSIs * gomgSs * ws * HSs / (alphaS * theta) / (bSs * RSLs)
%A44
vPs = 1 - RPIs * gomgPs * ws * HPs / (alphaP * theta) / (bPs * RPLs)

% //A41//////////////////////////////////////////////////////
 vs = bSs * vSs / (bPs + bSs) + bPs * vPs / (bPs + bSs)
%A40
gdps = (is + cs) / (1 - gc)
%A32
Gs = gdps * gc;
% A33
ys = cs + is + Gs + RSIs * mS * FomgSs * ws * HSs / (alphaS * theta) + RPIs * mP * FomgPs * ws * HPs / (alphaP * theta);
% //A34////////////////////////////////////////////////////
fun4 =ys - yWs 
lambdas = 1 / cs
% A5
ms = cs * psiM / (1 - beta / (pis * g))
%A31
m2s = ms + bSs / (1 - taos) + bPs / (1 - taos + taoTs)
Us = 1- Hs
format long
fun = [fun1,fun11, fun2, fun22, fun4]
eq = matlabFunction (fun, 'Vars', {[x(1), x(2), x(3), x(4), x(5)]})
options=optimset('MaxIter',1000000000000000, 'Maxfuneval', 10000000000000000);
xs = fsolve(eq, [0.12,0.1, 0.5, 0.5, 1], options)

HSs = xs(1)
omgSs = xs(2)
omgPs = xs(3)
cs = xs(4)
Hs = xs(5)
save SteadyState HSs omgSs omgPs cs Hs

