function run()

clear all
close all
clc

%%% Parameters' block
bbeta=0.99;
phi_b=0.1;
phi_s=0.01;
psi_n=219;
psi_d_c=0.068;
psi_d_p=0.048;
eta=4;
theta=0.5;
gamma=0.08;
Gamma_d=0.005;
Gamma_l=0.001;
Gamma_x=0.115;
delta=0.025;
phi=0.1;
kappa=0.04;
aalpha=0.3;
mu=0.3;
tao=0.13;
GY=0.20;
B_c_Y=0.60;
B_p_Y=0.85;
BY_max_bar=0.92;
ssigma=0.015;
gamma_b=0.8354;
%%% Given SS--bbetae_bar,bbetaes_bar,S_star_c_bar,S_star_p_bar,x_bar,h_bar,h_s_bar,eps_bar,Z_bar,Zs_bar
bbetae_bar=bbeta;
bbetaes_bar=bbeta;
S_star_c_bar=0.5;
S_star_p_bar=0.5;
x_bar=0;
LK=1/3;
h_bar=0.2;
h_s_bar=0.2;
eps_bar=0.00054;
Z_bar=1;
Z_s_bar=1;
%%% Compute the staedy state parameters b_bar,b_s_bar, L_bar,L_s_bar,I_bar ,I_s_bar,div_bar ,div_s_bar,G_bar, G_s_bar, T_bar,B_bar, T_s_bar, B_s_bar, BY_max_bar,D_bar, D_s_bar
S_c_bar=1-S_star_c_bar;
S_p_bar=1-S_star_p_bar;
R_b_bar=1/bbeta;
R_b_s_bar=1/(bbeta*(1-eps_bar));
R_l_bar=1/bbeta;
R_l_s_bar=1/bbeta;
R_d_bar=(1-Gamma_d)/bbeta;
Q_bar=(1-theta)/theta; %since C_b=C_b_s
phit_bar=Q_bar^((aalpha)/(2*aalpha-1));
phits_bar=1/phit_bar;
R_d_s_bar=(Q_bar-Gamma_d)/(bbeta*Q_bar);
w_bar=psi_n*h_bar^(eta-1);
w_s_bar=psi_n*h_s_bar^(eta-1);
q_bar=1;
q_s_bar=1;
r_bar=1/bbetae_bar-1+delta;
r_s_bar=1/bbetaes_bar-1+delta;
Y_bar=(w_bar*h_bar)/((1-mu)*phit_bar);
Y_s_bar=(w_s_bar*h_s_bar)/((1-mu)*phits_bar);
K_bar=(Y_bar*mu*phit_bar)/r_bar;
K_s_bar=(Y_s_bar*mu*phits_bar)/r_s_bar;
L_bar=LK*K_bar;
L_s_bar=LK*K_s_bar;
I_bar=delta*K_bar;
I_s_bar=delta*K_s_bar;
div_bar=L_bar+K_bar*r_bar-I_bar-L_bar*R_l_bar;
div_s_bar=L_s_bar+K_s_bar*r_s_bar-I_s_bar-L_s_bar*R_l_s_bar;
Rs_bar=(1+S_c_bar*phi_s)/bbeta;
Rs_s_bar=(1+S_p_bar*phi_s)/bbeta;
rho_bar=div_bar/(Rs_bar-1);
rho_s_bar=div_s_bar/(Rs_s_bar-1);
prof_bar=Y_bar*phit_bar-h_bar*w_bar-K_bar*r_bar;
prof_s_bar=Y_s_bar*phits_bar-h_s_bar*w_s_bar-K_bar*r_bar;
G_bar=GY*Y_bar;
G_s_bar=GY*Y_s_bar;
B_bar=B_c_Y*4*Y_bar;
B_s_bar=B_p_Y*4*Y_s_bar;
b_bar=0.33*B_bar;
b_s_bar=0.33*B_s_bar;
T_bar=G_bar+R_b_bar*B_bar-B_bar;
H_star_bar=eps_bar*R_b_s_bar*B_s_bar;
H_h_s_bar=eps_bar*R_b_s_bar*b_s_bar;
T_s_bar=G_s_bar+H_star_bar+R_b_s_bar*(1-eps_bar)*B_s_bar-B_s_bar;
C_bar=(h_bar*w_bar-(psi_d_c*(psi_n*h_bar^eta)/eta)/(1-bbeta*R_d_bar)*R_d_bar+b_bar*R_b_bar+(rho_bar+div_bar)*S_c_bar+Q_bar*(rho_s_bar+div_s_bar)*S_p_bar+prof_bar+(psi_d_c*((psi_n*h_bar^eta)/eta))*(1/(1-bbeta*R_d_bar))-T_bar-rho_bar*S_c_bar-Q_bar*rho_s_bar*S_p_bar-b_bar-rho_bar*phi_s/2*S_c_bar^2-Q_bar*rho_s_bar*phi_s/2*S_p_bar^2)/(psi_d_c/(1-bbeta*R_d_bar)+1-psi_d_c/(1-bbeta*R_d_bar)*R_d_bar);
D_bar=(psi_d_c*(C_bar-(psi_n*h_bar^eta)/eta))/(1-bbeta*R_d_bar);
C_s_bar=(h_s_bar*w_s_bar-(psi_d_p*(psi_n*h_s_bar^eta)/eta)/(1-bbeta*R_d_s_bar)*R_d_s_bar+b_s_bar*R_b_s_bar*(1-eps_bar)+(rho_s_bar+div_s_bar)*S_star_p_bar+(rho_bar+div_bar)*Q_bar^(-1)*S_star_c_bar+prof_s_bar+H_h_s_bar+(psi_d_p*(psi_n*h_s_bar^eta)/eta)/(1-bbeta*R_d_s_bar)-T_s_bar-rho_s_bar*S_star_p_bar-rho_bar*Q_bar^(-1)*S_star_c_bar-b_s_bar-rho_s_bar*phi_s/2*S_star_p_bar^2-Q_bar^(-1)*rho_bar*phi_s/2*S_star_c_bar^2)/(psi_d_p/(1-bbeta*R_d_s_bar)+1-psi_d_p/(1-bbeta*R_d_s_bar)*R_d_s_bar);
D_s_bar=(psi_d_p*(C_s_bar-(psi_n*h_s_bar^eta)/eta))/(1-bbeta*R_d_s_bar);

