var u gr theta rk rh Rk Rh kh R varho Rap v;
parameters gamma eps beta delta sigma alpha A;

gamma=3;
delta=0.06;
eps=0.5;
beta=0.96;
alpha=0.36;
sigma=0.107544;
A=0.302734;

%sigmay=sqrt(0.0313);
	 


model;

%%Equilibrium 

theta=(Rk/gamma)*((Rh-Rk)/((Rh-Rk)^2 + sigma^2));
theta=1/(1+kh);
Rk=1+rk-delta;
Rh=1+rh-delta;
rk=A*alpha*kh^(alpha-1);
rh=A*(1-alpha)*kh^(alpha);
%sigma=sigmay/theta;


%% Extended model 

u=1-(beta^eps)*varho^(eps-1);
varho=Rap^(1/(1-gamma));
Rap=Rk^(1-gamma)+ (1-gamma)*(Rk^(-gamma))*theta*(Rh-Rk)-0.5*(Rk^(-gamma-1))*gamma*(1-gamma)*theta^2*((Rh-Rk)^2 + sigma^2);
u=v^(1-eps);
R=Rk+theta*(Rh-Rk);
gr=R*(1-u)-1;
end;

initval;
kh=2;
theta=1/(1+kh);
Rk=1+rk-delta;
Rh=1+rh-delta;
rk=A*alpha*kh^(alpha-1);
rh=A*(1-alpha)*kh^(alpha);

%sigma=sigmay/theta;
R=Rk+theta*(Rh-Rk);
varho=1;
Rap=1;
u=1-(beta^eps)*varho^(eps-1);
Rap=Rk^(1-gamma)+ (1-gamma)*(Rk^(-gamma))*theta*(Rh-Rk)-0.5*(Rk^(-gamma-1))*gamma*(1-gamma)*theta^2*((Rh-Rk)^2 + sigma^2);
v=1;
gr=R*(1-u)-1;
end; 


steady(solve_algo=3,maxit=100000);