% clear all
% close all

% Calibrate or compute parameters in steady-state in Jianfeng Liu, and Zhiwei Xu (2017)
% Calibrate parameters
% Parameters related to utility function
bta = 0.995;                             % Subjective discount factor
am  = 0.06;                              % Utility weight on real money balances 

% Parameters related to household characterize and decision
muu = 1;                                 % Idiosyncratic net-wealth mean, "1" means no distorted resources constraint
sig = 0.954;                             %0.721 Idiosyncratic net-wealth uncertainty constant  
del = 0.025;                             % 0.10 for annual depreciation rate about real phycical capital

% Technology parameters
gA  = 1.02;                              % 1.08 for annual GDP per capita growth rate
A   = 1.00;                              % Technology level in steady-state

% Parameters for nominal ridigities
phi = 10;                                % Elasticity of substitution between differentiated products
Ome = 60;                                % Price adjustment cost parameter 

% Parameters for international trade
omega  = 1.5;                            % External final goods demands elasticity
re     = 1;                              % Real exchange rate or term of trade in steady-state
Rstar  = 1.0175;                         %1.0120 Foreign nominal interest rate in steady-state 
pistar = 1.005;                          % Foreign inflation rate in steady-state 

% Parameters related to monetary policy and capital controls
rhoR   = 0.00;                           % Persistence of original interest rate in Taylor rule is zero
rhopi  = 1.5;                            % Response intensity of inflation rate in Taylor rule
rhoGDP = 0.5;                            %0.125? Response intensity of GDP in Taylor rule
ge     = 1;                              %0.9870? Exchange rate growth in steady-state

% Constant
pipi    = pi;

% Compute the other parameters in steady-state
piss    = pistar;                        % Domestic inflation rate 
r       = gA/bta+del-1;                  % Real physical capital return rate 
R       = gA*piss/bta;                   % Domestic interest rate

CG      = 0.43;                          % Average private comsumpton-to-GDP ratio = 0.43 from 1990 to 2009
EG      = 0.23;                          % Average import-to-GDP ratio = 0.20 from 1990 to 2009
ImG     = 0.20;                          % Average export-to-GDP ratio = 0.23 from 1990 to 2009
IG      = 1-CG-EG+ImG;                   % Physical investment-to-GDP ratio = 0.54 in steady-state 
KG      = IG/(gA-(1-del));               % Physical capital-to-GDP ratio in steady-state
gam     = r*KG/(ImG+r*KG);               % Cost share of physical capital in intermediate goods production
v       = (phi-1)/phi;                   % Real marginal cost 
PIm    = (r/gam)^gam*(re/(1-gam))^(1-gam);                           % Intermediate goods price

Nss       = 0.40;                        % Targeted labor supply as a fraction of time endowment
thestar = Funccapifliendomini(bta,muu,sig,gA,pistar,Rstar,ge,CG,EG,ImG,Nss,del,phi,am); % Cutoff related to idiosyncratic net-wealth uncertainty
sigtil  = (log(1+sig^2/muu^2))^(1/2);    % Define sigtilt (lognormal distribution)
muutil  = log(muu)-1/2*sigtil^2;         % Define muutilt (lognormal distribution)
F1      = normcdf((log(thestar)-muutil)/sigtil);                      % Define F1t
F2      = normcdf((log(thestar)-(muutil+sigtil^2))/sigtil);           % Define F2t
Phi     = thestar*F1+muu*(1-F2);         % Unertainty premium related to idiosyncratic net-wealth 
Psi     = -Phi+thestar+muu;              % Define Psit
The     = muu-Psi;                       % Define Thet 
Prob    = F1;                            % Binding probabilility of foreign bonds holding contraint 

tau     = Phi*Rstar*ge/R;                % Capital controls degree (or tax) parameter
alf     = tau*Psi*(r*KG+ImG)/((tau-1)*The*v*CG+tau*Psi*v*(1+ImG));    % Cost share of intermediate goods in retail goods production
an      = (1-alf)/alf*Psi*Phi/(thestar*Nss)*(r*KG+ImG)/CG;            % Disutility weight on labor in utility function

w       = (1-alf)*A*v^(1/(1-alf))*(alf/PIm)^(alf/(1-alf));            % Real wage rate 
Y       = w*Nss/((1-alf)*v);             % Total final goods outputs 
GDP     = alf*v*Y/(r*KG+ImG);            % Real GDP
Ex      = GDP*EG;                        % Export in steady-state
Exstar  = Ex*(re)^(-omega)/A;            % External final goods demands parameter

%sig = 1.25*0.954;                        % 0.721 Idiosyncratic net-wealth uncertainty constant
%sig = 1.50*0.954;                        % 0.721 Idiosyncratic net-wealth uncertainty constant

%tau     = 1.025*Phi*Rstar*ge/R;           % Capital controls degree (or tax) parameter
%tau     = 1.050*Phi*Rstar*ge/R;           % Capital controls degree (or tax) parameter
%tau     = 1.075*Phi*Rstar*ge/R;           % Capital controls degree (or tax) parameter

%tau     = 1.10*Phi*Rstar*ge/R;            % Capital controls degree (or tax) parameter
%tau     = 5.00*Phi*Rstar*ge/R;            % Capital controls degree (or tax) parameter

save ParaSet bta an am muu sig del gA alf gam A phi Ome omega re Rstar Exstar pistar rhoR rhopi rhoGDP ge tau pipi
dynare Benchcapifliendom
%dynare Benchcapifliendom0
save SteadyStateValue GDPss Css Iss bgstarss bpstarss bstarss gess piss Rss gMss ress thestarss Phiss sigss cass Nss Ass Yss Kss rss wss Xss Psiss Thess Probss Exss Exstarss Rstarss Imss Gamss PImss mss vss WFstarss WFss

DelFiles

%set(gcf,'PaperPositionMode','auto')
%print  -dtiff IRF.tiff

% GDPt Ct It bgstart bpstart bstart get pit Rt gMt thestart Phit ret Imt Ext cat Nt WFt
figure (1)
T=1:20;
subplot(3,6,1)
plot(T,GDPt_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('GDP_{t}')
subplot(3,6,2)
plot(T,Ct_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('C_{t}')
subplot(3,6,3)
plot(T,It_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('I_{t}')
subplot(3,6,4)
plot(T,bgstart_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('b^{*}_{gt}')
subplot(3,6,5)
plot(T,bpstart_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('b^{*}_{pt}')
subplot(3,6,6)
plot(T,bstart_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('b^{*}_{t}')

subplot(3,6,7)
plot(T,get_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('g_{et}')
subplot(3,6,8)
plot(T,pit_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('\pi_{t}')
subplot(3,6,9)
plot(T,Rt_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('R_{t}')
subplot(3,6,10)
plot(T,gMt_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('g_{Mt}')
subplot(3,6,11)
plot(T,thestart_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('\theta^{*}_{t}')
subplot(3,6,12)
plot(T,Phit_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('\Phi_{t}')

subplot(3,6,13)
plot(T,ret_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('re_{t}')
subplot(3,6,14)
plot(T,Imt_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('Im_{t}')
subplot(3,6,15)
plot(T,Ext_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('Ex_{t}')
subplot(3,6,16)
plot(T,cat_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('ca_{t}')
subplot(3,6,17)
plot(T,Nt_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('N_{t}')
subplot(3,6,18)
plot(T,WFt_esigt,'-r','LineWidth',2.5); hold on; grid on;
title('WF_{t}')

%set(gcf,'PaperPositionMode','auto')
%print  -dtiff IRF1_esigt.tiff
