clc;clear all; close all;
%利用myfun函数求解yss
fun1 = @myfunc48;
options=optimset('Maxfuneval',50000, 'MaxIter', 3000, 'TolFun',1e-50);
[yss] = fsolve(fun1, [1, 1],options )
YLss = yss(1);
YHss = yss(2);


%%家庭部门参数
sigma=2;   % 朱军 
phic=0.7;   % 0.7朱军
phi=1.5; %王遥 劳动的跨期替代弹性
rho1=1; %王立勇(2018)% Carattini(2021)%
rho2=1;
omegaR=0.3;
beta=0.99;  %贴现率
delta = 0.025; 
chi=10;%%央行%%王博
%chih=2.395; %understanding%0.5882 王立勇（2018）

%%厂商有关的参数
%
epsilon=2;%中间品替代弹性--王博和徐飘扬
%epsilon=1.1;%中间品替代弹性--王博和徐飘扬
omega=0.58;%清洁中间品占比- 王博0.58 %央行0.6980 %
%omega=0.698;%清洁中间品占比- 王博0.58 %央行0.6980 %


alpha1 =0.33; %四个弹性参考王博
alpha2 =0.67;
alpha3=0.06;%刘
beta1=0.37;
beta2 =0.63;
epsilon1= 0.1810;%减排成本参数 %0.0326王博
epsilon2=2.8223;%减排成本参数 %2.6王博

%%环境有关的参数
%d0 = -0.0076;
%d1 = 6.176530605680263*10^(-05);
%d2 = 6.105320353459878*10^(-07);  %央行论文，王博论文
%d0=-0.0076;  参考 Gibson 和 Heutel（2020）；蓝天（2022）
%d1=8.10e-6 
%d2=1.05e-8;
d0=1.40*10^(-3);% Nordhus
d1=-6.67*10^(-6);
d2=1.46*10^(-8);
rho=0.15; % 产出的碳排放系数 武晓利（2017）
eta=0.0083; %环境自净能力 央行碳风险；朱军（2015）0.1
%Xss=2030; %王博


%%政府有关的参数
deltaG = 0.0125;
tau_crss=0.1506; 
tau_cnrss=0.1506; 
tau_kss= 0.3478;
tau_lss=0.1282;
tau_ehss=0.25;
tau_elss=0.25;
tau_xss=0.00758; %最优碳税王博和徐飘扬 0.00758%%央行0.0399

%%银行有关的参数
s=0.15;    %央行要求的存款准备金率 王遥0.15
%s=0.1125;    %22年4月份调整过后是11.25%
ksi=0.025;  %银行的营运成本占资金规模的比重，成蕾



Pss = 1;  %仅假设最终产品价格为1  
ALss = 1;
AHss=1;
BHss = 1; 
phiBss=0.159; %师姐
phiTRANSss =0.08; %师姐
phiIGss=0.123; %结构性减税
QUBAR=1000;
XBAR=5;



%%%%%%%%有关银行的变量稳态
RRss=1.0175; %碳减排支持工具的再贷款利率为1.75;;%蓝天0.0175 %王遥1.03
Mss=0.6; %碳减排支持工具的再贷款比例为60%
RDss = 1/beta;   %存款利率 家庭的一阶条件得出
thetass=((RDss-1)+ksi)/(1-s);
RHss=1+thetass;

theta2ss=1+thetass-RRss;
RLss=1+thetass-theta2ss*Mss;

%RLss=1+thetass+(RRss-1)*Mss-thetass*Mss;
Fss=0.2*RLss; %贴息率 王遥
RFss=RLss-Fss;





RKss = Pss*((1+tau_crss)/(1-tau_kss))*((1/beta)-(1-delta));  %稳态利率由家庭决定
Yss=((omega^(1/epsilon))*(YLss^((epsilon-1)/epsilon))+((1-omega)^(1/epsilon))*(YHss^((epsilon-1)/epsilon)))^(epsilon/(epsilon-1));
PLss=((omega*(Pss^epsilon)*Yss)/YLss)^( 1/epsilon);
PHss=(((1-omega)*(Pss^epsilon)*Yss)/ YHss) ^(1/epsilon);  
Zss=((tau_xss*rho)/(epsilon1*epsilon2*BHss)) ^(1/(epsilon2-1));
CEHss=((1-Zss)*YHss*rho)/BHss;      %棕色企业的碳排放量 改为线性了
TCREHss=epsilon1*(Zss^epsilon2)*YHss;
TCCEHss=tau_xss*CEHss;
INHss=CEHss/YHss;


IGss= phiIGss *Yss;
TRANSss=phiTRANSss*Yss;
KGss=IGss/deltaG;
Xss=(XBAR+CEHss)/eta;
QUss= QUBAR-Xss;
DXss=d0+d1*Xss+d2*(Xss^2);  %定义DX=d[Xt]