%%% FSolve block

% Parameters' aliases
P = [S_star_c_bar, S_star_p_bar,phit_bar,x_bar,h_bar,h_s_bar,eps_bar,q_bar,q_s_bar,bbeta, psi_d_c, psi_d_p ,psi_n ,eta, phi_s, phi_b, b_bar,b_s_bar, L_bar, Gamma_l, Gamma_x ,gamma, L_s_bar, Gamma_d, theta, delta, phi, I_bar ,I_s_bar, kappa ,div_bar ,div_s_bar ,Z_bar, Z_s_bar,mu, aalpha, G_bar, G_s_bar, T_bar, tao, B_bar, T_s_bar, B_s_bar, BY_max_bar, gamma_b ,ssigma, D_bar, D_s_bar]; % put the names (aliases) of parameters here
% Load initial values of model's variables
% Set the number of endogenous variables here 
Xorig = x0input(69);
%% Here goes a non-linear solver
% Define options
options = optimset('Display', 'iter', 'Tolfun', 1.0e-26, 'MaxFunEvals', 5000,'MaxIter', 5000);
% Solve the system of non-linear equations
[Xout, Fval, exitflag] = fsolve(@calibration, Xorig, options, P);      
% Check whether the routine was executed successfuly
if exitflag <= 0
    disp('Error: not calculated');
end
% Check whether the solution is not complex
if isreal(Xout) == 0
    disp('Warning: complex steady state')
