%% Ramsey growth model with search unemployment
%
% Last edited by Ben on 10 October 2015

% Endogenous labour supply (unit elastic case)
% Parameters as chosen below Table 13.6 in FoMM3

clear
clf
clc
close all

% Structural parameters
PAR.alpha   = 1/3;              % Efficiency parameter capital
PAR.epsilon = 0.183;            % Utility parameter of consumption
PAR.delta_k = (1.1)^(1/4)-1;    % Quarterly depreciation rate for capital
PAR.delta_m = 0.05;             % Quarterly job destruction rate
PAR.rho     = (1.065)^(1/4)-1;  % Quarterly rate of time preference Shi-Wen
PAR.phi     = 0.6;              % Efficiency parameter V in matching function Shi-Wen
PAR.beta    = 0.4;              % Bargaining power worker Shi-Wen
                                % Hosios beta = 1 - phi holds

% Exogenous variables
EXO.G     = 0.20;             % Government consumption
EXO.s_U   = 0.00;             % Set in calibration
EXO.Z_m   = 1.00;             % Constant in the matching function
EXO.time  = linspace(0,200,2001); % Post-shock time (for plotting) 

% Target variables
Y0        = 1.00;                       % Base line per capita output per quarter 
L0        = 0.20;                       % Base line employment
u0        = 0.06;                       % Base line unemployment rate
R0        = PAR.rho;
KYratio0  = PAR.alpha / (R0+PAR.delta_k); % Base line capital-output ratio
K0        = KYratio0 * Y0;
I0        = PAR.delta_k * K0;
PAR.Z_0   = K0^(-PAR.alpha) * L0^(PAR.alpha-1) * Y0;  
F_L0      = (1-PAR.alpha) * PAR.Z_0 * (K0/L0)^(PAR.alpha); 
EXO.s_U   = 0.4 * F_L0 ;                % Base line search subsidy
S0        = u0 * L0 /(1-u0);
f0        = PAR.delta_m * L0/S0;        % Base line job finding rate
w0        = (PAR.beta *(PAR.rho+PAR.delta_m+f0) * F_L0 ...
     +(1-PAR.beta) * (PAR.rho+PAR.delta_m) * EXO.s_U) ...
     / (PAR.rho+PAR.delta_m+PAR.beta * f0) ; 
MRS0      =  (PAR.beta * f0 * F_L0 ...
     + (PAR.rho+PAR.delta_m) * EXO.s_U) ...
     / (PAR.rho+PAR.delta_m+PAR.beta * f0) ;  
eta_L0    =  (PAR.beta * (F_L0 - EXO.s_U) ) ...
     / (PAR.rho+PAR.delta_m+PAR.beta * f0) ;  
mu_L0     =  ((1-PAR.beta) * (F_L0 - EXO.s_U) ) ...
     / (PAR.rho+PAR.delta_m+PAR.beta * f0) ;   
ctheta0   = f0 * mu_L0 ; 
C0        = Y0 - I0 - ctheta0 * S0 - EXO.G ;  
PAR.epsilon = C0 / (C0 + MRS0 *(1-S0-L0) );
theta0    = (f0/EXO.Z_m)^(1/PAR.phi);
PAR.c     = ctheta0 / theta0 ;


% dydx(1,1) = (PAR.alpha * PAR.Z_0 * K0^(PAR.alpha-1) * L0^(1-PAR.alpha) ...
%     - PAR.delta_k - PAR.rho) * C0 ;
% dydx(2,1) = PAR.Z_0 * K0^(PAR.alpha) * L0^(1-PAR.alpha) - C0 - EXO.G ...
%     - PAR.c * theta0 * S0 - PAR.delta_k * K0 ;
% dydx(3,1) = EXO.Z_m * theta0^PAR.phi * S0 - PAR.delta_m * L0 ;
% dydx(4,1) = (PAR.alpha * PAR.Z_0 * K0^(PAR.alpha-1) * L0^(1-PAR.alpha) - ...
%     PAR.delta_k + PAR.delta_m) * eta_L0 + PAR.beta ...
%     *(EXO.s_U + EXO.Z_m * theta0^PAR.phi * eta_L0 ...
%     - (1-PAR.alpha) * PAR.Z_0 * K0^(PAR.alpha) * L0^(-PAR.alpha) );
% dydx(5,1) = theta0 - ((1-PAR.beta) * EXO.Z_m * eta_L0 / (PAR.beta * PAR.c))^(1/(1-PAR.phi)) ;
% dydx(6,1) = -S0 + 1 - L0 - ((1-PAR.epsilon)/PAR.epsilon) * C0  ...
%       / (EXO.s_U + EXO.Z_m * theta0^PAR.phi * eta_L0) ;


