function [paths,ss0, elasticity]=transition_v09(Parameter, Calibration)

% - solver default options ------------------------------------------------
options=optimset('disp','off','LargeScale','off', ...
                 'TolFun',.00000001,'MaxFunEvals',10^7,'MaxIter',10^7);

%% - Parametrizations -----------------------------------------------------

global ssigma rrho ddelta p ...
       aawtd Lwtd kssf cssf mssf A_neutral A_neutral_0 A_neutral_f ...
       kss0 eeta_f eeta_0 aalpha_0 aalpha_f ell_z psi_L ...
       model time_path h start_year finish_year

% - feed in calibrated values ---------------------------------------------   
A_neutral_0=Calibration.A_0;
A_neutral_f=Calibration.A_f;
eeta_0     =Calibration.eta_0;
eeta_f     =Calibration.eta_f;
aalpha_0   =Calibration.alpha_0;
aalpha_f   =Calibration.alpha_f;
ell_z      =Calibration.ell_z;
psi_L      =Calibration.psi_L;
rk_target  =Calibration.rk_target;
Z=numel(ell_z);

% SWITCH the model: what happens to wealth upon death
% 1: wealth destruction upon death
% 2: one more day to leave: individuals consume all wealth upon death
model = 2; 

% SWITCH for time path of the shocks
%time_path = 'exponential';
%time_path = 'linear';
time_path = 'smoothed_linear';
   
% - model parametrization -------------------------------------------------

ssigma      = Parameter.sigma;   % inverse IES
ddelta      = Parameter.delta;   % depreciation
p           = Parameter.p;

% - calculate cap share
    aawtd_0  = aalpha_0*eeta_0';      
    aawtd_f  = aalpha_f*eeta_f';  

% - calibration of rho ----------------------------------------------------
    eq    = @(x) (-((1-aawtd_0)*ddelta-x-aawtd_0*p*ssigma)+(((1-aawtd_0)*ddelta-x-aawtd_0*p*ssigma)^2+4*(1-aawtd_0)*x*ddelta)^(0.5))/2 - rk_target;
    rrho     = fzero(eq,0.1, options);

