clear
close all
clc

global A_t A_n beta_n A_m alpha eta0 gamma kg_bar delta_g delta_n e n g ca_bar theta_g 


 A_t=1;
%A_t=2;

%population
n=1;

%modern agriculture
A_n=1.2;
A_n=2;
%A_n=2.1


beta_n=.6;

% manufacturing
A_m=1.2;
A_m=2;
%A_m=2.1;

%preference parameter
alpha=0.2;

%transportation costs function

eta0=1.2;

gamma=0.8;
kg_bar=1;
delta_g=0.03;

%education

delta_n=0.06;

e=2;

%Government spending
g=0.1;

%subsistence food
ca_bar=0.6;


 

ne_o=1;
eta_o=1;
p_o=1;
csM_o=0.5;
csA_o=0.5;
G_o=1;
theta_o=0.2;
z_o=1;
cmA_o=0.5;
cmM_o=0.5;


x0=[ne_o, eta_o, p_o, csM_o, csA_o, G_o, theta_o, z_o, cmM_o cmA_o];          

 
%[x,fval]  =   fsolve(@fun_no_land,x0,optimset('Display','off','MaxIter',1000,'TolFun',1e-10,'MaxFunEvals',1000));


% Solve using fmincon
       
j=0;
max=-100;
%for theta_g=0.77:0.001:0.82
for theta_g=0.76:0.0001:0.78
     j=j+1;

LB=0.01*ones(1, 10);  
%UB=5*ones(1,14);  
UB=[1, 1, 10, 10, 10, 10, 10, 10, 10, 1, 10, 10, 10];

% puts a lower constraint of zero (but not exactly zero to avoid divisions by zero and NaN problems)  and upper limit of 5 on all parameters.
options = optimset('MaxFunEvals',100000,'MaxIter',100000);
x =  fmincon(@fun_no_land_modif, x0,[],[],[],[],LB,UB,[],options);
x = x;
y=fun_no_land(x);
sum_of_abs_of_y_coordinates(j) = fun_no_land_modif(x);
csM(j)=x(4);
csA(j)=x(5);
cmM(j)=x(9);
cmA(j)=x(10);
mat1(j, :)=x;
wage_manuf(j)=x(3)*A_m;
plan_obj(j, :)= [(n-x(1))*(alpha*log(x(5)-ca_bar)+log(x(4)))+x(1)*(alpha*log(x(10)-ca_bar)+log(x(9))), theta_g];


%find maximum
if (max<plan_obj(j,1))&& (sum_of_abs_of_y_coordinates(j)<0.00001)
    max=plan_obj(j,1); theta_g_max=theta_g;
end
%theta_g,
%pause,
end;

 mat1 = horzcat(mat1, wage_manuf');

 % use Broyden Method to Solve the System of Equations
 %j=0;
 
theta_g=theta_g_max;
% 
 x0=[ne_o eta_o p_o csM_o csA_o G_o theta_o z_o cmM_o cmA_o]';          
 
  xx = broyden('fun_no_land_BR',x0);
  
  [y, xss] = fun_no_land_BR(xx);

 mat_xx=real(xx);
 
 
 
%% Define unknown variables
 
% ne_an=xx(1);
% eta_an=xx(2);
% p_an=xx(3);
% csM_an=xx(4);
% csA_an=xx(5);
% G=xx(6);
% theta_an=xx(7);
% z_an=xx(8);
% cmM_an=xx(9);
% cmA_an=xx(10);
% 
% YM=A_m*(1-theta_an)*ne_an;
% YtA=A_t*(n-ne_an);
% YnA=theta_an*ne_an*p_an*A_m/beta_n;

%% Evaluate steady state and stack values in Dynare vector 


 ne_init=xx(1);
eta_init=xx(2);
p_init=xx(3);
csM_init=xx(4);
csA_init=xx(5);
G_init=xx(6);
theta_init=xx(7);
z_init=xx(8);
cmM_init=xx(9);
cmA_init=xx(10);

theta_g_init=theta_g;

  save steady_state_values.mat ne_init eta_init p_init csM_init csA_init G_init theta_init z_init cmM_init cmA_init theta_g_init
 
 dynare StructTransf_no_land_D.mod noclearall