% Calibrated parameter values

clear;
beta=0.9;
alpha=0.75;
sigma=1;
delta=0.3;
rho=0.95;

% Calculate steady-state values

kbar=((1-(1-delta)*beta)/(alpha*beta))^(1/(alpha-1));
cbar=kbar^alpha-delta*kbar;
ybar=kbar^alpha;

% Define state space matrices

A0=zeros(3,3);
A0(1,1)=1-(1-delta)*beta;
A0(1,2)=(1-(1-delta)*beta)*(alpha-1);
A0(1,3)=-sigma;
A0(2,2)=kbar;
A0(3,1)=1;

A1=zeros(3,3);
A1(1,3)=-sigma;
A1(2,1)=ybar;
A1(2,2)=alpha*ybar+(1-delta)*kbar;
A1(2,3)=-cbar;
A1(3,1)=rho;

B0=zeros(3,1);
B0(3,1)=1;

% Calculate alternative state space matrices

A=inv(A0)*A1;
B=inv(A0)*B0;

% Jordan decomposition of A

[p,lambda] = eig(A);
pstar = inv(p);

% Sort eigenvalues and eigenvectors in ascending order

val=diag(lambda);
t=sortrows([val p'],1);
lambda=diag(t(:,1));
p=t(:,2:4)';
pstar=inv(p);