
clear all; 
clc; 
close all;
     
%% targets moments
piss=1;       % domestic inflation (IPF)-Excel
rss = (1+0.0376)*piss; % nominal interest rate 孔敏主编：《南开经济指数资料汇编》，中国社会科学出版社1988 年版1929-1936年的数据（1911-1936 excel "interest rate")
sp = 0.06;           % China lending rate from China_data 戴建兵 借款利率大萧期间大概平均是10%左右 用10%-0.0376=0.6左右

lss = 0.153;         % land 
phiss = 5;           % leverage 
hss = 1;             % hours
Abar = 0.61;         % 巫宝三的《中国国民所得(1933年)》第39页的内容，中国在那个时期确实是一个以农业为主的国家，其中约有61%的国民所得来自农业
bbar = 0.17;         % Goverment-gdp ratio
rrzss = 1+0.0508;      % foreign real interest rate （global financial database China_data "US UK_interest rates" 1929-1936

%% parameters
varphiH = 1;        % inverse of Frisch elasticity (IPF). ALternative: 3.79 (AQ19)
varphiL = 1.8;
alphak=0.33;        % elasticity of production wrt capital
alphah=0.33;        % elasticity of production wrt labor
delta=0.1;        % depreciation rate 
eps=10;             % elasticity of substitution btw differentiated goods 
alphaN = 0.2;       % trade openness 
alphaA	= 0.1;		% share of imported goods agricultural
%epsL=1.5;           % elasticity of substitution between means of payment (CPP)  
deltaL = 0.1;		% land decay rate  来自weather shock
phix = 1.5;         % Elasticity parameter for land utilization cost  calibration 来自weather shcok
tht1	= 0.12;		    % damage current
tht2	= 0;
tht3	= 0;
tht4	= 0;
chi=0.3;           % survival rate of bankers 倒闭率最高的时期，达74%（中国银行总管理处经济研究室编印：《全国银行年鉴》（1937），第7-8页）https://baijiahao.baidu.com/s?id=1776991918221720938&wfr=spider&for=pc
alphab = 0.4;
gammac = 0.6;
phi_pi =-0.1;
phi_fr=1.5;
phi_bG=1.5;
kappa_bG=1;
kappa_pi=1;
kappa_fr=1;

%% steady state of exogenous variables
gg = 0.14;           % public spending over GDP (IPF)% share of public spendings in gdp 资料来源: 数据来自刘䄑, 陭昭: 《近代中国 50 年 GDP 的估算与经济增长研究 (18871936 年), 经济科学出版社 2012 年版, 第 107-108 页; 财政岁入的数据来自表 2-1。和资料来源: 据 [美]阿瑟・恩・杨格: 《1927 至 1937 年中国财政经济情况》,陈泽宪、陈霞飞译，第 483-485 页。
eCss = 1;            % Climate shock
eZss = 1;            % Productivity
eEss = 1;
damage=1;
corrsg=0.8;
%% steady state of endogenous variables and parameters
beta = piss/rss;
rBss = rss/piss + sp;  % lending rate
rkss = rBss - (1-delta);% rental rate of capital
piHss = piss;           % PPI inflation
piXss = piss;           % export price inflation
mcss = (eps-1)/eps;     % real marginal cost
wss = alphah*(alphak/rkss)^(alphak/alphah)*(1-alphak-alphah)^((1-alphak-alphah)/alphah)*mcss^(1/alphah); % real wage
css = wss/hss^varphiH;  % home consumption
gdpss = (css + (delta*alphak*mcss/rkss + (1 - alphak - alphah)*mcss)*(damage - (1-deltaL)))/((1-Abar)*(1-delta*alphak*mcss/rkss - (1-alphak-alphah)*mcss) + Abar - gg);
Gss = gg*gdpss; 
yAss = Abar * gdpss;
yss = gdpss + damage - (1-deltaL);
yNss = yss - yAss; 
ass = (1-alphak-alphah)*mcss*yNss;
varphi = (yAss - ass)/css; 
kss = alphak*mcss/rkss*yNss;
iss = delta * kss; 
hNss = alphah*mcss/wss*yNss;
hAss = hss - hNss; 
omega = 1 - wss*hAss/yAss;
kappaA = yAss^(1/(1-omega))/((damage*lss)^(omega/(1-omega))*eZss*hAss);
lambdaLss = beta*omega/(1-beta)*yAss/lss;
tau = phix^(1-phix)*(damage-(1-deltaL))^(1-phix)*(beta*omega/(1-beta)*damage*yAss)^(-phix);
xss = (phix/tau*(damage-(1-deltaL)))^(1/phix);
nss = kss/phiss; 
iotab = (1-chi*rss/piss)/phiss - chi*sp;
Omegass = (1-chi)/(1-beta*chi*((rBss - rss/piss)*phiss + rss/piss));
muss = beta*Omegass * (rBss - rss/piss);
nuss = beta*Omegass * rss/piss;
theta = muss + nuss/phiss;
pizss = piss;
rzss = rrzss*pizss;
betaz = pizss/rzss;
imss = (varphi*alphaA + (1-varphi)*alphaN)*css;
xpss = (varphi*alphaA + (1-varphi)*alphaN)*css;
yXss = (1-varphi)*alphaN*css;
yHss = yNss - yXss;
%mss = ((1-1/rss)/css)^(-varphiL);
   mss= ((1 - (1 / rss)) / css)^(-1 / varphiL);
