% This is the Plain Vanilla RBC.


% The model implements exogenous technological progress. 

% Copyright: Peter A. Schmid, Catholic University of Eichstaett-Ingolstadt
% December 14th, 2009.


clc;
close all; 


var uc, ul, mpk, mpl, y, c, inv, k, l, re, rf, ep, z;
varexo e;

parameters a beta, eta, psi, delta, alpha, rho, stde;

a        = 1.000;
beta     = 0.99;
eta      = 2.5;
psi      = 2.0;

delta    = 0.025;
alpha    = 0.36;
rho      = 0.95;

stde     = 0.01;

% Help variables for steady state determination!
ykbar    = (1/(beta*a^(-eta))-1+delta)/(alpha*a);
ikbar    = (a+delta-1)/a;
ckbar    = ykbar-ikbar;

helpbar  = (1-alpha)*ykbar/ckbar;
lbar     = helpbar/(psi+helpbar);
kbar     = ykbar^(1/(alpha-1))*lbar*a^(alpha/(alpha-1));
cbar     = ckbar*kbar;
ibar     = ikbar*kbar;
ybar     = kbar^alpha*lbar^(1-alpha)*a^(-alpha);

ucbar    = cbar^(-eta)*(1-lbar)^(psi*(1-eta));
ulbar    = -psi*cbar^(1-eta)*(1-lbar)^(psi*(1-eta)-1);
mpkbar   = alpha * kbar^(alpha-1)*lbar^(1-alpha)*a^(1-alpha);
mplbar   = (1-alpha)*kbar^alpha*lbar^(-alpha)*a^(-alpha);


model;
% Marginal utilities, marginal products and definitions
uc     = c^(-eta)*(1-l)^(psi*(1-eta));
ul     = -psi*c^(1-eta)*(1-l)^(psi*(1-eta)-1);
mpk    = exp(z) * alpha * a^(1-alpha)* k(-1)^(alpha-1) * l^(1-alpha);
mpl    = exp(z) * (1-alpha) * a^(-alpha)* k(-1)^alpha * l^(-alpha); 
y      = exp(z) * a^(-alpha)*k(-1)^alpha * l^(1-alpha);

% First Order Conditions
-ul = uc * mpl;
uc  = beta*a^(-eta) * uc(+1) * rf;
uc  = beta*a^(-eta) * uc(+1) * re(+1);
re = mpk - delta + 1;

% Equations of Motion
a*k  = (1-delta) * k(-1) + a*inv;

z = rho * z(-1) + stde * e;

% Product Market Equilibrium
inv      = y - c;

% Definition
ep  = re - rf(-1);

end;

initval;
k = kbar;
c = cbar;
l = lbar;

z = 0.0;

uc  = ucbar;
ul  = ulbar;
mpk = mpkbar;
mpl = mplbar;

rf = 1/(beta*a^(-eta));

y = ybar;

inv = ibar;
re = mpk -delta+1;
ep = re - rf;
e = 0.0;
end; 

steady(solve_algo=0);

shocks;
var e = 1;
end;
 
stoch_simul(order=2,hp_filter=1600,irf=0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some Results: Standard Deviations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

statistic1 = 100*sqrt(diag(oo_.var(5:11,5:11)))./oo_.mean(5:11);
table('Relative standard deviations in %',strvcat('VARIABLE','REL. S.D.'),M_.endo_names(5:11,:),statistic1,10,8,4);