end
% Show the values of variables in the steady state
Xout;
% Create a cell array with values of vars in st. state
xCell = num2cell(Xout);
% Load steady state values of variables to aliases
[h,w,h_star,A,r, A_star, spread, spreadLD ,spreadLD_s, w_star, lampda, Y ,Y_s ,B ,T ,H_star, B_s, BY_max ,B_tilda_s , b_b, b_b_star, x, lampda_b, C_b, C_b_s , phit, phits , R_d ,lampda_star ,D_star ,R_d_s, C, D, Rs, Rs_s ,rho ,S_c ,Q ,rho_star ,S_p, b, R_b, R_b_s ,prof ,C_star ,T_star, b_s ,S_star_p ,S_star_c ,eps, prof_star, H_h_s, K, I, K_s, I_star, div, R_l, L, div_star , R_l_s, L_s ,r_star ,bbetae, bbetaes ,lampda_e ,lampda_es, q, q_star] = xCell{:}; % put the names (aliases) of variables here
%Clean up the environment
clear Xout xCell Pin Xorig t options exitflag F

%%% Values storage block

% Read the parameters' and variables' aliases
myvariables = whos;

for i=1:length(myvariables)
    str = sprintf('%30s;\n', myvariables(i).name);
    eval(str);
    fprintf(str);
end

for i=1:length(myvariables)
    str = sprintf('%30s_ss=%30s;\n', myvariables(i).name , myvariables(i).name);
    eval(str);
    fprintf(str)
end

save initial_value.mat;