% Compute the steady-state equilibrium

x = zeros(8,1);
[x(:,1),fval,exitflag,output] = fsolve(@StStMacroModelSearch,[C0,S0,L0,K0,eta_L0,mu_L0,theta0,R0]'...
    ,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-8,'TolFun',1.0e-8),PAR,EXO);
                                     
fval;
output;
x;
VAR.Cbase     = x(1,1);
VAR.Sbase     = x(2,1);
VAR.Lbase     = x(3,1);
VAR.Kbase     = x(4,1);
VAR.eta_Lbase = x(5,1);
VAR.mu_Lbase  = x(6,1);
VAR.thetabase = x(7,1);
VAR.Rbase     = x(8,1);
VAR.Ybase     = PAR.Z_0 * VAR.Kbase^PAR.alpha * VAR.Lbase^(1-PAR.alpha);
VAR.Wbase     = PAR.beta *(1-PAR.alpha) * VAR.Ybase / VAR.Lbase ...
    +(1-PAR.beta) * ((1-PAR.epsilon)/PAR.epsilon) * VAR.Cbase ...
    / (1 - VAR.Sbase - VAR.Lbase) ;
VAR.ubase     = 100* VAR.Sbase / (VAR.Sbase+VAR.Lbase); 
PAR.om_LLbase = (1 - VAR.Sbase - VAR.Lbase) / (VAR.Sbase+VAR.Lbase) ; 
VAR.fbase     = PAR.delta_m * VAR.Lbase/VAR.Sbase;
VAR.MRSbase   = ((1-PAR.epsilon)/PAR.epsilon) * VAR.Cbase / (1 - VAR.Sbase - VAR.Lbase);



% Jacobian matrix linearized system

% Delta      = zeros(5,5);
% Delta(1,2) = -(VAR.Cbase/VAR.Kbase)* (VAR.Ybase/VAR.Kbase) * PAR.alpha*(1-PAR.alpha)  ;
% Delta(1,3) =  (VAR.Cbase/VAR.Kbase)* (VAR.Ybase/VAR.Kbase) * PAR.alpha*(1-PAR.alpha) ;
% Delta(2,1) = -1 + PAR.om_LLbase * PAR.c * VAR.thetabase * VAR.Sbase ...
%     / (VAR.ubase * VAR.Cbase) ;
% Delta(2,2) = PAR.alpha * VAR.Ybase / VAR.Kbase - PAR.delta_k ;
% Delta(2,3) = (1-PAR.alpha) * VAR.Ybase / VAR.Lbase + (1-VAR.ubase) * PAR.c * VAR.thetabase * VAR.Sbase ...
%     / (VAR.ubase * VAR.Lbase) ;
% Delta(2,4) = -  PAR.om_LLbase * PAR.c * VAR.thetabase * VAR.Sbase * VAR.fbase ...
%     / (VAR.ubase * VAR.MRSbase) ;
% Delta(2,5) = -  PAR.phi * VAR.eta_Lbase * PAR.om_LLbase * PAR.c * VAR.thetabase * VAR.Sbase * VAR.fbase ...
%     / ((1-PAR.phi) * VAR.ubase * VAR.MRSbase * VAR.mu_Lbase) ...
%     - PAR.c * VAR.thetabase * VAR.Sbase /((1-PAR.phi) * VAR.mu_Lbase);
% Delta(3,1) = - PAR.delta_m * VAR.Lbase * PAR.om_LLbase / (VAR.ubase * VAR.Cbase) ;
% Delta(3,3) = - PAR.delta_m * (1-VAR.ubase) / VAR.ubase - PAR.delta_m ;
% Delta(3,4) = PAR.delta_m * VAR.Lbase * PAR.om_LLbase * VAR.fbase ...
%     / (VAR.ubase * VAR.MRSbase) ;
% Delta(3,5) = PAR.phi * PAR.delta_m * VAR.Lbase * PAR.om_LLbase * VAR.fbase * VAR.eta_Lbase ...
%     / ((1-PAR.phi) * VAR.ubase * VAR.MRSbase * VAR.mu_Lbase) ...
%     + PAR.delta_m * VAR.Lbase * PAR.phi /((1-PAR.phi) * VAR.mu_Lbase);
% Delta




