% Basic Dynare code for grad macro

% Eric Sims
% University of Notre Dame
% Spring 2011

var y i k a c;
varexo e;

parameters alpha beta delta rho sigma sigmae;

alpha = 0.33;
beta = 0.99;
delta = 0.025;
rho = 0.95;
sigma = 1;
sigmae = 0.01;

model;

exp(c)^(-sigma) = beta*(exp(c(+1))^(-sigma)*(alpha*exp(a(+1))*exp(k)^(alpha-1)+(1-delta)));
exp(y) = exp(a)*exp(k(-1))^(alpha);
exp(k) = exp(i) + (1-delta)*exp(k(-1));
exp(y) = exp(c) + exp(i);
a = rho*a(-1)+e;

end;

initval;

k = log(29);
y = log(3);
a = 0;
c = log(2.5);
i = log(1.5);

end;
shocks;
var e = sigmae^2;
end;

steady;

set_dynare_seed(92);
stoch_simul(hp_filter=1600,order=1,irf=200);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%SIMULATION%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Let s_t be a (m x 1) vector of states
% Let x_t be a (n x 1) vector of controls (including "static" controls like output)
% Here, m=2, n=3
% State space representation of the system:
%	s_t = A s_{t-1} + B epsilon_t
%	x_t = Phi s_t
% Where A is (m x m), epsilon_t is (w x 1), B is (m x w);
% Phi is (n x m) and is the policy function
% Dynare writes the system by substituting to get
%	x_t = Phi A s_{t-1} + Phi B epsilon_t
% Thus we can write the full system as
%	s_t = A s_{t-1} + B epsilon_t
%	x_t = C s_{t-1} + D epsilon_t
% Denoting Y_t = [s_t x_t]' and Psi = [A C]' and Omega = [B D]', finally we have
% Y_t = Psi s_{t-1} + Omega epsilon_t




m = 2; 				% number of states
n = 3;				% number of controls

%Using notation above for our system in matrix form:
psi = oo_.dr.ghx;		% psi
omega = oo_.dr.ghu;		% Omega
ss = oo_.dr.ys;			% steady state values

T = 200;			% Number of periods to simulate

es = sigmae*randn(T,1);
Xsim = zeros(n+m,T);

Xsim(:,1) = omega*es(1,1);	% 

for j = 2:T
	Xsim(:,j) = psi*Xsim(3:4,j-1) + omega*es(j,1);
end

% add back in steady state values so that these are expressed as log
% levels, not log deviations

for j = 1:T
	XSim(:,j) = Xsim(:,j) + ss;
end

%%%%%%PLOT SIMULATION%%%%%%%
figure
subplot(2,3,1)
plot(Xsim(1,:))
title('Output?')
subplot(2,3,2)
plot(Xsim(2,:))
title('Interest rate?')
subplot(2,3,3)
plot(Xsim(3,:))
title('Kapital?')
subplot(2,3,4)
plot(Xsim(4,:))
title('Productivity?')
subplot(2,3,5)
plot(Xsim(5,:))
title('Consumption?')