dynare peregopaper noclearall
 
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Functions definitions block %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Fsolve function definition
function F = calibration(X,P) 
 
    % Convert arrays to cell arrays
    xCell = num2cell(X);
    pCell = num2cell(P);

    % Parameters' values assignment
    [S_star_c_bar, S_star_p_bar,phit_bar,x_bar,h_bar,h_s_bar,eps_bar,q_bar,q_s_bar, bbeta, psi_d_c, psi_d_p ,psi_n ,eta, phi_s, phi_b, b_bar,b_s_bar, L_bar, Gamma_l, Gamma_x ,gamma, L_s_bar, Gamma_d, theta, delta, phi, I_bar ,I_s_bar, kappa ,div_bar ,div_s_bar ,Z_bar, Z_s_bar,mu, aalpha, G_bar, G_s_bar, T_bar, tao, B_bar, T_s_bar, B_s_bar, BY_max_bar, gamma_b ,ssigma, D_bar, D_s_bar] = pCell{:}; 
    % Variables' values assignment
    [h,w,h_star,A,r, A_star, spread, spreadLD ,spreadLD_s, w_star, lampda, Y ,Y_s ,B ,T ,H_star, B_s, BY_max ,B_tilda_s , b_b, b_b_star, x, lampda_b, C_b, C_b_s , phit, phits , R_d ,lampda_star ,D_star ,R_d_s, C, D, Rs, Rs_s ,rho ,S_c ,Q ,rho_star ,S_p, b, R_b, R_b_s ,prof ,C_star ,T_star, b_s ,S_star_p ,S_star_c ,eps, prof_star, H_h_s, K, I, K_s, I_star, div, R_l, L, div_star , R_l_s, L_s ,r_star ,bbetae, bbetaes ,lampda_e ,lampda_es, q, q_star] = xCell{:}; 
    
    F = [
       -(S_star_c)+S_star_c_bar;
       -(S_star_p)+S_star_p_bar;
       -(S_c)+1-S_star_c;
       -(S_p)+1-S_star_p;
       -(Q)+((1-theta)/theta);
       -(phit)+Q^(aalpha/(2*aalpha-1));
       -(phits)+1/phit;
       -(x)+x_bar;
       -(h)+h_bar;
       -(h_star)+h_s_bar;
       -(eps)+eps_bar;
       -(R_b)+1/bbeta;
       -(R_b_s)+1/(bbeta*(1-eps));
       -(R_l)+1/bbeta;
       -(R_l_s)+1/bbeta;
       -(R_d)+(1-Gamma_d)/bbeta;
       -(R_d_s)+(Q-Gamma_d)/(bbeta*Q);
       -(w)+(psi_n*h^(eta-1));
       -(w_star)+(psi_n*h_star^(eta-1));
       -(q)+q_bar;
       -(q_star)+q_s_bar;
       -(bbetae)+bbeta*(S_c+S_star_c);
       -(bbetaes)+bbeta*(S_star_p+S_p);
       -(r)+1/bbetae-1+delta;
       -(r_star)+1/bbetaes-1+delta;
       -(Y)+(w*h)/((1-mu)*phit);
       -(Y_s)+(w_star*h_star)/((1-mu)*phits);
       -(K)+(Y*mu*phit)/r;
       -(K_s)+(Y_s*mu*phits)/r_star;
       -(L)+L_bar;
       -(L_s)+L_s_bar;
       -(I)+I_bar;
       -(I_star)+I_s_bar;
       -(div)+div_bar;
       -(div_star)+div_s_bar;
       -(Rs)+(1+S_c*phi_s)/bbeta;
       -(Rs_s)+(1+S_p*phi_s)/bbeta;
       -(rho)+div/(Rs-1);
       -(rho_star)+div_star/(Rs_s-1);
       -(prof)+Y*phit-h*w-K*r;
       -(prof_star)+Y_s*phits-h_star*w_star-K_s*r_star;
       -(B)+B_bar;
       -(B_s)+B_s_bar;
       -(b)+b_bar;
       -(b_s)+b_s_bar;
       -(T)+T_bar;
       -(T_star)+T_s_bar;
       -(BY_max)+BY_max_bar;
       -(B_tilda_s)+G_s_bar+B_s*R_b_s-T_star;
       -(H_star)+eps*R_b_s*B_s;
       -(H_h_s)+eps*R_b_s*b_s;
       -(D)+D_bar;
       -(D_star)+D_s_bar;
       -(C)+h*w+D*R_d+b*R_b+(rho+div)*S_c+Q*(rho_star+div_star)*S_p+prof-D-T-rho*S_c-Q*rho_star*S_p-b-rho*phi_s/2*S_c^2-Q*rho_star*phi_s/2*S_p^2;
       -(C_star)+h_star*w_star+D_star*R_d_s+b_s*R_b_s*(1-eps)+(rho_star+div_star)*S_star_p+(rho+div)*Q^(-1)*S_star_c+prof_star+H_h_s-D_star-T_star-rho_star*S_star_p-rho*Q^(-1)*S_star_c-b_s-rho_star*phi_s/2*S_star_p^2-Q^(-1)*rho*phi_s/2*S_star_c^2;
       -(b_b)+B-b;
       -(b_b_star)+B_s-b_s;
       -(lampda)+(C-(psi_n*h^eta)/eta)^(-1);
       -(lampda_star)+(C_star-(psi_n*h_star^eta)/eta)^(-1);
       -(lampda_e)+1;
       -(lampda_es)+1;
       -(A)+(phit*Y-Y*phit*aalpha-Y_s*Q*aalpha*phits)/(1-2*aalpha);
       -(A_star)+(Y*phit+Y_s*Q*phits-A)/Q;
       -(C_b)+A-Q*rho_star*phi_s/2*S_p^2-rho*phi_s/2*S_c^2-G_bar-I-C;
       -(C_b_s)+C_b;
       -(lampda_b)+theta*(C_b^theta*C_b_s^(1-theta))/C_b;
       -(spread)+R_b_s-R_b;
       -(spreadLD)+R_l-R_d;
       -(spreadLD_s)+R_l_s-R_d_s;
       
    ];  
end

%%% Initial value generator function
function x0 = x0input(len) 
    x0 = ones(len, 1)*0.2
end