% - shock parametrization -------------------------------------------------
start_year = 1980; 
finish_year= 2018; 
h = 15;  % parameter for smoothing when smoothed_linear option chosen (width of the window from: [finish_year-h, finish_year+h] on which we smooth

% - weighted labor input --------------------------------------------------
Lwtd_0   = prod((psi_L.*ell_z./(eeta_0.*(1-aalpha_0))).^(eeta_0.*(1-aalpha_0)));
Lwtd_f   = prod((psi_L.*ell_z./(eeta_f.*(1-aalpha_f))).^(eeta_f.*(1-aalpha_f)));

% - check that etas sum up to 1 ----------------------------------------
if  abs(sum(eeta_0)-1)>0.00001 || abs(sum(eeta_f)-1)>0.00001
    disp('etas do not sum to 1. Cannot proceed')
    return
else 
end 


%% INITIAL STEADY STATE----------------------------------------------------

% set the effective alpha and L to the initial steady state values
aawtd = aawtd_0; 
Lwtd  = Lwtd_0;
A_neutral = A_neutral_0;

% Analytical calcs for initial guess of numerical procedure.
% As showed in the algebra, KSS can be written as a function of auxillary variable 
% X which solves the quadratic aX^2 + bX + c = 0, where X is the marginal
% product of capital. 
a = 1;
b = ((1-aawtd)*ddelta-rrho-aawtd*p*ssigma);
c = -(1-aawtd)*rrho*ddelta;
X = (-b+(b^2-4*a*c)^0.5)/2*a + ddelta; % quadratic formula, positive root
kss0_analytical = aawtd*((A_neutral*Lwtd)/X)^(1/(1-aawtd)); %KSS as a function of X
K = kss0_analytical; %temporary shorthand

% consumption differs between the two models
if model ==1 
css0_analytical = FK(K) - ddelta*K - p*K;
elseif model ==2
css0_analytical = FK(K) - ddelta*K;
end 

mss0_analytical = MPK(K) - ddelta + (1/ssigma)*(rrho+ddelta-MPK(K));
clear K

% Numerical initial SS [This is redundant really, but just in case and to
% check the algebra in the analytical solution. Plus left here for future
% if more complicated problems where no analytical solution possibly
% arise].

fun = @STEADY_STATE;
x0 = [kss0_analytical,css0_analytical,mss0_analytical]; %initial guess
ss_0 = fsolve(fun,x0, options);

kss0 = ss_0(1);
css0 = ss_0(2);
mss0 = ss_0(3);

if kss0 - kss0_analytical > 0.0001
    disp('Analytical and numerical initial steady states differ! Check calculations')
    pause
end 

%Other variables: 
yss0  = A_neutral*((kss0/aawtd).^aawtd)*Lwtd; %output
rss0  = MPK(kss0)-ddelta;
wzss0 = (1-aalpha_0).*yss0.*eeta_0./ell_z;
iss0  = yss0-css0;
kyss0 = kss0/yss0; 
mwss0 = mean(wzss0);
lsss0 = mwss0./yss0;
npv_ss0 = wzss0./rss0;
nkss0 = rss0*kss0 / (rss0*kss0 + mwss0);
ksharess0 = (rss0+ddelta)*kss0/yss0;

%% FINAL STEADY STATE------------------------------------------------------

aawtd = aawtd_f;
Lwtd  = Lwtd_f;
A_neutral = A_neutral_f;

% Analytical - same steps as before
a = 1;
b = ((1-aawtd)*ddelta-rrho-aawtd*p*ssigma);
c = -(1-aawtd)*rrho*ddelta;
X = (-b+(b^2-4*a*c)^0.5)/2*a + ddelta; % quadratic formula, positive root
kssf_analytical = aawtd*((A_neutral*Lwtd)/X)^(1/(1-aawtd)); %KSS as a function of X
K = kssf_analytical; %temporary shorthand

if model ==1 
cssf_analytical = FK(K) - ddelta*K - p*K;
elseif model ==2
cssf_analytical = FK(K) - ddelta*K;
end 

mssf_analytical = MPK(K) - ddelta + (1/ssigma)*(rrho+ddelta-MPK(K));
clear K

% Numerical steady state:
fun = @STEADY_STATE;
x0 = [kssf_analytical,cssf_analytical,mssf_analytical];
x_f = fsolve(fun,x0, options);

kssf = x_f(1);
cssf = x_f(2);
mssf = x_f(3);


if kssf - kssf_analytical > 0.0001
    disp('Analytical and numerical final steady states differ! Check calculations')
    %pause
end 

%Other variables: 
yssf  = A_neutral*((kssf/aawtd).^aawtd)*Lwtd; %output
rssf  = MPK(kssf) - ddelta;
wzssf = (1-aalpha_f).*yssf.*eeta_f./ell_z;
issf  = yssf-cssf;
kyssf = kssf/yssf; 
mwssf = mean(wzssf);
lsssf = mwssf./yssf;
npv_ssf = wzssf./rssf;
nkssf = rssf*kssf / (rssf*kssf + mwssf);
ksharessf = (rssf+ddelta)*kssf/yssf;

%% Slope of capital supply in the steady state-----------------------------

% Derived in closed form, this slope is:

if ddelta == 0
semi_elasticityKS_0 = (1/(1-aawtd_0))*(rrho/((rrho+p*ssigma*aawtd_0)*(p*ssigma*aawtd_0))+1/(p*ssigma*(1-aawtd_0)));
elasticityKS_0 = rss0*semi_elasticityKS_0;
semi_elasticityKS_f = (1/(1-aawtd_f))*(rrho/((rrho+p*ssigma*aawtd_f)*(p*ssigma*aawtd_f))+1/(p*ssigma*(1-aawtd_f)));
elasticityKS_f = rssf*semi_elasticityKS_f;

else 
% Case with positive depreciation
semi_elasticityKS_0 = (1/(1-aawtd_0))*((rrho/rss0)/(rss0-rrho)+1/(p*ssigma+rrho-rss0));
elasticityKS_0 = rss0*semi_elasticityKS_0;
semi_elasticityKS_f = (1/(1-aawtd_f))*((rrho/rssf)/(rssf-rrho)+1/(p*ssigma+rrho-rssf));
elasticityKS_f = rssf*semi_elasticityKS_f;
end 

elasticity.semi=semi_elasticityKS_0;
elasticity.full=elasticityKS_0;


%% Transitional dynamics---------------------------------------------------

if model == 1 

disp('Using relax_percentiles_1') 
cd code/relax_percentiles_1

elseif model == 2
    
disp('Using relax_percentiles_2') 
cd code/relax_percentiles_2 

end 

main

cd ..
cd ..

% convert the key paths and variables into the current notation / dimensions:
time = t';
time_save = time;

kt = k_path';
ct = c_path';
mt = m_path';

%% Construct timepaths for primitives and their change over time-----------

N = numel(time);
A_neutral_t = zeros(N,1);
aalpha_t    = zeros(N,Z);
eeta_t    = zeros(N,Z);

for i = 1:numel(time)
    t = time(i);
    if t==0        % At time -1 the parameters should have the values for the old steady state
    aalpha_t(i,:)   = aalpha_0;
    eeta_t(i,:)   = eeta_0;
    A_neutral_t(i,1)= A_neutral_0;
    
    elseif 0<=t
        
        if strcmpi(time_path, 'exponential')
        
            aalpha_t(i,:)    = aalpha_0'*exp(-aalpha_rate*(t+0.001))+aalpha_f'*(1-exp(-aalpha_rate*(t+0.001))); % the +0.001 here is needed to recreate an instantenous shock at t=0 when aalpha_rate = inf.
            eeta_t(i,:)    = eeta_0'*exp(-eeta_rate*(t+0.001))+eeta_f'*(1-exp(-eeta_rate*(t+0.001)));
            A_neutral_t(i,1) = A_neutral_0'*exp(-A_rate*(t+0.001))+A_neutral_f'*(1-exp(-A_rate*(t+0.001)));
            
        elseif strcmpi(time_path, 'linear')
            
            if (t<=finish_year-start_year)  
            slope_aalpha = (aalpha_f - aalpha_0)/(finish_year-start_year);
            slope_eeta = (eeta_f - eeta_0)/(finish_year-start_year);
            slope_A      = (A_neutral_f - A_neutral_0)/(finish_year-start_year);
            
            aalpha_t(i,:)    = aalpha_0 + slope_aalpha*t;
            eeta_t(i,:)    = eeta_0 + slope_eeta*t;
            A_neutral_t(i,1) = A_neutral_0 + slope_A*t;
            
            elseif t>finish_year-start_year     
            aalpha_t(i,:)=aalpha_f;
            eeta_t(i,:)=eeta_f;
            A_neutral_t(i,1)= A_neutral_f;
            end
            
         elseif strcmpi(time_path, 'smoothed_linear')
            
            period = finish_year-start_year;
            
            slope_aalpha = (aalpha_f - aalpha_0)/(finish_year-start_year);
            slope_eeta = (eeta_f - eeta_0)/(finish_year-start_year);
            slope_A      = (A_neutral_f - A_neutral_0)/(finish_year-start_year);
            
            if (t<=finish_year-h-start_year) 
             
            aalpha_t(i,:)    = aalpha_0 + slope_aalpha*t;
            eeta_t(i,:)    = eeta_0 + slope_eeta*t;
            A_neutral_t(i,1) = A_neutral_0 + slope_A*t;
            
            elseif t>finish_year-h-start_year && t<finish_year+h-start_year 
                
            aalpha_t(i,:)    = aalpha_f + slope_aalpha*(-(t - h - period)^2/(4*h));
            eeta_t(i,:)    = eeta_f + slope_eeta*(-(t - h - period)^2/(4*h));
            A_neutral_t(i,1) = A_neutral_f + slope_A*(-(t - h - period)^2/(4*h));
            
            elseif t>finish_year+h-start_year     
            aalpha_t(i,:)=aalpha_f;
            eeta_t(i,:)=eeta_f;
            A_neutral_t(i,1)= A_neutral_f;
            end
           

        end
    end 
end 

%% Construct time-path for other aggregates
L_t      = repmat(ell_z,N,1);
Lwtd_t   = prod((psi_L.*L_t./(eeta_t.*(1-aalpha_t))).^(eeta_t.*(1-aalpha_t)),2);
aawtd_t  = diag(aalpha_t*eeta_t');

% - aggregates ------------------------------------------------------------
yt  = A_neutral_t.*((kt./aawtd_t).^aawtd_t).*Lwtd_t; 
kzt = kt.*(eeta_t.*aalpha_t./aawtd_t);
rt  = eeta_t(:,1).*aalpha_t(:,1).*yt./kzt(:,1) - ddelta;
yzt = ((aalpha_t.^(-aalpha_t)).*(1-aalpha_t).^(-(1-aalpha_t))).*(kzt.^aalpha_t).*((psi_L.*L_t).^(1-aalpha_t));
pzt = ((repmat(rt,1,Z)+ddelta)./aalpha_t).*kzt./yzt; 
wzt = (1-aalpha_t).*yzt.*pzt./L_t;
it  = (yt-ct)./yt;
kyt = kt./yt; 
mwt = sum(wzt.*ell_z,2);
lst = mwt./yt;
nkt = rt.*kt ./ (rt.*kt + mwt);

% - collect answers--------------------------------------------------------

paths.c=ct;         % consumptions
paths.k=kt;         % capital
paths.i=it;         % investment rate
paths.ktoy=kyt;     % capital-output ratio
paths.w_ave=mwt;    % average wage
paths.wz=wzt;       % vector of wages
paths.t=time;       % time
paths.ls=lst;       % labor share
paths.nk=nkt;       % net capital share
paths.r=rt;         % return to wealth
paths.y=yt;         % GDP per hour
paths.N=N;

% - initial ss values -----------------------------------------------------

ss0.c=css0;
ss0.i=iss0/yss0;
ss0.r=rss0;
ss0.wz=wzss0;
ss0.nk=nkss0;


end 

