clear
clc
%%%%%%%%%
% NOTE: for this program to run properly you also need to have:
%		FSectio4.m
%		Csolve.m
%%%%%%%%%
disp('_____________________________________________________')
disp('Solving Nonlinear Rational Expectations Models')
disp('by Eigenvalue-Eigenvector Decompositions')
%
disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('A. Novales, E. Domínguez, J.J. Pérez and J. Ruiz')
%
disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Section 4: Solving a standard, stochastic growth model')
disp('This program generates a single realization of the variables in the economy')
%
disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause
%
%%%%%%%%%
%
clc
disp('_____________________________________________________')
disp('These are the structural PARAMETERS for this economy...')
disp('_____________________________________________________')
%
sigma		= 1.5;		% Relative risk aversion coefficient
beta		= 0.99;		% Utility discount factor
rho		= 0.95;		% Autorregressive coeficient to technology
alpha		= 0.36;		% Capital share
delta		= 0.025;	% Depreciation rate
sigma_eps	= 0.00721;	% Standard deviation of the noise
disp('   ')
disp('[sigma, beta, rho, alpha, delta, sigma_eps]=')
disp('   ')
disp(['[ ' num2str(sigma) ', ' num2str(beta) ', ' num2str(rho) ', ' num2str(alpha) ', ' num2str(delta) ', ' num2str(sigma_eps) ' ]'])
disp('   ')
disp('and the sample size...')
T		= 150		% Time series lenght 
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause
%
%%%%%%%%%
%
clc
disp('_____________________________________________________')
disp('The DETERMINISTIC STEADY STATE is...')
disp('_____________________________________________________')
%
theta_ss 	= 1;				% Autoregressive
y_k_ss 		= (1/alpha) * ( 1/beta - 1 + delta );% Output / Capital stock 
N_ss 		= 1/3;				% Employment
k_ss 		= ( (  y_k_ss / theta_ss  ) ^ ( 1 / (alpha-1) ) ) * N_ss;
						% Capital stock 
y_ss 		= y_k_ss * k_ss;		% Output
x_ss 		= k_ss * delta;			% Investment
c_ss 		= y_ss - x_ss;			% Consumption
AN		= ( c_ss^(-sigma) ) * (1-alpha) * y_ss / N_ss;% Labor utility parameter
W_ss		= c_ss^(-sigma)/beta;		% Expectation
disp('   ')
disp('[c_ss N_ss W_ss k_ss]=')
disp('   ')
disp(['[ ' num2str(c_ss) ', ' num2str(N_ss) ', ' num2str(W_ss) ', ' num2str(k_ss) ' ]'])
disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause
%
%%%%%%%%%
%
% Linearized System. The following notation:
%	GAMMA0 * y(t) = GAMMA1 * y(t-1) + PSI * z(t) + PI * eta(t)
%%%%%%%%%
%
clc
disp('_____________________________________________________')
disp('The Matrices for the Linear Approximation...')
disp('_____________________________________________________')
%
% where
%
GAMMA0 = [
% 		c(t),							N(t), 							W(t),	 k(t), 		log(theta(t))
  	 	1,						-k_ss^(alpha)*(1-alpha)*N_ss^(-alpha)				0,	1,	-k_ss^(alpha)*N_ss^(1-alpha)
sigma*c_ss^(-sigma-1)*k_ss^(alpha)*(1-alpha)*N_ss^(-alpha),	c_ss^(-sigma)*k_ss^(alpha)*(1-alpha)*alpha*N_ss^(-alpha-1),	0,	0,-c_ss^(-sigma)*k_ss^(alpha)*(1-alpha)*N_ss^(-alpha)
	-sigma*c_ss^(-sigma-1),						0 							-beta, 	0,			0,
sigma*c_ss^(-sigma-1)*(alpha*k_ss^(alpha-1)*N_ss^(1-alpha)+1-delta), -c_ss^(-sigma)*alpha*k_ss^(alpha-1)*(1-alpha)*N_ss^(-alpha),0,	0,-alpha*c_ss^(-sigma)*k_ss^(alpha-1)*N_ss^(1-alpha) 
		0, 							0 							0, 	0, 			1,

																								]