% Fiscal policy shock in the macro model: G from 0.2 to 0.25
EXO.G  = 0.25 ;

x = zeros(8,1);
[x(:,1),fval,exitflag,output] = fsolve(@StStMacroModelSearch,[C0,S0,L0,K0,eta_L0,mu_L0,theta0,R0]'...
    ,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-8,'TolFun',1.0e-8),PAR,EXO);
                                     
fval;
output;
x;
VAR.CGinc     = x(1,1);
VAR.SGinc     = x(2,1);
VAR.LGinc     = x(3,1);
VAR.KGinc     = x(4,1);
VAR.eta_LGinc = x(5,1);
VAR.mu_LGinc  = x(6,1);
VAR.thetaGinc = x(7,1);
VAR.RGinc     = x(8,1);
VAR.YGinc     = PAR.Z_0 * VAR.KGinc^PAR.alpha * VAR.LGinc^(1-PAR.alpha);
VAR.WGinc     = PAR.beta *(1-PAR.alpha) * VAR.YGinc / VAR.LGinc ...
    +(1-PAR.beta) * ((1-PAR.epsilon)/PAR.epsilon) * VAR.CGinc ...
    / (1 - VAR.SGinc - VAR.LGinc) ;
VAR.uGinc     = 100* VAR.SGinc / (VAR.SGinc+VAR.LGinc); 

EXO.G = 0.2 ;

% Search subsidy shock in the macro model: s_U up by 10 percent
EXO.s_U  = 1.10 * EXO.s_U ;

x = zeros(8,1);
[x(:,1),fval,exitflag,output] = fsolve(@StStMacroModelSearch,[C0,S0,L0,K0,eta_L0,mu_L0,theta0,R0]'...
    ,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-8,'TolFun',1.0e-8),PAR,EXO);
                                     
fval;
output;
x;
VAR.CSinc     = x(1,1);
VAR.SSinc     = x(2,1);
VAR.LSinc     = x(3,1);
VAR.KSinc     = x(4,1);
VAR.eta_LSinc = x(5,1);
VAR.mu_LSinc  = x(6,1);
VAR.thetaSinc = x(7,1);
VAR.RSinc     = x(8,1);
VAR.YSinc     = PAR.Z_0 * VAR.KSinc^PAR.alpha * VAR.LSinc^(1-PAR.alpha);
VAR.WSinc     = PAR.beta *(1-PAR.alpha) * VAR.YSinc / VAR.LSinc ...
    +(1-PAR.beta) * ((1-PAR.epsilon)/PAR.epsilon) * VAR.CSinc ...
    / (1 - VAR.SSinc - VAR.LSinc) ;
VAR.uSinc     = 100* VAR.SSinc / (VAR.SSinc+VAR.LSinc); 

EXO.s_U  = EXO.s_U / 1.1 ;

% Matching productivity shock in the macro model: Z_m from 1.0 to 1.1

EXO.Z_m  = 1.10 * EXO.Z_m ;

x = zeros(8,1);
[x(:,1),fval,exitflag,output] = fsolve(@StStMacroModelSearch,[C0,S0,L0,K0,eta_L0,mu_L0,theta0,R0]'...
    ,optimset('MaxFunEvals',1e6,'MaxIter',1000,'TolX',1.0e-8,'TolFun',1.0e-8),PAR,EXO);
                                     
