eta        =1.013;
delta_p    =0.028;
betta      =0.9817; 
theta        =5.83;
psi_l      =0.4285;
eta        =1.013;
betta      =0.9817;  
delta_p    =0.028; 
psi_m      =0.4625; 
alpha_kp   =0.42; 
alpha_kgp  =0.23;
alpha_lp   =0.5; 
alpha_kg   =0.65; 
alpha_lg   =0.5; 
delta_g    =0.03; 
phi        =40; 

phi_ig     =0.2;
alpha_nd   =0.05; 
alpha_f    =0.14; 
r          =0.15; 
rd         =0.1; 
phi_f      =0.2;
theta      =5.83;
phi_p       =3.15;


% the steady state computation follows 
pi =1.15;
A =1;
Y_oil =1;
G_c =1;
mu =eta*pi;
i =((pi*eta)/betta)-1;
r_kp =((1+i)/pi)+ delta_p -1;
mc =(1/theta)*((1-betta*eta^2)* phi_p *(pi-1)* pi + theta-1);

options = optimset('Display','iter','TolFun',1e-16); % set options for numerical solver
%options = optimoptions('fmincon','Display','Final');

[W,fval]=fzero(@(W) (1-(eta+delta_p-1)*(alpha_kp*mc/r_kp)-(phi_p/2)*(pi-1)^2)*(alpha_lp/W)^((alpha_lp-alpha_kp)/alpha_kgp)...
+((alpha_lg/W)^(alpha_lg/alpha_kg)-(eta+delta_p-1)*(1+phi*(delta_g/2)))*((W/(alpha_lp*mc))^(alpha_kp/alpha_kgp)*(r_kp/(alpha_kp*mc))^(alpha_kp/alpha_kgp)/mc^((alpha_lp-alpha_kp)/alpha_kgp))-(W/psi_l) ,0.2,options);

C = W / psi_l;
m =((1+i)/i) * C* psi_m;
K_g = (W/(alpha_lp*mc))^(alpha_kp/alpha_kgp)*(r_kp/(alpha_kp*mc))^(alpha_kp/alpha_kgp)/mc^((alpha_lp-alpha_kp)/alpha_kgp);
Y_g =  (alpha_lg/W)^(alpha_lg/alpha_kg)* K_g;
Y_p = (alpha_lp*mc/W)^(alpha_lp/alpha_kgp)*(alpha_kp*mc/r_kp)^(alpha_kp/alpha_kgp)* K_g;
K_p = ( alpha_kp * mc / r_kp )* Y_p;
I_ge = (eta + delta_g - 1)* K_g;
I_g = I_ge * ( 1 + phi * ( delta_g / 2 ));



