% Solving the Competitive equilibirum 
% ==================================== 

var psi     $\psi$             
    varho   $\varrho$            
    theta   $\theta$              
    v       $\nu$                 
    Rap     $Rap$                   
    R       $R$                  
    Rk      $Rk$                  
    Rh      $Rh$                             
    rk      $rk$                
    rh      $rh$                 
    exr     $expi$               
    u       $u$                 
    kh      $kh$                 
    gr      $gr$                  
    A       $A$                   
    sigmat  $\sigma$;            

varexo e    $e$                  
       eta  $\eta$;              

predetermined_variables kh theta;

parameters alpha    $\alpha$             
           gamma    $\gamma$           
           delta    $\delta$            
           rho      $\rho$              
           z        $z$                
           sigma    $\sigma$              
           eps      $\epsilon$          
           beta     $\beta$            
           Abar     $A$;                

gamma=3;
delta=0.06;
eps=0.5; 
beta=0.96;
alpha=0.36;
sigma=0.107544;
Abar=0.302734;
rho=0.95;
z=0.02;

model;

% Asset markets 

psi=v^(1-gamma);                                                                                    
[dynamic] varho = (psi*Rap)^(1/(1-gamma));                                       
[static]  varho = Rap^(1/(1-gamma));
[dynamic] Rap(+1)=Rk(+1)^(1-gamma)+(1-gamma)*Rk(+1)^(-gamma)*theta(+1)*exr(+1)*(1-0.5*gamma*Rk(+1)^(-1)*theta(+1)*exr(+1)) ; 
[static]  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);
[dynamic] exr=Rk-Rh+exp(sigmat)*eta;                                     
[static]  exr=0;
[dynamic] theta=(1/gamma)*((psi*Rk^(-gamma)*exr)/(psi*Rk^(-gamma-1)*(exr)^2));                    
[static]  kh=0.640833;

% theta=(Rk/gamma)*((Rh-Rk)/((Rh-Rk)^2 + sigma^2));
%optimal choices 

[dynamic] v=(1+(beta^eps)*(varho(+1))^(eps-1))^(1/(eps-1)); 
[static]  v=u^(1/(1-eps));
[dynamic] u=v^(1-eps);
[static]  u=1-(beta^eps)*varho^(eps-1);
theta=1/(1+kh);

%supply side

Rk=1+rk-delta;
Rh=1+rh-delta;
R=Rk+theta*(Rh-Rk);
rk=exp(A)*alpha*kh^(alpha-1);
rh=exp(A)*(1-alpha)*kh^(alpha);


%growth rate
[dynamic] gr=R*(1-u(-1))-1;
[static] gr=R*(1-u)-1;

%Exogenous Processes
A=(1-rho)*Abar+ rho*A(-1)+e;
[dynamic] sigmat=sigma*(1-z*ln(A))^2;
[static] sigmat=sigma;
end;

steady_state_model;
 kh=0.640833;      % ==> The value was computed once I solved the system seperately 
 theta=1/(1+kh);
 A=Abar;
 rk=A*alpha*kh^(alpha-1);
 rh=A*(1-alpha)*kh^(alpha);
 Rk=1+rk-delta;
 Rh=1+rh-delta;
 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);
 varho=Rap^(1/(1-gamma));
 u=1-(beta^eps)*varho^(eps-1);
 v=u^(1/(1-eps));
 psi=v^(1-gamma);
 R=Rk+theta*(Rh-Rk);
 gr=R*(1-u)-1;
 sigmat=sigma;
 exr=0;
end;

steady(solve_algo=3,maxit=1000000);
check;

shocks;
var e; stderr 1;
end;
write_latex_dynamic_model;
write_latex_static_model;

stoch_simul(order=3,periods=100000,drop=1000,irf=0);