fval;
output;
x;
VAR.CZinc     = x(1,1);
VAR.SZinc     = x(2,1);
VAR.LZinc     = x(3,1);
VAR.KZinc     = x(4,1);
VAR.eta_LZinc = x(5,1);
VAR.mu_LZinc  = x(6,1);
VAR.thetaZinc = x(7,1);
VAR.RZinc     = x(8,1);
VAR.YZinc     = PAR.Z_0 * VAR.KZinc^PAR.alpha * VAR.LZinc^(1-PAR.alpha);
VAR.WZinc     = PAR.beta *(1-PAR.alpha) * VAR.YZinc / VAR.LZinc ...
    +(1-PAR.beta) * ((1-PAR.epsilon)/PAR.epsilon) * VAR.CZinc ...
    / (1 - VAR.SZinc - VAR.LZinc) ;
VAR.uZinc     = 100* VAR.SZinc / (VAR.SZinc+VAR.LZinc); 


EXO.Z_m  = EXO.Z_m / 1.10 ;



% Model parameters
    disp(' ')
    disp('Model parameters')
    disp(sprintf('      Efficiency parameter capital (alpha):                   %12.8f',PAR.alpha))
    disp(sprintf('      Scale parameter production function (Z_0):              %12.8f',PAR.Z_0))
    disp(sprintf('      Quarterly capital depreciation rate (delta_k):          %12.8f',PAR.delta_k))
    disp(sprintf('      Preference for consumption (epsilon):                   %12.8f',PAR.epsilon))
    disp(sprintf('      Vacancy posting cost (c):                               %12.8f',PAR.c))
    disp(sprintf('      Vacancy parameter in matching (phi):                    %12.8f',PAR.phi))
    disp(sprintf('      Job destruction rate (delta_m):                         %12.8f',PAR.delta_m))
    disp(sprintf('      Bargaining power of the worker (beta):                  %12.8f',PAR.beta))


    
    disp(' ')    
    disp('Exogenous variables')
    disp(sprintf('      Government consumption (G):                             %12.8f',EXO.G))
    disp(sprintf('      Search subsidy (s_U):                                   %12.8f',EXO.s_U))


