close all;
clear all;
clc;

%%%% Code for solving PBAR Model using the Blanchard-Kahn method %%%%%
%%%% There are two forward looking variables and three predetermined
%%%% variables, hence the partition of matrices, wherever applicable, is
%%%% done so as to have the upper block of three rows for predetermined variables and the
%%%% lower block for of two rows for non predetermined (jump) variables %%%

% Parameters

sigma   =   1;
b       = -1/sigma;
phi     = 0.89;
kappa   = 1;
rhoa    = 0.9;
mu1     = 0.5;
mu2     = 0.5;

%%%% The Predetermined variables are x = [delp a m] and the non
%%%% predetermined variables are y = [y R]

%% 1. y = y(+1) + b(r-infl) ** inflation is not (+1) as it is predetermined
%% 2. y = phi*y(-1) - phi*kappa*a(-1) + a
%% 3. R = mu1*infl(-1) + mu2*y - mu*kappa*a
%% 4. a = rho*a + e
%% 5. m+1 = y ** the dummy for y(-1). This is the additional eqn compared to dynare

%      infl a   m   y+1 r+1 
B   = [-b  0    0   1   0;
        0   kappa 0 0   0;
        0   mu2*kappa   0   0   0;
        0   1   0   0   0;
        0   0   1   0   0]; % LHS B = QTZ'

A   = [0    0   0   1   -b;
       0    phi*kappa   -phi    1   0;
       mu1  0   0   mu2     1;
       0    rhoa    0   0   0;
       0    0   0   1   0]; % RHS A = QSZ'

G   = [0;
       0;
       0;
       1;
       0];

[BB,AA,Q,Z] = qz(B,A);

% here B is A and A is B from book.  

[T,S,QS,ZS] = ordqz(BB,AA,Q,Z,'udo');

% in the output, BBS is T, AAS is S, Q' is Q, and Z is Z from book.
% eigenvalues are diagonal elements of AAS/BBS

QQ        = QS';

SS         = diag(S);
TT         = diag(T);

for i = 1:5;
    
    eigenvals(i) = SS(i)/TT(i);
    
end

eigenvals;

%% ** The problem when eigenvalues are checked seems to be that (1) if mu1 
%% is taken to be the usual 1.5, then the order of eigenvalues is wrong, as
%% well as there are more explosive eigenvalues than forward looking variables
%% and
%% (2) if mu1 is brought below 1, then also we have issues with ordering ** %%

ZST = ZS';
 
ZS11(1:3,1:3)=ZST(1:3,1:3);

ZS12 = ZST(1:3,4:5);

ZS21 = ZST(4:5,1:3);

ZS22 = ZST(4:5,4:5);

      %% N = inv(ZS22)*ZS21;

N   = ZS22\ZS21

B11 = B(1:3,1:3);

B12 = B(1:3,4:5);

A11 = A(1:3,1:3);

A12 = A(1:3,4:5);

Omega = inv(B11-B12*N)*(A11-A12*N);

G1 = G(1:3);

G2 = G(4:5);

S11 = S(1:3,1:3);

S12 = S(1:3,4:5);

S21 = S(4:5,1:3);

S22 = S(4:5,4:5);

T11 = T(1:3,1:3);

T12 = T(1:3,4:5);

T21 = T(4:5,1:3);

T22 = T(4:5,4:5);

QS11 = QS(1:3,1:3);

QS12 = QS(1:3,4:5);

QS21 = QS(4:5,1:3);

QS22 = QS(4:5,4:5);

L = inv(ZS22)*inv(S22)*(QS21*G1+QS22*G2);

D = inv(B11-B12*N)*(G1-A12*L);