% ===================================================
% This code solves the basic RCK model
% (...)
% ====================================================

% ====================================================
% PREFACE: AS IN MATLAB
	close all;
	clc;
	format short;
	delete('RCK_Ex_1.txt');
	diary on;
	diary('RCK_Ex_1.txt');
% ====================================================
	
% =======================================================	
% FROM HERE ONWARDS, DYNARE CODE
% =======================================================

% Initial Steady State
display('Inital Steady State');
A = 1;
beta = 0.95;
delta = 0.02;
alpha = 0.25;
n = 0.0;
kss_0 = ((alpha*A)/((1/beta) - (1 - delta)))^(1/(1 - alpha))
css_0 = A*(kss_0^(alpha)) - kss_0*(n + delta)


% -------------------------------------------------------
% DECLARE ENDOGENOUS VARIABLES 
% WITHOUT TIME, AS IN THE STEADY STATE
% -------------------------------------------------------
% Number of variables: 5
    var  k,  c ; 

% -------------------------------------------------------
% LIST OF PARAMETERS 
% -------------------------------------------------------
  parameters beta, sigma, delta, alpha, n, A;
	beta = 0.95;
	sigma = 1.5;	
	delta = 0.02;
	alpha = 0.25;
	n = 0.0;
	A = 1.5;

% Final Steady State	
kss = ((alpha*A)/((1/beta) - (1 - delta)))^(1/(1 - alpha));
css = A*(kss^(alpha)) - kss*(n + delta);
 

	
% -------------------------------------------------------
% MODEL DESCRIPTION
% -------------------------------------------------------
model;  

	% i) Resource constraint
		A*(k(-1)^(alpha)) = k*(1 + n) - (1 - delta)*k(-1) + c;
	
	% ii) Euler equation
	    # aux = 1 - delta + alpha*A*(k^(alpha - 1));
		# u_1 = 1/c^(sigma);
		# u_2 = 1/c(+1)^(sigma);
		u_1 = beta*u_2*aux;
		
end;


% -------------------------------------------------------
% COMPUTING THE STEADY STATE
% -------------------------------------------------------

initval; 
	k = kss;
	c = css;
end;
steady (solve_algo = 0);

 
check;   


% -------------------------------------------------------
% COMPUTING THE DYNAMICS
% -------------------------------------------------------

	initval; 
		% Initial Steady state
		k = ((alpha*1)/((1/beta) - (1 - delta)))^(1/(1 - alpha));
	end;
	
simul (stack_solve_algo=0, periods = 90);



	XY = oo_.endo_simul';
	save RCK_Dynamics_Results  XY;

	


% -------------------------------------------------------
% GRAPHING THE RESULTS: 
% -------------------------------------------------------
vec_time = 0:1:95; % 95


	load RCK_Dynamics_Results;
	A = 1;
	kss_0 = ((alpha*A)/((1/beta) - (1 - delta)))^(1/(1 - alpha)) 
	css_0 = A*(kss_0^(alpha)) - kss_0*(n + delta)
    vec_k = [kss_0*ones(11,1); XY(2:86,1)];
	vec_c = [css_0*ones(10,1); XY(1:86,2)];
	 
	
    
	vec_k_1 = vec_k;
	vec_c_1 = vec_c;
	
	
	

	figure;
	subplot(2,1,1);
	z = plot(vec_time, vec_k, 'b');
	set(z,'LineWidth',2.5);
	z = legend('Capital', 'Location', 'Best');
	set(z, 'FontSize', 12);
	vline(10);
	axis tight;
	ylabel('Stock of capital');
    title('RCK Model: Unexpected change in TFP', 'FontSize', 12);

	subplot(2,1,2);
	z = plot(vec_time, vec_c, 'r');
	set(z,'LineWidth',2.5);
	z = legend('Consumption', 'Location', 'Best');
	set(z, 'FontSize', 12);
	axis tight;
	vline(10);
	ylabel('Consumption');
	xlabel('Time', 'FontSize', 12);
	
    print -depsc2 -tiff RCK_ProblemSet_1.eps;



    


	
XY_1 = [vec_time', vec_k_1, vec_c_1];

save RCK_Dynamics_Results_1  XY_1;

aux_1 = [vec_time', vec_k_1, vec_c_1]
	
dynatype (RCK_Ex_1_OUTPUT_txt);
diary off;