% Implied calibration
    disp(' ')
    disp('Steady-state equilibrium')
    disp(sprintf('      Consumption (C):                                        %12.8f',VAR.Cbase))
    disp(sprintf('      Capital stock (K):                                      %12.8f',VAR.Kbase))
    disp(sprintf('      Employment (L):                                         %12.8f',VAR.Lbase))
    disp(sprintf('      Search time (S):                                        %12.8f',VAR.Sbase))
    disp(sprintf('      Tightness indicator (theta=V/S):                        %12.8f',VAR.thetabase))
    disp(sprintf('      Output (Y):                                             %12.8f',VAR.Ybase))
    disp(sprintf('      Investment (I):                                         %12.8f',PAR.delta_k*VAR.Kbase))
    disp(sprintf('      Interest rate (R):                                      %12.8f',VAR.Rbase))
    disp(sprintf('      Wage rate (W):                                          %12.8f',VAR.Wbase))
    disp(sprintf('      Unemployment rate (S/(S+L)):                            %12.8f',VAR.ubase))
    disp(sprintf('      Worker job value (eta_L):                               %12.8f',VAR.eta_Lbase))
    disp(sprintf('      Worker job value (mu_L):                                %12.8f',VAR.mu_Lbase))
    disp(sprintf('      Relative job value (mu_L/eta_L):                        %12.8f',VAR.mu_Lbase/VAR.eta_Lbase))

    disp(' ')
    disp('Government consumption increase from G=0.2 to G=0.25')
    disp(sprintf('      Consumption (C):                                        %12.8f',VAR.CGinc))
    disp(sprintf('      Capital stock (K):                                      %12.8f',VAR.KGinc))
    disp(sprintf('      Employment (L):                                         %12.8f',VAR.LGinc))
    disp(sprintf('      Search time (S):                                        %12.8f',VAR.SGinc))
    disp(sprintf('      Tightness indicator (theta=V/S):                        %12.8f',VAR.thetaGinc))
    disp(sprintf('      Output (Y):                                             %12.8f',VAR.YGinc))
    disp(sprintf('      Investment (I):                                         %12.8f',PAR.delta_k*VAR.KGinc))
    disp(sprintf('      Interest rate (R):                                      %12.8f',VAR.RGinc))
    disp(sprintf('      Wage rate (W):                                          %12.8f',VAR.WGinc))
    disp(sprintf('      Unemployment rate (S/(S+L)):                            %12.8f',VAR.uGinc))
    disp(sprintf('      Worker job value (eta_L):                               %12.8f',VAR.eta_LGinc))
    disp(sprintf('      Worker job value (mu_L):                                %12.8f',VAR.mu_LGinc))
    disp(sprintf('      Relative job value (mu_L/eta_L):                        %12.8f',VAR.mu_LGinc/VAR.eta_LGinc))


    disp(' ')
    disp('Search subsidy increase by 10 percent')
    disp(sprintf('      Consumption (C):                                        %12.8f',VAR.CSinc))
    disp(sprintf('      Capital stock (K):                                      %12.8f',VAR.KSinc))
    disp(sprintf('      Employment (L):                                         %12.8f',VAR.LSinc))
    disp(sprintf('      Search time (S):                                        %12.8f',VAR.SSinc))
    disp(sprintf('      Tightness indicator (theta=V/S):                        %12.8f',VAR.thetaSinc))
    disp(sprintf('      Output (Y):                                             %12.8f',VAR.YSinc))
    disp(sprintf('      Investment (I):                                         %12.8f',PAR.delta_k*VAR.KSinc))
    disp(sprintf('      Interest rate (R):                                      %12.8f',VAR.RSinc))
    disp(sprintf('      Wage rate (W):                                          %12.8f',VAR.WSinc))
    disp(sprintf('      Unemployment rate (S/(S+L)):                            %12.8f',VAR.uSinc))
    disp(sprintf('      Worker job value (eta_L):                               %12.8f',VAR.eta_LSinc))
    disp(sprintf('      Worker job value (mu_L):                                %12.8f',VAR.mu_LSinc)) 
    disp(sprintf('      Relative job value (mu_L/eta_L):                        %12.8f',VAR.mu_LSinc/VAR.eta_LSinc))
    
    
    disp(' ')
    disp('Matching productivity increase by 10 percent')
    disp(sprintf('      Consumption (C):                                        %12.8f',VAR.CZinc))
    disp(sprintf('      Capital stock (K):                                      %12.8f',VAR.KZinc))
    disp(sprintf('      Employment (L):                                         %12.8f',VAR.LZinc))
    disp(sprintf('      Search time (S):                                        %12.8f',VAR.SZinc))
    disp(sprintf('      Tightness indicator (theta=V/S):                        %12.8f',VAR.thetaZinc))
    disp(sprintf('      Output (Y):                                             %12.8f',VAR.YZinc))
    disp(sprintf('      Investment (I):                                         %12.8f',PAR.delta_k*VAR.KZinc))
    disp(sprintf('      Interest rate (R):                                      %12.8f',VAR.RZinc))
    disp(sprintf('      Wage rate (W):                                          %12.8f',VAR.WZinc))
    disp(sprintf('      Unemployment rate (S/(S+L)):                            %12.8f',VAR.uZinc))
    disp(sprintf('      Worker job value (eta_L):                               %12.8f',VAR.eta_LZinc))
    disp(sprintf('      Worker job value (mu_L):                                %12.8f',VAR.mu_LZinc))   
    disp(sprintf('      Relative job value (mu_L/eta_L):                        %12.8f',VAR.mu_LZinc/VAR.eta_LZinc))
    

    
    

% Solve the model in continuous time using BVP4C 
    
% Shock time:            t_A = 0
% Implementation time:   t_I
% End time:              t_E
% Unanticipated / permanent shock t_I = 0, t_E = inf
% Anticipated   / permanent shock t_I > 0, t_E = inf
% Unanticipated / permanent shock t_I = 0, t_E = finite
% Anticipated   / permanent shock t_I > 0, t_E = finite
    
    CaseVal = 1;
    
    switch CaseVal
        case 1
        % Convergence to the steady state following a five percent
        % reduction in the initial capital stock
        
        EXO.G_new     = EXO.G;
        EXO.t_I       = 0;
        EXO.t_E       = 1e4;
        VAR.K_initial = 0.95 * VAR.Kbase;
        VAR.Keq0      = VAR.K_initial ;           
        VAR.Keq1      = VAR.Kbase ; 
        VAR.Leq0      = VAR.Lbase ;           
        VAR.Leq1      = VAR.Lbase ; 
        VAR.Ceq1      = VAR.Cbase ; 


