clear all;

% benchmark model

%parameters
alpha=0.4;
beta1=0.99;
beta2=0.5;
n=0.5;
ktotal=1;
rho=0.95;
sig=2;

%steady state
brat=(beta1/beta2)^(1/(1-alpha))
kbar2=ktotal/((1-n)*brat+n)
kbar1=brat*kbar2
ybar1=kbar1^alpha
ybar2=kbar2^alpha
ybartot=(1-n)*ybar1+n*ybar2
qbar=alpha*beta1*kbar1^(alpha-1)/(1-beta1)
rbar=1/beta1
ucbar=(1-beta1)*qbar
cbar1=(1+(n/(1-n))*alpha*beta1/brat)*ybar1
cbar2=(1-alpha*beta1)*ybar2
mubar=(beta1-beta2)/cbar2^sig
gamma=n*kbar2/((1-n)*kbar1)

disp('Press enter to keep going');
pause;

amat=zeros(4);
bmat=zeros(4);
nstate=2;

zz=1;
kk2=2;
cc2=3;
cc1=4;


amat(zz,zz)=1;
amat(kk2,zz)=(n/(1-n))*alpha*beta2*ybar2;
amat(kk2,kk2)=(n/(1-n))*alpha^2*beta2*ybar2;
amat(kk2,cc2)=-(n/(1-n))*alpha*beta2*ybar2*sig;
amat(cc2,zz)=alpha*beta2*ybar2;
amat(cc2,kk2)=(alpha^2)*beta2*ybar2;
amat(cc2,cc2)=-alpha*beta2*ybar2*sig;
amat(cc1,kk2)=(1-alpha)*(1+gamma);
amat(cc1,cc2)=sig;
amat(cc1,cc1)=-sig;

bmat(zz,zz)=rho;
bmat(kk2,zz)=-ybar1;
bmat(kk2,kk2)=alpha*gamma*ybar1;
bmat(kk2,cc2)=-alpha*beta2*ybar2*sig*n/(1-n);
bmat(kk2,cc1)=cbar1;
bmat(cc2,zz)=ybar2;
bmat(cc2,kk2)=alpha*ybar2;
bmat(cc2,cc2)=-cbar2-alpha*beta2*ybar2*sig;
bmat(cc1,cc2)=sig;
bmat(cc1,cc1)=-sig;


[fmat,pmat]=solab(amat, bmat, nstate)

format compact;
disp(' ')
disp('      Date     khat1     khat2     yhat1     yhat2      yhat      zhat')

statevar=[0.01
          0];

tt=1; 
while tt<21; 
    

zhat=statevar(zz,1);

khat2=statevar(kk2,1);

endvar = fmat*statevar;

chat2=endvar(1,1);
chat1=endvar(2,1);

khat1=-gamma*khat2;

yhat1=zhat+alpha*khat1;
yhat2=zhat+alpha*khat2;

yhattot=(1-n)*(ybar1/ybartot)*yhat1+n*(ybar2/ybartot)*yhat2;

disp([tt,khat1*100,khat2*100,yhat1*100,yhat2*100,yhattot*100,zhat*100])

statevar = pmat*statevar;
tt=tt+1;
end;