clear all
close all
clc

% -----------------------------------------------------------------------
% 1) RUN first "benchmark model''
% -----------------------------------------------------------------------
dynare Model

% -----------------------------------------------------------------------
% 2) Record parameter value and steady state value
%       they will be used to initalize the next model but when the 
%       tax will change. I define new param like "param_"
% -----------------------------------------------------------------------

% PARAMETERS 
betas_=betas;betau_=betau;etas_=etas;etau_=etau;alpha_=alpha;r_=r;
ds_=ds;du_=du;Ns_=Ns;Nu_=Nu;yfu_=yfu;yfs_=yfs;tauws_=tauws;tauwu_=tauwu;
tauf_=tauf;Ai_=Ai;gammas_=gammas;gammau_=gammau;
rho_=rho;xs_=xs;xu_=xu;%z_=z;

% STEADY STATE VARIABLES 
us_bar_=us_bar;  uu_bar_=uu_bar;
uLi_bar_=uLi_bar;thetas_sq_bar_=thetas_sq_bar;
thetau_uq_bar_=thetau_uq_bar;sq_bar_=sq_bar;
uq_bar_=uq_bar;thetas_bar_=thetas_bar;
thetau_bar_=thetau_bar;Lfs_bar_=Lfs_bar;
Lfu_bar_=Lfu_bar;Ms_bar_=Ms_bar;
Mu_bar_=Mu_bar;vs_bar_=vs_bar;
vu_bar_=vu_bar;sw_bar_=sw_bar;
uw_bar_=uw_bar;wi_bar_=wi_bar;
yi_bar_=yi_bar;
zs_bar=zs_bar;zu_bar=zu_bar;

% return
% -----------------------------------------------------------------------
% 3) Simulate the model for different value of the taxes
% Here I use the parameters and the steady state saved above
% Note that here I change the tax tauw 
% -----------------------------------------------------------------------

% a) Here are the vector of taxes (I just consider the first)
%z_vec= [0.20:0.01:0.40];	%benchmark is z=0.30
tauwu_vec =[0.25:0.01:0.35]; %benchmark is tauwu =0.15;
tauws_vec =[0.30:0.01:0.40];%benchmark is tauwu =0.35;
%tauf_vec =[0.25:0.01:0.35];	%benchmark is tauf =0.30;
% b) I initialize the steady state for unemployment which is now
%    a vector whose dimension will be the number of different value
%    for the tax
%us_simul=zeros(length(tauwu_vec),1);
uu_simul=zeros(length(tauwu_vec),1);

%Lfs_simul=zeros(length(tauwu_vec),1);
Lfu_simul=zeros(length(tauwu_vec),1);

%sw_simul=zeros(length(tauwu_vec),1);
uw_simul=zeros(length(tauwu_vec),1);

uLi_simul=zeros(length(tauwu_vec),1);
wi_simul=zeros(length(tauwu_vec),1);

zs_simul=zeros(length(tauwu_vec),1);


%---------------------------------------
us_simul=zeros(length(tauws_vec),1);
%uu_simul=zeros(length(tauws_vec),1);

Lfs_simul=zeros(length(tauws_vec),1);
%Lfu_simul=zeros(length(tauws_vec),1);

sw_simul=zeros(length(tauws_vec),1);
%uw_simul=zeros(length(tauws_vec),1);

%uLi_simul=zeros(length(tauws_vec),1);
%wi_simul=zeros(length(tauws_vec),1);
zu_simul=zeros(length(tauwu_vec),1);

%---------------------------------------

disp('simulation')




for i=1:length(tauwu_vec)
    tauwu_=tauwu_vec(i);
    
    dynare Model1
    
    uu_simul(i)=oo_.steady_state(4);
    
  
    uw_simul(i)=oo_.steady_state(10);
    
    uLi_simul(i)=oo_.steady_state(15);
  
    Lfu_simul(i)=oo_.steady_state(18);
    
    wi_simul(i)=oo_.steady_state(11);
    
    zu_simul(i)=oo_.steady_state(20);
end

 for ii=1:length(tauws_vec)
     tauws_=tauws_vec(ii);
     
     dynare Model1
     
    us_simul(ii)=oo_.steady_state(3);

    
    sw_simul(ii)=oo_.steady_state(9);
  
    
    Lfs_simul(ii)=oo_.steady_state(17);
    zs_simul(ii)=oo_.steady_state(19);
 end
 
figure(1)

%subplot(331), plot(tauwu_vec,us_simul,'o-k'), grid on, title('skilled unemployement');
subplot(331), plot(tauwu_vec,uu_simul,'o-k'), grid on, title('unskilled unemployement');
%subplot(33), plot(tauwu_vec,sw_simul,'o-k'), grid on, title('Formal skilled wage');
subplot(332), plot(tauwu_vec,uw_simul,'o-k'), grid on, title('Formal unskilled wage');
subplot(333), plot(tauwu_vec,wi_simul,'o-k'), grid on, title('Informal wage');
subplot(334), plot(tauwu_vec,uLi_simul,'o-k'), grid on, title('Informal Employment');
%subplot(337), plot(tauwu_vec,Lfs_simul,'o-k'), grid on, title('Skilled Employment');
subplot(335), plot(tauwu_vec,Lfu_simul,'o-k'), grid on, title('Unskilled Employment');
subplot(336), plot(tauwu_vec,zu_simul,'o-k'), grid on, title('Unskilled unemployment insurance');


figure(2)

subplot(221), plot(tauws_vec,us_simul,'o-k'), grid on, title('skilled unemployement');
%subplot(332), plot(tauws_vec,uu_simul,'o-k'), grid on, title('unskilled unemployement');
subplot(222), plot(tauws_vec,sw_simul,'o-k'), grid on, title('Formal skilled wage');
%subplot(334), plot(tauws_vec,uw_simul,'o-k'), grid on, title('Formal unskilled wage');
%subplot(335), plot(tauws_vec,wi_simul,'o-k'), grid on, title('Informal wage');
%subplot(336), plot(tauws_vec,uLi_simul,'o-k'), grid on, title('Informal Employment');
subplot(223), plot(tauws_vec,Lfs_simul,'o-k'), grid on, title('Skilled Employment');
%subplot(338), plot(tauws_vec,Lfu_simul,'o-k'), grid on, title('Unskilled Employment');
subplot(224), plot(tauws_vec,zs_simul,'o-k'), grid on, title('skilled unemployment insurance');

% NOTE: the order of steady state variable in oo_.steady_state is
% given by the order when the variable are declared in the .mod file
% 1 Ms 
% 2 Mu 
% 3 us 
% 4 uu
% 5 vs 
% 6 vu 
% 7 sq 
% 8 uq 
% 9 sw  
% 10 uw  
% 11 wi 
% 12 thetas 
% 13 thetau 
% 14 yi 
% 15 uLi 
% 16 a 
% 17 Lfs 
% 18 Lfu