pwgss = 3;
pwsss = 1;
pdsss = 1;
ess = pwgss/pwsss;
bGss = bbar * gdpss;
Tax = Gss+bGss*(rss/piss - 1);
%FRss = (mss/gammac-bGss)/pdsss;
FRinss = xpss/pwsss;
FRoutss = imss/pwsss;
FRss= 1;
buffer=0.8;
phi_pds =1.5;
sss =1;
phi_s=1.5;
rho_FR=0.9;
%% Other parameters
  kappaI = 2;        %Smets and Wouters (2007):1.5; Christiano, Eichenbaum, and Evans (2005):1.6
kappaZ = 1;
kappaP = 10;
muN = 1.5;          % elasticity of substitution between domestic and foreign goods 
muA	= 0.2;		% substitution agricultural
muw	= 1.2;		% substitution non-agriculture vs agriculture

% shock persistence 
rho_c	= 0.8;	
rho_g   = 0.95; 
rho_eC   = 0.38; 
rho_eZ = 0.95; 
rho_eE   = 0.31; 
rho_pds = 0.95;
rho_pws = 0.95;
rho_pwg = 0.95;

% shock std
sig_c	= 5; 
sig_g	= 5; 
sig_eC	= 5;
sig_eZ	= 5;
sig_eE	= 5;
sig_pds = 0.95;
sig_pws = 0.95;
sig_pwg = 0.95;
sig_FR =0.95;
% measure error
eps_y_ME= 0.1;
eps_c_ME= 0.1;
eps_im_ME= 0.1;
eps_xp_ME= 0.1;
eps_exch_ME= 0.1;
eps_smdi_ME= 0.1;

% 设置初始条件
target_reduction = 0.6; % 目标剩余产量百分比，即减少40%
shock_value = 2.215;     % shock 大于1的值 

% 计算 tht1 使得在给定 shock 值下的 damage 减少至 0.6
tht1 = -log(target_reduction) / log(shock_value);

% 假设 deltaL, yAss, gdpss, yNss, yss 已经定义并赋值
% 打印 deltaL
fprintf('deltaL = %.4f\n', deltaL);

% 计算并打印 yAss/gdpss
yAss_gdpss_ratio = yAss / gdpss;
fprintf('yAss / gdpss = %.4f\n', yAss_gdpss_ratio);

% 计算并打印 yNss/gdpss
yNss_gdpss_ratio = yNss / gdpss;
fprintf('yNss / gdpss = %.4f\n', yNss_gdpss_ratio);

% 计算并打印 yAss/yss
yAss_yss_ratio = yAss / yss;
fprintf('yAss / yss = %.4f\n', yAss_yss_ratio);

% 计算并打印 yNss/yss
yNss_yss_ratio = yNss / yss;
fprintf('yNss / yss = %.4f\n', yNss_yss_ratio);

% 输出结果
fprintf('To achieve a 40%% reduction in yield with eC = %.2f, tht1 should be set to approximately %.2f.\n', shock_value, tht1);
save param_and_ss alphak alphah beta betaz delta theta iotab chi ...
    varphiH varphiL eps omega phix  deltaL varphi tht1 tht2 tht3 tht4 alphaN muN ...
    alphaA muA muw kappaI kappaZ kappaP kappaA tau alphab gammac ...
    rho_c rho_g rho_eC rho_eZ rho_eE rho_pds rho_pws rho_pwg ...
    sig_c sig_g sig_eC sig_eZ sig_eE sig_pds sig_pws sig_pwg ...
    eps_y_ME eps_c_ME eps_im_ME eps_xp_ME eps_exch_ME eps_smdi_ME ...
    piss rss lss phiss hss rrzss gg eCss eZss eEss rBss rkss piHss ...
    piXss mcss wss css gdpss Gss yAss yss yNss ass kss iss hAss hNss ...
    lambdaLss xss nss Omegass muss nuss pizss rzss imss xpss yHss ...
    yXss mss pwgss pwsss pdsss ess bGss Tax  FRinss FRoutss phi_pi FRss buffer  phi_fr kappa_pi kappa_fr phi_bG kappa_bG sig_FR phi_pds sss phi_s rho_FR bbar corrsg;