[T,CpadT,KpadT,LpadT,etaLpadT,time,Cpad,Kpad,Lpad,etaLpad] = UnemploymentRamseyCT(PAR,EXO,VAR) ;
        
        Ypad     = PAR.Z_0 * Kpad.^(PAR.alpha) .* Lpad.^(1-PAR.alpha) ;
        thetapad = ((1-PAR.beta) * EXO.Z_m * etaLpad / (PAR.beta * PAR.c)).^(1/(1-PAR.phi)) ;
        Spad     = 1 - Lpad - ((1-PAR.epsilon)/PAR.epsilon) * Cpad  ...
            ./ (EXO.s_U + EXO.Z_m * thetapad.^PAR.phi .* etaLpad) ;
        Vpad     = thetapad .* Spad ;
        Rpad     = PAR.alpha * (Ypad ./ Kpad) - PAR.delta_k ;
        Ipad     = Ypad - Cpad - PAR.c * thetapad .* Spad - EXO.G_new ;
        Wpad     = PAR.beta * (1-PAR.alpha) * Ypad ./Lpad ...
            + ((1-PAR.beta)*(1-PAR.epsilon)/PAR.epsilon) * Cpad ./ ...
            (1 - Spad - Lpad) ;
        upad     = Spad ./ (Spad + Lpad);
        
        
        figure(1)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('capital')
        plot(time,Kpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_K') 
        pause(1)
        
        figure(2)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('employment')
        plot(time,Lpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_L')
        pause(1)
        
        figure(3)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('Consumption')
        plot(time,Cpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_C')
        pause(1)
             
        figure(4)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('search time')
        plot(time,Spad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_S')
        pause(1)
     
        figure(5)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('job value')
        plot(time,etaLpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_etaL')
        pause(1)
       
        figure(6)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('labour market tightness')
        plot(time,thetapad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_theta')
        pause(1)
        
        figure(7)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('output')
        plot(time,Ypad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_Y')
        pause(1)
        
        figure(8)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('investment')
        plot(time,Ipad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_I')
        pause(1)

        figure(9)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('wage rate')
        plot(time,Wpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_W')
        pause(1)

        figure(10)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('interest rate')
        plot(time,Rpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_R')
        pause(1)
        
        figure(11)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('unemployment rate')
        plot(time,upad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_U')
        pause(1)

        figure(12)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('vacancies')
        plot(time,Vpad,'k','LineWidth',3);
        hold off
        print('-deps2','-loose','F13_11_V')
        pause(1)
        
        
%         clf;
%         subplot(4,2,1)
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('capital')
%         plot(time,Kpad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_K') 
%         
%         subplot(4,2,2)
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('employment')
%         plot(time,Lpad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_L')
%  
%         
%         subplot(4,2,3)
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('Consumption')
%         plot(time,Cpad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_C')
%                 
%        
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('wage')
%         plot(time,Wpad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_W')
%         
%         subplot(4,2,4)
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('interest rate')
%         plot(time,Rpad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_R')
%         
%         subplot(4,2,5)
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('output')
%         plot(time,Ypad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_Y')
%         
%         subplot(4,2,6)
%         hold on
%         box off
%         xlim([0,120])
%         xlabel('post-shock time')
%         ylabel('investment')
%         plot(time,Ipad,'k','LineWidth',3);
%         hold off
%         %print('-deps2','-loose','Figure_UnemploymentRamsey__Case1_I')
%         
%         subplot(4,2,7)
%         

        case 2
        % Unanticipated and permanent increase in G        
        EXO.G_new  = 0.25 ;
        EXO.t_I    = 0;
        EXO.t_E    = 1e4;
        VAR.Keq0   = VAR.Kbase ;    
        VAR.Keq1   = VAR.KGinc ;  
        VAR.Leq0   = VAR.Lbase ;    
        VAR.Leq1   = VAR.LGinc ;  
        VAR.Ceq1   = VAR.CGinc ;  
        

[T,CpadT,KpadT,LpadT,etaLpadT,time,Cpad,Kpad,Lpad,etaLpad] = UnemploymentRamseyCT(PAR,EXO,VAR) ;
        
        Ypad     = PAR.Z_0 * Kpad.^(PAR.alpha) .* Lpad.^(1-PAR.alpha) ;
        thetapad = ((1-PAR.beta) * EXO.Z_m * etaLpad / (PAR.beta * PAR.c)).^(1/(1-PAR.phi)) ;
        Spad     = 1 - Lpad - ((1-PAR.epsilon)/PAR.epsilon) * Cpad  ...
            ./ (EXO.s_U + EXO.Z_m * thetapad.^PAR.phi .* etaLpad) ;
        Vpad     = thetapad .* Spad ;
        Rpad     = PAR.alpha * (Ypad ./ Kpad) - PAR.delta_k ;
        Ipad     = Ypad - Cpad - PAR.c * thetapad .* Spad - EXO.G_new ;
        Wpad     = PAR.beta * (1-PAR.alpha) * Ypad ./Lpad ...
            + ((1-PAR.beta)*(1-PAR.epsilon)/PAR.epsilon) * Cpad ./ ...
            (1 - Spad - Lpad) ;
        upad     = Spad ./ (Spad + Lpad);

        
        figure(1)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('capital')
        plot(time,Kpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_K') 
        pause(1)
        
        figure(2)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('employment')
        plot(time,Lpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_L')
        pause(1)
        
        figure(3)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('Consumption')
        plot(time,Cpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_C')
        pause(1)
             
        figure(4)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('search time')
        plot(time,Spad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_S')
        pause(1)
     
        figure(5)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('job value')
        plot(time,etaLpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_etaL')
        pause(1)
       
        figure(6)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('labour market tightness')
        plot(time,thetapad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_theta')
        pause(1)
        
        figure(7)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('output')
        plot(time,Ypad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_Y')
        pause(1)
        
        figure(8)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('investment')
        plot(time,Ipad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_I')
        pause(1)

        figure(9)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('wage rate')
        plot(time,Wpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_W')
        pause(1)

        figure(10)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('interest rate')
        plot(time,Rpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_R')
        pause(1)

        figure(11)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('unemployment rate')
        plot(time,upad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_U')
        pause(1)     
       
        figure(12)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('vacancies')
        plot(time,Vpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case2_V')
        pause(1)
        
        
        
        case 3
        % Un anticipated and permanent increase in s_U    
        EXO.G_new  = 0.20 ;
        EXO.s_U    = 1.10 * EXO.s_U ;
        EXO.t_I    = 0;
        EXO.t_E    = 1e4;
        VAR.Keq0   = VAR.Kbase ;    
        VAR.Keq1   = VAR.KSinc ;  
        VAR.Leq0   = VAR.Lbase ;    
        VAR.Leq1   = VAR.LSinc ;  
        VAR.Ceq1   = VAR.CSinc ;  
 
[T,CpadT,KpadT,LpadT,etaLpadT,time,Cpad,Kpad,Lpad,etaLpad] = UnemploymentRamseyCT(PAR,EXO,VAR) ;
        
        Ypad     = PAR.Z_0 * Kpad.^(PAR.alpha) .* Lpad.^(1-PAR.alpha) ;
        thetapad = ((1-PAR.beta) * EXO.Z_m * etaLpad / (PAR.beta * PAR.c)).^(1/(1-PAR.phi)) ;
        Spad     = 1 - Lpad - ((1-PAR.epsilon)/PAR.epsilon) * Cpad  ...
            ./ (EXO.s_U + EXO.Z_m * thetapad.^PAR.phi .* etaLpad) ;
        Vpad     = thetapad .* Spad ;
        Rpad     = PAR.alpha * (Ypad ./ Kpad) - PAR.delta_k ;
        Ipad     = Ypad - Cpad - PAR.c * thetapad .* Spad - EXO.G_new ;
        Wpad     = PAR.beta * (1-PAR.alpha) * Ypad ./Lpad ...
            + ((1-PAR.beta)*(1-PAR.epsilon)/PAR.epsilon) * Cpad ./ ...
            (1 - Spad - Lpad) ;       
        upad     = Spad ./ (Spad + Lpad);

        figure(1)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('capital')
        plot(time,Kpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_K') 
        pause(1)
        
        figure(2)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('employment')
        plot(time,Lpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_L')
        pause(1)
        
        figure(3)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('Consumption')
        plot(time,Cpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_C')
        pause(1)
             
        figure(4)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('search time')
        plot(time,Spad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_S')
        pause(1)
     
        figure(5)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('job value')
        plot(time,etaLpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_etaL')
        pause(1)
       
        figure(6)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('labour market tightness')
        plot(time,thetapad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_theta')
        pause(1)
        
        figure(7)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('output')
        plot(time,Ypad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_Y')
        pause(1)
        
        figure(8)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('investment')
        plot(time,Ipad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_I')
        pause(1)

        figure(9)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('wage rate')
        plot(time,Wpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_W')
        pause(1)

        figure(10)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('interest rate')
        plot(time,Rpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_R')
        pause(1)

        figure(11)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('unemployment rate')
        plot(time,upad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_U')
        pause(1)     
       
        figure(12)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('vacancies')
        plot(time,Vpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case3_V')
        pause(1)
        
        case 4
        % Un anticipated and permanent increase in s_U    
        EXO.G_new  = 0.20 ;
        EXO.Z_m    = 1.10 * EXO.Z_m ;
        EXO.t_I    = 0;
        EXO.t_E    = 1e4;
        VAR.Keq0   = VAR.Kbase ;    
        VAR.Keq1   = VAR.KZinc ;  
        VAR.Leq0   = VAR.Lbase ;    
        VAR.Leq1   = VAR.LZinc ;  
        VAR.Ceq1   = VAR.CZinc ;  
 
[T,CpadT,KpadT,LpadT,etaLpadT,time,Cpad,Kpad,Lpad,etaLpad] = UnemploymentRamseyCT(PAR,EXO,VAR) ;
        
        Ypad     = PAR.Z_0 * Kpad.^(PAR.alpha) .* Lpad.^(1-PAR.alpha) ;
        thetapad = ((1-PAR.beta) * EXO.Z_m * etaLpad / (PAR.beta * PAR.c)).^(1/(1-PAR.phi)) ;
        Spad     = 1 - Lpad - ((1-PAR.epsilon)/PAR.epsilon) * Cpad  ...
            ./ (EXO.s_U + EXO.Z_m * thetapad.^PAR.phi .* etaLpad) ;
        Vpad     = thetapad .* Spad ;
        Rpad     = PAR.alpha * (Ypad ./ Kpad) - PAR.delta_k ;
        Ipad     = Ypad - Cpad - PAR.c * thetapad .* Spad - EXO.G_new ;
        Wpad     = PAR.beta * (1-PAR.alpha) * Ypad ./Lpad ...
            + ((1-PAR.beta)*(1-PAR.epsilon)/PAR.epsilon) * Cpad ./ ...
            (1 - Spad - Lpad) ;       
        upad     = Spad ./ (Spad + Lpad);

        figure(1)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('capital')
        plot(time,Kpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_K') 
        pause(1)
        
        figure(2)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('employment')
        plot(time,Lpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_L')
        pause(1)
        
        figure(3)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('Consumption')
        plot(time,Cpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_C')
        pause(1)
             
        figure(4)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('search time')
        plot(time,Spad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_S')
        pause(1)
     
        figure(5)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('job value')
        plot(time,etaLpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_etaL')
        pause(1)
       
        figure(6)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('labour market tightness')
        plot(time,thetapad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_theta')
        pause(1)
        
        figure(7)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('output')
        plot(time,Ypad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_Y')
        pause(1)
        
        figure(8)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('investment')
        plot(time,Ipad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_I')
        pause(1)

        figure(9)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('wage rate')
        plot(time,Wpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_W')
        pause(1)

        figure(10)
        hold on
        box off
        xlim([0,120])
        xlabel('post-shock time')
        ylabel('interest rate')
        plot(time,Rpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_R')
        pause(1)

        figure(11)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('unemployment rate')
        plot(time,upad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_U')
        pause(1)     
       
        figure(12)
        hold on
        box off
        xlim([0,10])
        xlabel('post-shock time')
        ylabel('vacancies')
        plot(time,Vpad,'k','LineWidth',3);
        hold off
        %print('-deps2','-loose','Figure_UnemploymentRamsey__Case4_V')
        pause(1)
       
     end
      
    