KLss=(alpha1*PLss*(1-tau_elss)*YLss)/(RKss*RFss);
KHss=((beta1*YHss)/(RHss*RKss))*((1-tau_ehss)*PHss-(tau_xss*(1-Zss)*rho)/BHss-epsilon1*(Zss^epsilon2));
KPss=KHss+KLss;
IPss=delta*KPss;
IHss=delta*KHss;
ILss=delta*KLss;
LHss=(YHss/((1-DXss)*AHss*(KHss^beta1)))^(1/beta2);
LLss=(YLss/((1-DXss)*ALss*(KLss^alpha1)*(KGss^alpha3)))^(1/alpha2);
WLss=(alpha2*YLss*(1-tau_elss)*PLss)/(LLss*RFss);
WHss=((beta2*YHss)/(LHss*RHss))*((1-tau_ehss)*PHss-tau_xss*(1-Zss)*rho-epsilon1*(Zss^epsilon2));    
LRHss=LHss;
LRLss=LLss;
LNRHss=LHss;
LNRLss=LLss;
LRss=(LRHss^(1+rho1)+ LRLss^(1+rho1))^(1/(1+rho1));
LNRss=(LNRHss^(1+rho2)+ LNRLss^(1+rho2))^(1/(1+rho2));
LAMRss=((LRss^(phi-rho1))*(LRHss^rho1))/(WHss*(1-tau_lss));
LAMNRss=((LNRss^(phi-rho2))*(LNRHss^rho2))/(WHss*(1-tau_lss));


HHss=WHss*LHss+RKss*KHss;
HLss=WLss*LLss+RKss*KLss;
Hss=HLss+HHss;
HRss=Mss*HLss;
Dss=(Hss-HRss)/(1-s);
Css= (Yss-TCREHss+TRANSss+ Fss*HLss- tau_lss*WHss*LHss-tau_lss*WLss*LLss-tau_kss*RKss*KPss-tau_elss*PLss*YLss-tau_ehss*PHss*YHss-tau_xss*CEHss-(1+tau_crss)*IPss)/(1+tau_crss);
Tss=tau_crss*Pss*(Css+IPss)+tau_lss*WHss*LHss+tau_lss*WLss*LLss+tau_kss*RKss*KPss+tau_ehss*PHss*YHss+tau_elss*PLss*YLss +tau_xss*CEHss;
CNRss=Css;
CRss=Css;
Gss=Tss-IGss-TRANSss- Fss*HLss;
Qss=LAMRss*Pss*(1+tau_crss);

INRss= WHss*LRHss*(1-tau_lss)+ WLss*LRLss*(1-tau_lss)+RKss*KPss*(1-tau_kss)+RDss*Dss;
INNRss= WHss*LNRHss*(1-tau_lss)+ WLss*LNRLss*(1-tau_lss)+Pss*TRANSss; 
INDss=INRss/INNRss;
inrss=INRss/Pss;
innrss=INNRss/Pss;
indss=inrss/innrss;


shareKss=RKss*KPss/(RKss*KPss+WLss*LLss+ WHss*LHss);
shareLss=(WLss*LLss+ WHss*LHss)/(RKss*KPss+WLss*LLss+ WHss*LHss);
PIss=1;


F1= LAMRss-((LRss^(phi-rho1))*(LRLss^rho1))/(WLss*(1-tau_lss));
F2=CRss-((LAMRss*Pss*(1+tau_crss))/(phic*(QUss^((1-sigma)*(1-phic)))))^(1/(phic*(1-sigma)-1));
F=[F1;F2];
Yss
B1=(1+tau_cnrss)*CNRss-(1-tau_lss)*WHss*LNRHss-(1-tau_lss)*WLss*LNRLss-Pss*TRANSss;
B2=(1+tau_crss)*(CRss+IPss)+Dss-(1-tau_lss)*WHss*LRHss-(1-tau_lss)*WLss*LRLss-(1-tau_kss)*RKss*KPss-RDss*Dss;


save _paramc48 theta2ss TCCEHss INHss  PIss  INDss inrss innrss indss XBAR Dss RDss sigma phic phi rho1 rho2 beta delta deltaG chi epsilon omega omegaR  alpha1 alpha2 alpha3 beta1 beta2 epsilon2 epsilon1 rho d0 d1 d2 eta QUBAR tau_crss tau_cnrss tau_lss tau_kss tau_ehss tau_elss tau_xss   phiIGss  phiTRANSss PHss PLss Pss ALss AHss    RKss Yss YHss KHss LHss WHss WLss BHss Zss CEHss TCREHss YLss KLss IGss KGss LLss KPss IPss   Xss QUss DXss Css Tss Gss LAMRss LAMNRss Qss IHss ILss CNRss CRss LNRss LRss LRHss LRLss  LNRHss LNRLss  INRss INNRss shareKss shareLss TRANSss RRss RHss RLss RFss HHss HLss Hss HRss thetass Fss s ksi Mss
