%set parameter values
beta=0.99;eta=2;gamma=0;alpha1=0.36;alpha2=0.36;
delta=0.025; b=0.8; rho1=0.95; rho2=0.95;
param=[beta;eta;gamma;alpha1;alpha2;delta;b];
xint=ones(8,1); xint(7:8)=5;
ys=broyden('ss_tsjhrbc',xint,param);
if abs(ys)-abs(real(ys))>1;
    warning('Try Different Initial Valuse');
else;
    ys=real(ys);
end;
xint=ys;
c1ss=xint(1);
lam1ss=xint(2);
dss=xint(3);
lam2ss=xint(4);
h1ss=xint(5);
h2ss=xint(6);
k1ss=xint(7);
k2ss=xint(8);
kss=k1ss+k2ss;
Iss=delta*kss;
ome2ss=alpha1*lam1ss*(k1ss/h1ss)^(alpha1-1);
ome1ss=beta*ome2ss/(1-beta*(1-delta));
%loglinearization
%f0(t)=[c h1 h2 k1 k2 ome2 I] we have c=-lam1 lam2=ome1
%s0(t)=[ome1 k2];
%bigA is 5 by 5 and bigB is 5 by 2 and C is 
bigA=zeros(7,7);bigB=zeros(7,2);bigC=zeros(7,2);
bigA(1,1)=1;big(1,2)=gamma+alpha1;bigA(1,4)=-alpha1;bigC(1,1)=1;

bigA(2,3)=gamma+alpha2;bigA(2,5)=-alpha2;bigB(2,1)=0;bigC(2,2)=1;

bigA(3,1)=-1;bigA(3,2)=1-alpha1; bigA(3,4)=alpha1-1;big(3,6)=-1;bigC(3,1)=-1;

bigA(4,3)=1-alpha2;bigA(4,5)=alpha2-1;bigA(4,6)=-1;bigB(4,1)=-1;bigC(4,2)=-1;

bigA(5,1)=1; bigA(5,2)=alpha1-1; bigA(5,4)=-alpha1;bigC(5,1)=1;

bigA(6,3)=alpha2-1;bigA(6,5)=-alpha2;bigA(6,7)=1; bigC(6,2)=1;

bigA(7,4)=k1ss; bigA(7,5)=k2ss; bigB(7,2)=kss;


%Dynamic System;
bigD=zeros(2,2);bigF=zeros(2,7);bigG=zeros(2,2);bigH=zeros(2,7);bigJ=zeros(2,2);
bigD(1,1)=beta*(1-delta)*ome1ss; bigF(1,6)=beta*ome2ss; bigG(1,1)=ome1ss;
bigD(2,2)=kss;bigF(2,7)=Iss; bigG(2,2)=(1-delta)*kss;

bigP=[rho1 0;0 rho2];
%substitute f0 by f0=A^(-1)*B*s0+A^(-1)*B*v;

%bigX1*Ec(s0(t+1))+bigY1*v(t)=bigX2*s0(t)+bigY2*v(t);
bigX1=bigD+bigF*inv(bigA)*bigB;
bigY1=bigF*inv(bigA)*bigC*bigP;
bigX2=bigG+bigH*inv(bigA)*bigB;
bigY2=bigJ+bigH*inv(bigA)*bigC;

%Then we have Et{s0(t+1)}=bigK*s0(t)+bigL*v(t);
bigK=inv(bigX1)*bigX2;
bigL=inv(bigX1)*(bigY2-bigY1);

[V_K,eigK]=eig(bigK);
eigK=diag(eigK)