clear
clc
%%%%%%%%%
% NOTE: for this program to run properly you also need to have:
%		Qz.m	  (this is under MATHWORKS copyright)
%		Qzdiv.m   (this is a C. Sims program) 
%		Qzswitch.m(this is a C. Sims program)
%%%%%%%%%
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 6.2: Identifying capital stock and inventories apart from each other')
%
disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause
%
%%%%%%%%%
%
clc
disp('_____________________________________________________')
disp('For the parameterization considered in the subsection, the')
disp('matrices for the about the steady state linear approximation are')
disp('_____________________________________________________')
%
%
%%%%%%%%%
%
% 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
%
% NOTE!!!!!: as we are simply copying the four-digits matrices in the paper,
% 	some degree of accuracy is lost in the LAMBDA,... matrices. To
%	have the very accuracy that is shown in subsection 6.2, use the
%	analitic expressions for the steady-state, GAMMA0 and GAMMA1.
%
GAMMA0 = [
%c(t),	k(t),	i(t),	Wk(t),	Wi(t),	log(theta(t))
-0.1222	0	0	-0.99	0	0
-0.1222 0	0	0	-0.99	0
1	1	1	0	0	-3.6313
-0.1235	0	0	0	0	0.0078
-0.1235	0	0	0	0	0.0022
0	0	0	0	0	1		]
%
% and,
%
GAMMA1 = [
%c(t-1),k(t-1),	i(t-1),	Wk(t-1),Wi(t-1),log(theta(t-1))
0	0	0	0	0	0
0	0	0	0	0	0
0	1.0101	1.0101	0	0	0
0	0.0002	-0.0003	1	0	0
0	-0.0003	0.003	0	1	0
0	0	0	0	0	0.95  	]

disp('   ')
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause
%
%%%%%%%%%
clc
disp('   ')
disp('_____________________________________________________')
disp('To study the stability of the linear system we need')
disp('to perform a QZ decomposition. These are LAMBDA and OMEGA...')
disp('_____________________________________________________')
%
%
[LAMBDA,OMEGA,Q,Z,V] = qz(GAMMA0,GAMMA1);
%
% Order the matrices such that eigenvalues are increasing along the
% main diagonal (with the QZDIV.M program),
%
[LAMBDA,OMEGA,Q,Z] = qzdiv(0.99^(-1/2),LAMBDA,OMEGA,Q,Z);
[LAMBDA,OMEGA,Q,Z] = qzdiv(0.96^(-1/2),LAMBDA,OMEGA,Q,Z);
LAMBDA
OMEGA
EIGENVALUES = diag(OMEGA)./diag(LAMBDA);
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause

clc
disp('_____________________________________________________')
disp('   ')
disp('From the pair (LAMBDA,OMEGA) we can extract the generalized eigenvalues')
disp('_____________________________________________________')
disp('   ')
disp('On the one hand, the finite eigenvalues are:  ')
disp(EIGENVALUES(1:5))
disp('   ')
disp('On the other hand, the infinite eigenvalue is:  ')
disp(EIGENVALUES(6))
disp('_____________________________________________________')
disp('   ')
disp('Hit any key to continue')
pause

clc
disp('_____________________________________________________')
disp('   ')
disp('The eigenvectors associated with the unstable eigenvalues')
disp('give us the stability conditions')
disp('_____________________________________________________')
disp('   ')
disp('For the finite eigenvalue is:  ')
disp(Z(:,5))
disp('   ')
disp('For the infinite eigenvalue is:  ')
disp(Z(:,6))
disp('_____________________________________________________')
disp('   ')
disp('And we are done. In the next step you can use these two')
disp('conditions to generate time series as described in the paper.')
