var yt yh c x it ih xt xh k h r lt lh l p w;   % 12 variables;

varexo a;

parameters alphat alphah beta delta phi psi yt_ss yh_ss c_ss x_ss it_ss ih_ss r_ss l_ss lt_ss lh_ss h_ss;

alphat  =   0.35;
alphah  =   0.45;  
beta    =   0.99;
delta   =   0.025;
phi     =   1/3;
psi     =   -10;

%% SS values 

yt_ss   = 5;
yh_ss   = 1;
x_ss    = 1;
c_ss    = 3.4363;
h_ss    = 12.6905;
it_ss   = 1.2464;
ih_ss   = 0.3173;
r_ss    = 1/beta; 
l_ss    = 1.6922;
lt_ss   = 1.4494;
lh_ss   = 0.2428;


model(linear);

k = (1-delta)*k(-1)+delta*it;                           
xt = -psi*delta*(it-k(-1));                                     

yt = alphat*k(-1)+(1-alphat)*lt;
w = yt-lt;
r = ((r_ss-1+delta)/r_ss)*(yt-k(-1))+((1-delta)/r_ss)*xt-xt(-1);

h = (1-delta)*h(-1)+delta*ih;
xh = -psi*delta*(ih-h(-1));

yh = a+alphah*h(-1)+(1-alphah)*lh;                                    
w = p+yh-lh;                                                        
r = ((r_ss-1+delta)/r_ss)*(p+yh-h(-1))+((1-delta)/r_ss)*xh-xh(-1); 

w = c+phi*l;
x(+1)+p(+1) = r(+1)+x+p;
c(+1) = r(+1)+c;

yh = x; 
yt_ss*yt = c_ss*c+it_ss*it+ih_ss*ih;
l_ss*l = lt_ss*lt+lh_ss*lh;
end;


initval;
yt = 0;
yh = 0;
c  = 0;
x  = 0;
it = 0;
ih = 0;
xt = 0;
xh = 0;
k  = 0;
h  = 0;
r  = 0;
a  = 0;
lt = 0;
lh = 0;
l  = 0;
p  = 0;
w  = 0;
end;
steady;
check;

endval;
yt = 0;
yh = 0;
c  = 0;
x  = 0;
it = 0;
ih = 0;
xt = 0;
xh = 0;
k  = 0;
h  = 0;
r  = 0;
a  = 0.1;
lt = 0;
lh = 0;
l  = 0;
p  = 0;
w  = 0;
end;
steady;
check;

simul(periods=2100);


figure(1)
subplot(3,2,1);
hold on;
 plot(yt(1:40), '--');
subplot(3,2,2);
hold on;
plot(yh(1:40) , '--');
subplot(3,2,3);
hold on;
 plot(c(1:40), '--');
subplot(3,2,4);
hold on;
plot(x(1:40), '--');
subplot(3,2,5);
hold on;
 plot(k(1:40), '--');
 subplot(3,2,6);
hold on;
plot(h(1:40), '--');

figure(2)
subplot(3,2,1);
hold on;
 plot(it(1:40), '--');
subplot(3,2,2);
hold on;
plot(ih(1:40), '--');
 subplot(3,2,3);
hold on;
 plot(l(1:40), '--');
subplot(3,2,4);
hold on;
plot(lt(1:40), '--');
subplot(3,2,5);
hold on;
 plot(lh(1:40), '--');
 subplot(3,2,6);
hold on;
plot(w(1:40), '--');

figure(3)
subplot(3,2,1);
hold on;
 plot(p(1:40), '--');
subplot(3,2,2);
hold on;
plot(xt(1:40), '--');
 subplot(3,2,3);
hold on;
 plot(xh(1:40), '--');
subplot(3,2,4);
hold on;
plot(i,r_ss*10000*r(2:40),'--');