%% set inputs for solution below
%  The program produces responses to the shocks selected below under
%  irfshock. Paths for all the endogenous variables in the model selected
%  are produced. VariableName_difference holds the piece-wise linear solution for
%  VariableName.  VariableName_uncdifference holds the linear solution for
%  VariableName.
clear
clc

global M_ oo_

% modname below chooses model
% directory. But simple param choices are made from paramfile_borrcon in current
% directory.
modnam = 'model_benchmark';
modnamstar = 'model_benchmark_notbinding';

% Enter parameter values and save them in PARAM_EXTRA.mat file.
% This choices overwrite any choice that are made in the paramfile_borrcon file.
betap = 0.9943;                     %patient hh discount factor (GM 2023)
etap = 1;                           %patient hh Frisch labour elasticity (SW 2021 - GM 2023)
betae = 0.975;                      %entrepreneur discount factor (GM 2023)
deltak = 0.025;                     %depreciation of capital (SW 2021 - GM 2023)
alphae = 0.2;                       %capital share (GM 2023)
mc = 0.3;                           %entrepreneur LTV ratio CB
gammap = 0;
phip = 0.87;                        %price stickiness (estimated in GM 2023)
epsilonp = 6;                       %price markup (GM 2023)
kappai = 3.98;                      %investment adjustment costs (estimated in GM 2023)
rhor = 0.88;                        %TR interest rate persistence (estimated in GM 2023)
phipi = 1.87;                       %TR response to inflation (estimated in GM 2023)
phiy = 0.24;                        %TR response to output (estimated in GM 2023)
rhoa = 0.95;                        %persistence productivity shock (SW 2021)
sa = 0.0065;                        %stdev productivity shock (SW 2021)
sr = 0.01;                          %stdev mon policy shock 
betab = 0.965;
gammab = 0.9; 
save PARAM_EXTRA betap etap betae deltak alphae mc gammap ...
    phip epsilonp kappai rhor phipi phiy rhoa sa sr betab gammab

% see notes 1 and 2 in README. Note that 
% lb_ss is the steady-state value of the multiplier lb (see borrcon_steadystate.m).
%    lb_ss is calculated automatically as additional auxiliary parameter
%    around line 49 of solve_one_constraint
constraint = 'deltab < -deltab_ss';
constraint_relax ='dpc > gammab*bce';

% Pick innovation for IRFs
irfshock =char('ea');      % label for innovation for IRFs

shockssequence(3,1) = -0.03;

nperiods = 60;

    
[zdatalinear zdatapiecewise zdatass oobase_ Mbase_  ] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  shockssequence,irfshock,nperiods);


% unpack the IRFs
for i=1:Mbase_.endo_nbr
  eval([deblank(Mbase_.endo_names(i,:)),'_l=zdatalinear(:,i);']);
  eval([deblank(Mbase_.endo_names(i,:)),'_p=zdatapiecewise(:,i);']);
  eval([deblank(Mbase_.endo_names(i,:)),'_ss=zdatass(i);']);
end
   

titlelist = char('c (consumption)','b (borrowing)','y (income)','lb (multiplier)');
percent = 'Percent';
level = 'Level';
ylabels = char(percent,percent,percent,level);
figtitle = 'Simulated variables';
legendlist = cellstr(char('Piecewise Linear','Linear'));

line1=100*[c_p/c_ss,b_p/b_ss,y_p/y_ss,(lb_p+lb_ss)/100];
line2=100*[c_l/c_ss,b_l/b_ss,y_l/y_ss,(lb_l+lb_ss)/100 ];

makechart(titlelist,legendlist,figtitle,ylabels,line1,line2);

 
