function [T,CpadT,KpadT,LpadT,etaLpadT,time,Cpad,Kpad,Lpad,etaLpad] = UnemploymentRamseyCT(PAR,EXO,VAR)

% Temporary shocks (unanticipated or anticipated)
% Announcement time of the policy    t_A = 0
% Implementation time of the policy  t_I >= 0
% End time of the policy  t_E >= t_I

% Parameter values:

alpha   = PAR.alpha;           
epsilon = PAR.epsilon;
delta_k = PAR.delta_k;    
delta_m = PAR.delta_m;             
rho     = PAR.rho; 
phi     = PAR.phi;             
beta    = PAR.beta;  
Z_0     = PAR.Z_0;
c       = PAR.c;


% Exogenous variables:
G_old   = EXO.G;
time    = EXO.time;
G_new   = EXO.G_new ;
Z_m     = EXO.Z_m;
s_U     = EXO.s_U; 

% Capital and employment
Keq0    = VAR.Keq0;
Keq1    = VAR.Keq1;
Leq0    = VAR.Leq0;
Leq1    = VAR.Lbase;
Ceq1    = VAR.Ceq1;


% Initial consumption and eta_L
Ceq0    = VAR.Cbase;
etaL0   = VAR.eta_Lbase;

% Policy parameters

t_I     = EXO.t_I;
t_E     = EXO.t_E;


solinit  = bvpinit(linspace(0,400),[Ceq0 Keq0 Leq0 etaL0]);
sol      = bvp4c(@capck,@bcfun,solinit);
T        = sol.x;
CpadT    = sol.y(1,:);
KpadT    = sol.y(2,:);
LpadT    = sol.y(3,:);
etaLpadT = sol.y(4,:);

y       = deval(sol,time);
Cpad    = y(1,:);
Kpad    = y(2,:);
Lpad    = y(3,:);
etaLpad = y(4,:);

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


% Dynamics of C, K, L, and eta_L

function dydx = capck(x,y)
        if x >= t_I & x < t_E 
            G = G_new;
        else
            G = G_old;
        end
        
theta = ((1-beta) * Z_m * y(4) / (beta * c))^(1/(1-phi)) ;
S     =  1 - y(3) - ((1-epsilon)/epsilon) * y(1)  ...
      / (s_U + Z_m * theta^phi * y(4)) ;        
dydx(1,1) = (alpha * Z_0 * y(2)^(alpha-1) * y(3)^(1-alpha) ...
    - delta_k - rho) * y(1) ;
dydx(2,1) = Z_0 * y(2)^(alpha) * y(3)^(1-alpha) - y(1) - G ...
    - c * theta * S - delta_k * y(2) ;
dydx(3,1) = Z_m * theta^phi * S - delta_m * y(3) ;
dydx(4,1) = (alpha * Z_0 * y(2)^(alpha-1) * y(3)^(1-alpha) - ...
    delta_k + delta_m) * y(4) + beta ...
    *(s_U + Z_m * theta^phi * y(4) ...
    - (1-alpha) * Z_0 * y(2)^(alpha) * y(3)^(-alpha) );        

end


% Boundary conditions

function res = bcfun(ya,yb)
       res = [ya(2) - Keq0
              yb(2) - Keq1
              ya(3) - Leq0
              yb(1) - Ceq1];
end


end