%% 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 
 global M_ oo_

modnam = 'rbc_testing';
modnamstar = 'rbc_testing_1';

% Constraint (see notes 1 and 2 in readme.pdf file)

%constraint = '(kappa * q * k - (z + (1-delta)*q)* theta * k(-1) + r*d(-1) - D) < 0';
constraint = ('(1-kappa).* (q / k.^(-1)) - d < 0');
constraint_relax = 'mu < 0';

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

%SHOCKS = [ zeros(5,2)...
%   0.01,  0.05...
%  zeros(94,2) ];


shockssequence = zeros(60,1);
shockssequence(10) = 0.02; 
shockssequence(30) = -0.02; 
%shockssequence = SHOCKS;

maxiter = 20;

nperiods=100;
    
[zdatalinear, zdatapiecewise, zdatass, oobase_, Mbase_  ] = ...
  solve_one_constraint(modnam,modnamstar,...
  constraint, constraint_relax,...
  shockssequence,irfshock,nperiods,maxiter);
    
% 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


%define inputs for plotting
titlelist = char('i (investment)','k (capital)','y (output)','l (labor)','c (consumption)','d (deposits)', 'w (wage)', 'D (dividends)', 'r (risk_free)', 'omega (multiplier1)', 'mu (multiplier2)', 're (ret_capital)', 'q (asset_price)', 'z (mpk)', 'lamb (SDF)','ul (mar_ul)', 'uc (mar_uc)', 'ud (mar_ud)', 'A (tfp)', 'theta (fin_shock)');
percent = 'Percent';
level = 'Level';
ylabels = char(percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent,percent);
figtitle = '';
legendlist = cellstr(char('Piecewise Linear','Linear'));
line1=100*[i_p,k_p,y_p,l_p,c_p,d_p,w_p,D_p,r_p,omega_p,mu_p,re_p,q_p,z_p,lamb_p,ul_p,uc_p,ud_p,A_p,theta_p];
line2=100*[i_l,k_l,y_l,l_l,c_l,d_l,w_l,D_l,r_l,omega_l,mu_l,re_l,q_l,z_l,lamb_l,ul_l,uc_l,ud_l,A_l,theta_l];

%create plots
makechart(titlelist,legendlist,figtitle,ylabels,line1,line2)