clc;
clear;
cd('/Users/mehreganameri/Dropbox/PhD/Macro & Inequality/2-agent/matlab')
 addpath /Applications/Dynare/4.4.3/matlab;
%addpath /Applications/Dynare/4.4_2013-10-23/matlab;

%% Set parameters
delta=0.025;        %   Capital depreciation rate
eta=3;              %   EoS between differentiated labour types
zeta=6;             %   EoS between differentiated goods
lambda=0.71;        %   ???????
theta_p=0.9;      %   Degree of price stickyness
theta_w=0.73;       %   Degree of wage stickyness
PHI=1.62;           %   ???
alpha=0.31;        %   Capitla share of production
beta=0.9;         %   Discount Factor
rho_R=0.62;         %   coeff of persistence, ??????????????
rho_a=0.95;         %   coeff of persistence, Technology shock
rho_r=0.42;         %   coeff of persistence, MP shock
rho_l=0.952;        %   coeff of persistence, labour supply shock
rho_c=0.882;        %   coeff of persistence, consumption preference shock
phi_pi=1.6;        %   MP rule response to inflation
phi_y=0.5;          %   MP rule response to output
sigma_a=1;      %   Sd. of technology shock
sigma_g=1;     %   Sd. of shock to government spending
sigma_r=1;       %   Sd. of shock to consumption preference
sigma_l=1;      %   Sd. of shock to consumption preference
sigma_c=1;      %   Sd. of shock to consumption preference
GY_ss = 0.2;       %   Steady state Goverment spending to GDP ratio
psi = 0.182;        %   Measure of transfer to the poor
k = 200;            %   ???????
upsilon = 1.25;        %   Labour suply elasticity
tau_w = 0;          %   Labour income tax rate
rho_b = -0.77;      %   Government spending rule response to Gov debt
rho_g = 0.8713;     %   Government spending shock coeff of persistence
tau_c = 0;      %   Consumption tax rate
tau_k = 0;      %   Capital tax rate
tau_p = 0;		%	Property tax rate 									
nu = 0.35;          %   Share of non-Ricardian HH

parameters = [delta eta zeta lambda theta_p theta_w PHI alpha beta rho_R phi_pi phi_y rho_a rho_r rho_l rho_c ...
                    sigma_a sigma_g sigma_r sigma_l sigma_c GY_ss psi k upsilon tau_w rho_b rho_g tau_c tau_k tau_p nu];
                  

%% Solving for steady state
X0=ones(21,1);
X=fsolve(@steady_stateHH,X0,[],parameters); 
X=fsolve(@steady_stateHH,X0,[],parameters); 
X=fsolve(@steady_stateHH,X0,[],parameters); 



%% Extracting SS values

Cr_ss = X(1); 
Cn_ss = X(2);
C_ss = X(3);
Nr_ss = X(4);
Nn_ss = X(5);
L_ss = X(6);
W_ss = X(7);
R_ss = X(8);
Rk_ss = X(9);
Lambda_ss = X(10);
K_ss = X(11);
MC_ss = X(12);
Gamma_1_ss = X(13);
Gamma_2_ss = X(14);
Y_ss = X(15);
Pstar_P_ss = X(16);
Pi_ss = X(17);
sp_ss = X(18);
G_ss = X(19);
B_P_ss = X(20);
Pr_ss = X(21);



save parameters delta eta zeta lambda theta_p theta_w PHI alpha beta rho_R phi_pi phi_y rho_a rho_r rho_l rho_c ...
                    sigma_a sigma_g sigma_r sigma_l sigma_c GY_ss psi k upsilon tau_w rho_b rho_g tau_c tau_k tau_p nu ...
                    Cr_ss Cn_ss C_ss Nr_ss Nn_ss L_ss W_ss R_ss Rk_ss Lambda_ss K_ss MC_ss Gamma_1_ss Gamma_2_ss Y_ss... 
                    Pstar_P_ss Pi_ss sp_ss G_ss B_P_ss Pr_ss;


%% after having run the mod-file:
%% first order:
%set_param_value('delta',d);
dynare nk3HH;
%%
%plot(sp_eta_R);
%print('-dpdf','sp_eta_R');
%plot(sw_eta_R);
%print('-dpdf','sw_eta_R');

%% second order with pruning:

%dynare newkeynesian_pruning;
%%
%plot(sp_eta_R);
%print('-dpdf','sp_eta_R_pruning');
%plot(sw_eta_R);
%print('-dpdf','sw_eta_R_pruning');

varD_all = oo_.gamma_y{7};
var_all = oo_.var;
mean_all = oo_.mean;

save('varD_all', 'varD_all');
save('var_all', 'var_all');
save('mean_all', 'mean_all');