%
% and,
%
GAMMA1 = [
% c(t-1),N(t-1),W(t-1),	 			k(t-1), 				log(theta(t))
   0, 	 0, 	0, 	alpha*k_ss^(alpha-1)*N_ss^(1-alpha)+1-delta, 			0
   0, 	 0, 	0, 	c_ss^(-sigma)*alpha*k_ss^(alpha-1)*(1-alpha)*N_ss^(-alpha),	0
   0, 	 0, 	0, 				0, 					0
   0, 	 0, 	-1, 	alpha*(alpha-1)*c_ss^(-sigma)*k_ss^(alpha-2)*N_ss^(1-alpha) 	0	 
   0, 	 0, 	0, 				0, 					rho
												]

disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause
%
% and also,
%
PSI	= [ 0, 0, 0, 0, 1]';
PI	= [ 0, 0, 0, -1, 0]';
%
%%%%%%%%%
clc
disp('   ')
disp('_____________________________________________________')
disp('To study the stability of the linear system we need the eigenvalues')
disp('_____________________________________________________')
disp('   ')
%
%
[  P , LAMBDA  ]   =   eig(inv(GAMMA0)*GAMMA1);
%
% Order the matrices such that eigenvalues are increasing along the
% main diagonal. For QZ decompositions you have the program QZDIV.M
%
[diag_LAMBDA,Index] 	= sort(diag(LAMBDA));
LAMBDA			= diag(diag_LAMBDA)
P 			= P(:,Index)
disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause

% In P, the eigenvectors and in LAMBDA the eigenvalues. 
% The inverse of P:
clc
disp('_____________________________________________________')
disp('then, the unstable eigenvalue (in absolute value) is:  ')
lambda = [max(abs(diag(LAMBDA)))]
disp('and the matrix with the left eigenvectors...')
Pinv = inv(P)

% so the stability condition, P_5_bullet, arises from the fifth row of Pinv,
%	the one associated with the unstable eigenvalue
%
disp('  ')
disp('so the stability condition arises from the fifth row, being...')
disp('[ c(t)-c_ss,   N(t)-N_ss,   W(t)-W_ss, k(t)-k_ss, log(theta(t)) ]=')
disp('  ')
P_5_bullet = Pinv(5,:);
disp([ P_5_bullet ])
%
disp('  ')
disp('_____________________________________________________')
disp('Hit any key to generate time series. Then wait a minute, please...')
pause
clc
%
%%%%%%%%%
% And finally one can generate time series for all the variables 
%	in the economy
%

% Exogenous technology...
epsilon 	= randn(1,T+1000)*sigma_eps;	% This is the innovation
epsilon		= epsilon(1,1001:T+1000);
% For IRF's, skip the two previous lines and let them be...
%epsilon=zeros(1,T);
%epsilon(2)=1;
theta(1) 	= 1;	% Set the initial value to the steady state
for t=2:T
	theta(t) = exp( rho * log(theta(t-1)) + epsilon(t) );
end

% and for consumption, labor, capital, 
%	the expectation and the expectation error...
c(1)	= c_ss;
N(1)	= N_ss;
k(1)	= k_ss;
W(1)	= W_ss;

for t=2:T
% for consumption one have to solve the non-linear equation...

  c(t) = csolve('FSectio4',c(t-1),[W_ss c_ss k_ss AN beta sigma alpha k(t-1) theta(t) P_5_bullet delta]);

% then labor and capital arise from the optimality conditions...
  N(t) = ( AN *c(t)^(sigma) /...
			(theta(t)*k(t-1)^(alpha)*(1-alpha)) )^(-1/alpha);

  k(t) = -c(t) + theta(t) * k(t-1)^(alpha) * N(t)^(1-alpha) +...
							(1-delta) * k(t-1);

% and also the expectation...

  W(t) = (1/beta) * c(t)^(-sigma);

% Then, the expectation error from its definition...

  eta(t) = -W(t-1) + c(t)^(-sigma) * ( 1 - delta + ...
			alpha*theta(t) * k(t-1)^(alpha-1) * N(t)^(1-alpha) );

end

clc
disp('_____________________________________________________')
plot_or_not	= input('Do you want to plot the time series (y/n)?  ','s');
disp('  ')
disp('_____________________________________________________')
disp('Hit any key to generate time series...')

if plot_or_not == 'y'
	plot(c),grid,title('Consumption')
	figure
	plot(k),grid,title('Capital')
	figure
	plot(N),grid,title('Labor')
	figure
	plot(W),grid,title('Expectation W')
	figure
	plot(eta),grid,title('Expectation Error')
end

clc
disp('_____________________________________________________')
disp('The process has finished.')
disp('  ')
disp('_____________________________________________________')
