%% Script to compute steady state of the model:

function [ys,check] = oilpol_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs:
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output:
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
    paramname = deblank(M_.param_names(ii,:));
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%% Enter model equations here

% Get the element PIcF in ys
ramsey_flag=M_.endo_nbr>=2*M_.orig_endo_nbr-1;
if ramsey_flag
    PIcF_loc=strmatch('PIcF',M_.endo_names);
    PIcF=ys(PIcF_loc);
else
    PIcF = 1.000;
end

%tic

% CALIBRATION:

Numsec = 2;                                 % # of industries

betta = 0.99;                            % Time discount factor
sigma = 1;                             % Inverse of IES
varphi = 2;                              % Inverse of Frisch
phi = [0.60,0.40];                     % CD-production materials Home
phiF = [0.60,0.40];                    % CD-production materials Foreign
psi = 1 - phi;                              % CD-production labor Home
psiF = 1 - phiF;                            % CD-production labor Foreign
xi = [0.35,0.65];                      % SS Consumption share (of total consumption) Home
xiF = [0.35,0.65];                     % SS Consumption share (of total consumption) Foreign
zeta =  [0.70, 0.30;                  % SS Material share (of total material) Home (Note: each column sums to one)
    0.30, 0.70];
zetaF = [0.70, 0.30;                  % SS Material share (of total material) Foreign (Note: each column sums to one)
    0.30, 0.70];
gammaex = [0.60, 0.15];                 % SS export-GDP ratio
gammaim = [0.60, 0.15];                 % SS import-GDP ratio
L = 1/3;                                    % Normalization of hours Home
eps_p = 0.2;                          % Markup goods markets
eps_w = 0.2;                          % Markup labor markets
theta_pj = [0.25, 0.75];

% Oil sector:
gammao = 0.2;    % Oil share of total GDP
alphao = 0.2;      % Non-land share of oil output
phio = 0.7;%.16;  % Use of mainland intermediates
psio = 1 - phio;    % Use of mainland labor
Oilshare_ss = gammao/(1 - gammao);  % Oil share of mainland GDP.
zetao = [0.4, 0.6]; % Oil input shares from goods and services.

if length(zeta) ~= Numsec
    error('Error: Input-output matrix must be square and of size equal to # sectors!');
end

%% Recursive steady state solution:

% REST OF THE WORLD:

RF = PIcF/betta;
BarPrjF = ((1 - theta_pj)./(1 - theta_pj*PIcF^(1/eps_p))).^eps_p;
VjF = (1 - theta_pj).*BarPrjF.^(-(1 + eps_p)/eps_p)./(1 - theta_pj*PIcF^((1 + eps_p)/eps_p));
RMCjF = (BarPrjF/(1 + eps_p)) .* (1 - betta*theta_pj*PIcF^((1 + eps_p)/eps_p)) ./ (1 - betta*theta_pj*PIcF^(1/eps_p));
% Normalization
CF = 1;
% C_j = xi_j*C
CjF = xiF*CF;
% Y_j = C_j + M_{j1} + ... + M_{jJ}
iomatrixF = [zetaF(1,1)*phiF(1), zetaF(1,2)*phiF(2);
    zetaF(2,1)*phiF(1), zetaF(2,2)*phiF(2)];
