%% stephane [DOT] adjemian [AT] ens [DOT] fr

function [ys, check] = steadystate(ys, exe)
global M_

%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%

NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                   %    Get the name of parameter i. 
  eval([ paramname ' = M_.params(' int2str(i) ');']);         %    Get the value of parameter i.
end                                                           % End of the loop.  

NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
for i = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
  eval([ varname ' = 0;']);                					  %    Start with 0 as default value.
end                                                           % End of the loop.

check = 0;

%%
%% END OF THE FIRST MODEL INDEPENDENT BLOCK.

%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.

% Demography & external data
J  = 80;

% Loading demography - to be modified so that it is taken from the mod file
mort_dt	= xlsread('Demography_EA.xlsx', 'Dynare_mort');
n_1_dt	= xlsread('Demography_EA.xlsx', 'Dynare_n_y');
ret_dt  = xlsread('Demography_EA.xlsx', 'Dynare_JR');

mort = mort_dt(dt,:);
n_1  = n_1_dt(dt);
ret  = ret_dt(dt,:);

N_rel_all = zeros(J,1);
N_rel_all(1) = 1;
for j = 2:J
    N_rel_all(j) = (1-mort(j-1)) * N_rel_all(j-1)/(1+n_1);
end
N_rel = sum(N_rel_all);
n = n_1;

%% Solving for steady state

global start options
start.r   = 0.01;
start.beq = 0.01;

options	= optimset('MaxIter', 10000, 'MaxFunEvals', 1000, 'TolFun', 1e-08);

fun     = @(x)solve_rbeqpen2(x, alffa, betta, delta, ret, J, mort, n, N_rel_all, N_rel);
[x,~]	= fsolve(fun, [start.r, start.beq], options);

r   = x(1);
beq = x(2);

w	= (1-alffa) * (alffa/(r+delta))^(alffa/(1-alffa));

%% Solve the problem of OLG households given prices, taxes and transfers
% Calculate the sum of lifetime discount factors
surv = cumprod(1-mort);
disc = 1;
for j = 2:J
    disc = disc + betta^(j-1) * surv(j-1);
end

% Calculate the NPV of lifetime income
lti = 0;
for j = 2:J
    lti = lti + ret*((w + beq)/(1+r)^(j-1));
end


% Find c1 that satisfies the lifetime budget constraint
c1 = lti/disc;

a_0 = 0;

c_all = zeros(J, 1);
h_all = zeros(J, 1);
a_all = zeros(J, 1);

c_all(1) = c1;
h_all(1) = 1;
a_all(1) = w + beq + (1+r)*a_0 - c_all(1);

for j = 2:J-1
    c_all(j) = betta*(1-mort(j-1))*c_all(j-1)*(1+r);
    h_all(j) = ret;
    a_all(j) = w*ret + beq + (1+r)*a_all(j-1) - c_all(j);
end



c = sum(c_all.*N_rel_all)/N_rel;
h = sum(h_all.*N_rel_all)/N_rel;
a = sum(a_all.*N_rel_all)/N_rel/(1+n);

% Creating output
for j = 1:J
    eval(['c_' int2str(j) ' = ' num2str(c_all(j),16) ';']);
    eval(['h_' int2str(j) ' = ' num2str(h_all(j),16) ';']);
    eval(['a_' int2str(j) ' = ' num2str(a_all(j),16) ';']);
    eval(['N_rel_' int2str(j) ' = ' num2str(N_rel_all(j),16) ';']);
end

%% Remaining variables
k	= (alffa/(r+delta))^(1/(1-alffa))*h;
inv	= (delta+n)*k;
y	= k^alffa*h^(1-alffa);
feas = 0;
    
%%
%% END OF THE MODEL SPECIFIC BLOCK.

%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state values of the endogenous variables of
%% the model.
%%
NumberOfEndogenousVariables = M_.endo_nbr;                    % Number of endogenous variables.
ys = zeros(NumberOfEndogenousVariables,1);                    % Initialization of ys (steady state).
for i = 1:NumberOfEndogenousVariables                         % Loop...
  varname = deblank(M_.endo_names(i,:));                      %    Get the name of endogenous variable i.                     
  eval(['ys(' int2str(i) ') = ' varname ';']);                %    Get the steady state value of this variable.
end                                                           % End of the loop.
%%
%% END OF THE SECOND MODEL INDEPENDENT BLOCK.