XjF = ((eye(Numsec) - iomatrixF)\CjF')';
YjF = XjF.*VjF;
% 1 + eps_p = 1 + Phi_j/Y_j
PhijF = (1./RMCjF - VjF).*XjF;%eps_p*YjF;
% M_j/Y_j = phi_j
MjF = phiF.*XjF;%phiF.*YjF;
MjjF = bsxfun(@times, zetaF, MjF);
% GDP = Omega*L
OmegaF = CF/L;
% Omega = (1 + eps_w)*MRS
chi_NF = OmegaF/((L^varphi)*(CF^sigma)*(1 + eps_w));
% RMC_j = Z_j^(-1) * (1/phi_j)^(phi_j) * (Omega/psi_j)^(psi_j)
ZjF = (1./RMCjF) .* (1./phiF).^(phiF) .* (OmegaF./psiF).^(psiF);
% Omega*N_j/Y_j = psi_j
NjF = psiF.*XjF/OmegaF;%psiF.*YjF/OmegaF;
% L = N_j/mu_j
mujF = NjF/L;
% Value added
GDPjF = (1 - phiF).*XjF;%YjF - MjF;
GDPF = sum(GDPjF);
LambdaF = CF^(-sigma);
GF = LambdaF*XjF.*RMCjF./(1 - betta*theta_pj*PIcF^((1 + eps_p)/eps_p));
HF = LambdaF*XjF.*BarPrjF./(1 - betta*theta_pj*PIcF^(1/eps_p));
TildeGF = PIcF^((1 + eps_p)/eps_p)*GF;
TildeHF = PIcF^(1/eps_p)*HF;

%gdpF = log(GDPF);
VAF = CF;
%cF = log(CF);
%lambdaF = log(LambdaF);
LambdaauxF = LambdaF/betta;
%PIcF = PIcF;
RRF = 1/betta;
%rF = log(RF);
NF = L;
RWF = OmegaF;
UtilityF = log(CF) - chi_NF*(L)^(1 + varphi)/(1 + varphi);
WelfareF = UtilityF/(1 - betta);
for i=1:Numsec
    % Variables
    evalc(sprintf('C%dF = CjF(i)', i));
    evalc(sprintf('GDP%dF = GDPjF(i)', i));
    evalc(sprintf('Y%dF = YjF(i)', i));
    evalc(sprintf('X%dF = XjF(i)', i));
    evalc(sprintf('M1%dF = MjjF(1,i)', i));
    evalc(sprintf('M2%dF = MjjF(2,i)', i));
    evalc(sprintf('M%dF = MjF(i)', i));
    evalc(sprintf('N%dF = mujF(i)*L', i));
    evalc(sprintf('L%dF = L', i));
    evalc(sprintf('PIw%dF = PIcF', i));
    evalc(sprintf('RW%dF = OmegaF', i));
    evalc(sprintf('PI%dF = PIcF', i));
    evalc(sprintf('Pr%dF = 1', i));
    evalc(sprintf('barPr%dF = BarPrjF(i)', i));
    evalc(sprintf('G%dF = GF(i)', i));
    evalc(sprintf('tildeG%dF = TildeGF(i)', i));
    evalc(sprintf('H%dF = HF(i)', i));
    evalc(sprintf('tildeH%dF = TildeHF(i)', i));
    evalc(sprintf('V%dF = VjF(i)', i));
    evalc(sprintf('Prm%dF = 1', i));
    evalc(sprintf('RMC%dF = RMCjF(i)', i));
    evalc(sprintf('Z%dF = ZjF(i)', i));
    % Parameters
    evalc(sprintf('theta_p%d = theta_pj(i)', i));
    evalc(sprintf('Phi%dF = PhijF(i)', i));
    evalc(sprintf('phi%dF = phiF(i)', i));
    evalc(sprintf('psi%dF = 1 - phiF(i)', i));
    evalc(sprintf('xi%dF = xiF(i)', i));
    evalc(sprintf('mu%dF = mujF(i)', i));
    evalc(sprintf('zeta1%dF = zetaF(1,i)', i));
    evalc(sprintf('zeta2%dF = zetaF(2,i)', i));
end



% CHECK SOLUTION:

if abs(sum(mujF) - 1) >= 1e-10
    warning('Warning: Employment shares in labor markets do not add up to total labor force');
end
if any(sum(zetaF,1)) ~= 1
    warning('Warning: Material shares do not sum to 1');
end

%% end own model equations

for iter = 1:length(M_.params) %update parameters set in the file
    eval([ 'M_.params(' num2str(iter) ') = ' M_.param_names(iter,:) ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
    varname = deblank(M_.endo_names(ii,:));
    eval(['ys(' int2str(ii) ') = ' varname ';']);
end

%keyboard

% if ramsey_flag && any(isnan(ys))
%     ys(108)=ys(85)*(-(exp(exo(1))*1/(ys(32))*getPowerDeriv(ys(32)/(ys(32)),M_.params(10),1)));
%     ys(109)=ys(107)*(-(exp(exo(2))*1/(ys(54))*getPowerDeriv(ys(54)/(ys(54)),M_.params(10),1)));
% end
% if any(isnan(ys))
%     check=1;
% end