%Steady state prices
%Parameters
beta=.99;
delta=.05;
sigmava1=.5;
sigmaiy1=10;
sigmadm1=1.59;
sigmava2=.5;
sigmaiy2=10;
sigmadm2=1.59;
sigmavi2=.5;
sigmava3=.5;
sigmavi3=.5;
sigmaic=10;
c1=1; 
c2=1;
n=.04;
nstar=1-n;
x=2*n^(1/sigmadm1)/(n^(1/sigmadm1)+(1-n)^(1/sigmadm1));
xstar=2*nstar^(1/sigmadm1)/(nstar^(1/sigmadm1)+(1-nstar)^(1/sigmadm1));
theta1=.69;
theta1star=.69;
omega1=.5;
omega1star=.5;
theta2=.69;
theta2star=.69;
omega2=.5;
omega2star=.5;
gamma2=.5;
gamma2star=.5;
theta3=.69;
theta3star=.69;
gamma3=.5;
gamma3star=.5;
kappac=1;
kappah=1.02;
sigmac=.2;
sigmah=.1;
sigmahh=6.36;
psiw=.25;
psihab=.7;
pdotss=0;
pdotssstar=0;
chid1=127.8;
chid1star=127.8; 
chid2=127.8; 
chid2star=127.8; 
chix1=180; 
chix1star=180;
chix2=180; 
chix2star=180; 
chi3=127.8; 
chi3star=127.8; 
chib=.01;
epsilon1d=.25; 
epsilon1dstar=.25; 
epsilon2d=.25; 
epsilon2dstar=.25; 
epsilon1x=.5; 
epsilon1xstar=.5; 
epsilon2x=.5; 
epsilon2xstar=.5; 
epsilon3=.25; 
epsilon3star=.25; 
epsilonk=.7; 
epsilonw=.9;
sigmaz=.1; 
chik=120; 
chiz=.02; 
thetar=.5; 
thetap1=.5; 
thetap2=.25;
thetap3=.125;
thetap4=.0625
thetay=.0625; 
shareh1=.6;
shareh2=.6;
shareh3=.6;
shareh1star=.6;
shareh2star=.6;
shareh3star=.6;
rho1=.9;
rho2=.9;
rho3=.9;
rho1star=.9;
rho2star=.9;
rho3star=.9;
sigma1=1;
sigma2=1; 
sigma3=1;
sigma1star=1;
sigma2star=1;
sigma3star=1;


%Find expenditure shares, use Russian doll approach
%note, the way the code is written now is OK as long as value added
%technology does not change between stages, this is because these shares 
%were calculated under the assumption that marginal cost of the value 
%added portion is constant across stages, if the value added technology
%does change then the marginal cost of value added should differ across
%stages, so these shares need to be recalculated (add another layer to the
%Russian doll to account for changing technology across stages but costant
%wage and rental rates.
%Create a grid for mcdiff
mcdiffmin=.01;
mcdiffmax=100;
int=.001;
mcdiffgrid=[mcdiffmin:int:mcdiffmax]';
i=1
diff=10;

while diff>.001;
mcdiff=mcdiffgrid(i,1);
fx=sigmaiy1/(sigmaiy1-1);
fy=(1+c1)*sigmaiy1/(sigmaiy1-1)*mcdiff;
fxstar=sigmaiy1/(sigmaiy1-1);
fystar=(1+c1)*sigmaiy1/(sigmaiy1-1)/mcdiff;
a1=((x*omega1)^sigmadm1*fx^(1-sigmadm1)+(1-x*omega1)^sigmadm1*fy^(1-sigmadm1))^(1/(1-sigmadm1));
a1star=((xstar*omega1star)^sigmadm1*fxstar^(1-sigmadm1)+(1-xstar*omega1star)^sigmadm1*fystar^(1-sigmadm1))^(1/(1-sigmadm1));
shareii2=(1-gamma2)^sigmavi2*a1^(1-sigmavi2)/(gamma2^sigmavi2+(1-gamma2)^sigmavi2*a1^(1-sigmavi2));
shareii2star=(1-gamma2star)^sigmavi2*a1star^(1-sigmavi2)/(gamma2star^sigmavi2+(1-gamma2star)^sigmavi2*a1star^(1-sigmavi2));
shareva2=1-shareii2;
shareva2star=1-shareii2star;
shared1=(x*omega1)^sigmadm1/((x*omega1)^sigmadm1+(1-x*omega1)^sigmadm1*((1+c1)*mcdiff)^(1-sigmadm1));
shared1star=(xstar*omega1star)^sigmadm1/((xstar*omega1star)^sigmadm1+(1-xstar*omega1star)^sigmadm1*((1+c1)/mcdiff)^(1-sigmadm1));
sharem1=1-shared1;
sharem1star=1-shared1star;
fx=1;
fy=a1;
fxstar=1;
fystar=a1star;
b2=(gamma2^sigmavi2*fx^(1-sigmavi2)+(1-gamma2)^sigmavi2*fy^(1-sigmavi2))^(1/(1-sigmavi2));
b2star=(gamma2star^sigmavi2*fxstar^(1-sigmavi2)+(1-gamma2star)^sigmavi2*fystar^(1-sigmavi2))^(1/(1-sigmavi2));
fx=sigmaiy1/(sigmaiy1-1)*b2;
fy=(1+c2)*sigmaiy1/(sigmaiy1-1)*mcdiff*b2star;
fxstar=sigmaiy1/(sigmaiy1-1)*b2star;
fystar=(1+c2)*sigmaiy1/(sigmaiy1-1)*b2/mcdiff;
a2=((x*omega2)^sigmadm2*fx^(1-sigmadm2)+(1-x*omega2)^sigmadm2*fy^(1-sigmadm2))^(1/(1-sigmadm2));
a2star=((xstar*omega2star)^sigmadm2*fxstar^(1-sigmadm2)+(1-xstar*omega2star)^sigmadm2*fystar^(1-sigmadm2))^(1/(1-sigmadm2));
shareii3=(1-gamma3)^sigmavi3*a2^(1-sigmavi3)/(gamma3^sigmavi3+(1-gamma3)^sigmavi3*a2^(1-sigmavi3));
shareii3star=(1-gamma3star)^sigmavi3*a2star^(1-sigmavi3)/(gamma3star^sigmavi3+(1-gamma3star)^sigmavi3*a2star^(1-sigmavi3));
shareva3=1-shareii3;
shareva3star=1-shareii3star;
shared2=(x*omega2)^sigmadm2*b2^(1-sigmadm2)/((x*omega2)^sigmadm2*b2^(1-sigmadm2)+(1-x*omega2)^sigmadm2*((1+c2)*mcdiff*b2star)^(1-sigmadm2));
shared2star=(xstar*omega2star)^sigmadm2*b2star^(1-sigmadm2)/((xstar*omega2star)^sigmadm2*b2star^(1-sigmadm2)+(1-xstar*omega2star)^sigmadm2*((1+c2)*b2/mcdiff)^(1-sigmadm2));
sharem2=1-shared2;
sharem2star=1-shared2star;

%Now calculate the steady state prices 
p3=1;
p3star=1;
Psi3=sigmaic/(sigmaic-1);
Psi3star=sigmaic/(sigmaic-1);
mc3=p3/Psi3;
mc3star=p3star/Psi3star;
p2=(1/(1-gamma3))^(sigmava3/(1-sigmava3))*(shareii3)^(1/(1-sigmava3))*mc3;
p2star=(1/(1-gamma3star))^(sigmava3/(1-sigmava3))*(shareii3star)^(1/(1-sigmava3))*mc3star;
mc3va=(1/(gamma3))^(sigmava3/(1-sigmava3))*(shareva3)^(1/(1-sigmava3))*mc3;
mc3vastar=(1/(gamma3star))^(sigmava3/(1-sigmava3))*(shareva3star)^(1/(1-sigmava3))*mc3star;
pm2=(1/(1-x*omega2))^(sigmadm2/(1-sigmadm2))*(sharem2)^(1/(1-sigmadm2))*p2;
pm2star=(1/(1-xstar*omega2star))^(sigmadm2/(1-sigmadm2))*(sharem2star)^(1/(1-sigmadm2))*p2star;
pd2=(1/(x*omega2))^(sigmadm2/(1-sigmadm2))*(shared2)^(1/(1-sigmadm2))*p2;
pd2star=(1/(xstar*omega2star))^(sigmadm2/(1-sigmadm2))*(shared2star)^(1/(1-sigmadm2))*p2star;
px2star=pm2/(1+c2);
px2=pm2star/(1+c2);
Psi2d=sigmaiy2/(sigmaiy2-1);
Psi2x=sigmaiy2/(sigmaiy2-1);
Psi2dstar=sigmaiy2/(sigmaiy2-1);
Psi2xstar=sigmaiy2/(sigmaiy2-1);
mc2=pd2/Psi2d;
mc2star=pd2star/Psi2dstar;
p1=(1/(1-gamma2))^(sigmava2/(1-sigmava2))*(shareii2)^(1/(1-sigmava2))*mc2;
p1star=(1/(1-gamma2star))^(sigmava2/(1-sigmava2))*(shareii2star)^(1/(1-sigmava2))*mc2star;
mc2va=(1/(gamma2))^(sigmava2/(1-sigmava2))*(shareva2)^(1/(1-sigmava2))*mc2;
mc2vastar=(1/(gamma2star))^(sigmava2/(1-sigmava2))*(shareva2star)^(1/(1-sigmava2))*mc2star;
pm1=(1/(1-x*omega1))^(sigmadm1/(1-sigmadm1))*(sharem1)^(1/(1-sigmadm1))*p1;
pm1star=(1/(1-xstar*omega1star))^(sigmadm1/(1-sigmadm1))*(sharem1star)^(1/(1-sigmadm1))*p1star;
pd1=(1/(x*omega1))^(sigmadm1/(1-sigmadm1))*(shared1)^(1/(1-sigmadm1))*p1;
pd1star=(1/(xstar*omega1star))^(sigmadm1/(1-sigmadm1))*(shared1star)^(1/(1-sigmadm1))*p1star;
px1star=pm1/(1+c1);
px1=pm1star/(1+c1);
Psi1d=sigmaiy1/(sigmaiy1-1);
Psi1x=sigmaiy1/(sigmaiy1-1);
Psi1dstar=sigmaiy1/(sigmaiy1-1);
Psi1xstar=sigmaiy1/(sigmaiy1-1);
mc1=pd1/Psi1d;
mc1star=pd1star/Psi1dstar;
r=1/beta-1+delta;
rstar=1/beta-1+delta;
w=(1/theta1)^(sigmava1/(1-sigmava1))*(shareh1)^(1/(1-sigmava1))*mc1;
wstar=(1/theta1star)^(sigmava1/(1-sigmava1))*(shareh1star)^(1/(1-sigmava1))*mc1star;
diff=abs(mc1/mc1star-mcdiff)
i=i+1
end;


%Calculate steady state quantities, for this we need to solve for the
%appropriate capital stock k by a numerical solution. 

%Grid for k and kstar
kmin=1;
kmax=100;
int=.001;
kgrid=[kmin:int:kmax]';
kgridstar=[kmin:int:kmax]';
i=1;
istar=1;
diff=10;
diffstar=10;
pop=1;
popstar=1;

while diff>.001 | diffstar>.001; 
k=kgrid(i,1);
kstar=kgridstar(istar,1);
c=kappac^(1/(1+psihab*sigmac-psihab));
cstar=kappac^(1/(1+psihab*sigmac-psihab));
I=delta*k;
Istar=delta*kstar;
z=1;
zstar=1;
ks=z*k;
ksstar=zstar*kstar;
y3=c+I;
y3star=cstar+Istar;
VA3=gamma3^sigmavi3*(mc3va/mc3)^(-sigmavi3)*y3;
y2bar=(1-gamma3)^sigmavi3*(p2/mc3)^(-sigmavi3)*y3;
VA3star=gamma3star^sigmavi3*(mc3vastar/mc3star)^(-sigmavi3)*y3star;
y2barstar=(1-gamma3star)^sigmavi3*(p2star/mc3star)^(-sigmavi3)*y3star;
yD2bar=(x*omega2)^sigmadm2*(pd2/p2)^(-sigmadm2)*y2bar;
yM2bar=(1-x*omega2)^sigmadm2*(pm2/p2)^(-sigmadm2)*y2bar;
yD2barstar=(xstar*omega2star)^sigmadm2*(pd2star/p2star)^(-sigmadm2)*y2barstar;
yM2barstar=(1-xstar*omega2star)^sigmadm2*(pm2star/p2star)^(-sigmadm2)*y2barstar;
yX2bar=(1+c2)*yM2barstar;
yX2barstar=(1+c2)*yM2bar;
y2=yD2bar+yX2bar;
y2star=yD2barstar+yX2barstar;
VA2=gamma2^sigmavi2*(mc2va/mc2)^(-sigmavi2)*y2;
y1bar=(1-gamma2)^sigmavi2*(p1/mc2)^(-sigmavi2)*y2;
VA2star=gamma2star^sigmavi2*(mc2vastar/mc2star)^(-sigmavi2)*y2star;
y1barstar=(1-gamma2star)^sigmavi2*(p1star/mc2star)^(-sigmavi2)*y2star;
yD1bar=(x*omega1)^sigmadm1*(pd1/p1)^(-sigmadm1)*y1bar;
yM1bar=(1-x*omega1)^sigmadm1*(pm1/p1)^(-sigmadm1)*y1bar;
yD1barstar=(xstar*omega1star)^sigmadm1*(pd1star/p1star)^(-sigmadm1)*y1barstar;
yM1barstar=(1-xstar*omega1star)^sigmadm1*(pm1star/p1star)^(-sigmadm1)*y1barstar;
yX1bar=(1+c1)*yM1barstar;
yX1barstar=(1+c1)*yM1bar;
y1=yD1bar+yX1bar;
y1star=yD1barstar+yX1barstar;
h3d=theta3^sigmava3*(w/mc3va)^(-sigmava3)*VA3;
k3d=(1-theta3)^sigmava3*(r/mc3va)^(-sigmava3)*VA3;
h3dstar=theta3star^sigmava3*(wstar/mc3vastar)^(-sigmava3)*VA3star;
k3dstar=(1-theta3star)^sigmava3*(rstar/mc3vastar)^(-sigmava3)*VA3star;
h2d=theta2^sigmava2*(w/mc2va)^(-sigmava2)*VA2;
k2d=(1-theta2)^sigmava2*(r/mc2va)^(-sigmava2)*VA2;
h2dstar=theta2star^sigmava2*(wstar/mc2vastar)^(-sigmava2)*VA2star;
k2dstar=(1-theta2star)^sigmava2*(rstar/mc2vastar)^(-sigmava2)*VA2star;
h1d=theta1^sigmava1*(w/mc1)^(-sigmava1)*y1;
k1d=(1-theta1)^sigmava1*(r/mc1)^(-sigmava1)*y1;
h1dstar=theta1star^sigmava1*(wstar/mc1star)^(-sigmava1)*y1star;
k1dstar=(1-theta1star)^sigmava1*(rstar/mc1star)^(-sigmava1)*y1star;
hs=kappah*((sigmahh-1)/sigmahh)^sigmah*(w^sigmah);
hsstar=kappah*((sigmahh-1)/sigmahh)^sigmah*(wstar^sigmah);
pop=(h1d+h2d+h3d);
popstar=(h1dstar+h2dstar+h3dstar);
hs=(h1d+h2d+h3d)/pop;
hsstar=(h1dstar+h2dstar+h3dstar)/popstar;
diffstar=abs(k1dstar+k2dstar+k3dstar-kstar);
if diffstar<.001;
istar=istar;
else;
istar=istar+1;
end;
diff=abs(k1d+k2d+k3d-k);
if diff<.001;
i=i;
else;
i=i+1;
end;
end;

Lambda=1;
Lambdastar=1;
rrf=(1/beta)-1;
rrfstar=(1/beta)-1;
Q=1;
Xi=hs^((sigmah+1)/sigmah)/(1-beta*(1-psiw));
Xistar=hsstar^((sigmah+1)/sigmah)/(1-beta*(1-psiw));
WTheta=Lambda*hs/(1-beta*(1-psiw));
WThetastar=Lambdastar*hsstar/(1-beta*(1-psiw));
kappah=(sigmahh/(w*(sigmahh-1)))^sigmah;
kappahstar=(sigmahh/(wstar*(sigmahh-1)))^sigmah;
wsquig=(w^(sigmahh/sigmah)*kappah^(-1/sigmah)*sigmahh/(sigmahh-1)*Xi/WTheta)^(1/(1+sigmahh/sigmah));
wsquigstar=(wstar^(sigmahh/sigmah)*kappahstar^(-1/sigmah)*sigmahh/(sigmahh-1)*Xistar/WThetastar)^(1/(1+sigmahh/sigmah));
xiw=(1+pdotss);
xiwstar=(1+pdotssstar);
Upsilon1d=0;
xi1d=0;
Upsilon1x=0;
xi1x=0;
Upsilon2d=0;
xi2d=0;
Upsilon2x=0;
xi2x=0;
Upsilon3=0;
xi3=0;
Upsilon1dstar=0;
xi1dstar=0;
Upsilon1xstar=0;
xi1xstar=0;
Upsilon2dstar=0;
xi2dstar=0;
Upsilon2xstar=0;
xi2xstar=0;
Upsilon3star=0;
xi3star=0;
b=0;
bf=0;
bfstar=0;
pi1=pd1*yD1bar+px1*yX1bar-w*h1d-r*k1d;
pi2=pd2*yD2bar+px2*yX2bar-w*h2d-r*k2d-pd1*yD1bar-pm1*yM1bar;
pi3=p3*y3-w*h3d-r*k3d-pd2*yD2bar-pm2*yM2bar;
pi1star=pd1star*yD1barstar+px1star*yX1barstar-wstar*h1dstar-rstar*k1dstar;
pi2star=pd2star*yD2barstar+px2star*yX2barstar-wstar*h2dstar-rstar*k2dstar-pd1star*yD1barstar-pm1star*yM1barstar;
pi3star=p3star*y3star-wstar*h3dstar-rstar*k3dstar-pd2star*yD2barstar-pm2star*yM2barstar;
pi1=(Psi1x-1)*mc1*yX1bar+(Psi1d-1)*mc1*yD1bar;
pi2=(Psi2x-1)*mc2*yX2bar+(Psi2d-1)*mc2*yD2bar;
pi3=(Psi3-1)*mc3*y3;
pi1star=(Psi1xstar-1)*mc1star*yX1barstar+(Psi1dstar-1)*mc1star*yD1barstar;
pi2star=(Psi2xstar-1)*mc2star*yX2barstar+(Psi2dstar-1)*mc2star*yD2barstar;
pi3star=(Psi3star-1)*mc3star*y3star;

OG1=1;
OG2=1;
OG3=1;
OG1star=1;
OG2star=1;
OG3star=1;
tfp1=1;
tfp2=1;
tfp3=1;
tfp1star=1;
tfp2star=1;
tfp3star=1;


%More shares
sharek1=k1d/ks;
sharek2=k2d/ks;
sharek3=k3d/ks;
sharek1star=k1dstar/ksstar;
sharek2star=k2dstar/ksstar;
sharek3star=k3dstar/ksstar;
shareva1h=w*h1d/(w*h1d+r*k1d);
shareva2h=w*h2d/(w*h2d+r*k2d);
shareva3h=w*h3d/(w*h3d+r*k3d);
shareva1hstar=wstar*h1dstar/(wstar*h1dstar+rstar*k1dstar);
shareva2hstar=wstar*h2dstar/(wstar*h2dstar+rstar*k2dstar);
shareva3hstar=wstar*h3dstar/(wstar*h3dstar+rstar*k3dstar);
shareva1k=1-shareva1h;
shareva2k=1-shareva2h;
shareva3k=1-shareva3h;
shareva1kstar=1-shareva1hstar;
shareva2kstar=1-shareva2hstar;
shareva3kstar=1-shareva3hstar;

vecss=[bf            
bfstar                  
c            
cstar        
h1d          
h1dstar      
h2d          
h2dstar      
h3d          
h3dstar      
hs            
hsstar       
I            
Istar        
k            
k1d          
k1dstar      
k2d          
k2dstar      
k3d          
k3dstar      
ks           
ksstar       
kstar        
Lambda       
Lambdastar   
mc1          
mc1star      
mc2          
mc2star      
mc2va        
mc2vastar    
mc3          
mc3star      
mc3va        
mc3vastar    
OG1          
OG1star      
OG2          
OG2star      
OG3          
OG3star      
p1           
p1star       
p2           
p2star       
p3           
p3star       
pd1          
pd1star      
pd2          
pd2star      
pi1          
pi1star      
pi2          
pi2star      
pi3          
pi3star      
pm1          
pm1star      
pm2          
pm2star      
Psi1d        
Psi1dstar    
Psi1x        
Psi1xstar    
Psi2d        
Psi2dstar    
Psi2x        
Psi2xstar    
Psi3         
Psi3star     
px1          
px1star      
px2          
px2star      
Q            
r            
rrf          
rrfstar      
rstar        
tfp1         
tfp1star     
tfp2         
tfp2star     
tfp3         
tfp3star     
Upsilon1d    
Upsilon1dstar
Upsilon1x    
Upsilon1xstar
Upsilon2d    
Upsilon2dstar
Upsilon2x    
Upsilon2xstar
Upsilon3     
Upsilon3star 
VA2          
VA2star      
VA3          
VA3star      
w            
wsquig       
wsquigstar   
wstar        
WTheta       
WThetastar   
Xi           
xi1d         
xi1dstar     
xi1x         
xi1xstar     
xi2d         
xi2dstar     
xi2x         
xi2xstar     
xi3          
xi3star      
Xistar       
xiw          
xiwstar      
y1           
y1bar        
y1barstar    
y1star       
y2           
y2bar        
y2barstar    
y2star       
y3           
y3star       
yD1bar       
yD1barstar   
yD2bar       
yD2barstar   
yM1bar       
yM1barstar   
yM2bar       
yM2barstar   
yX1bar       
yX1barstar   
yX2bar       
yX2barstar   
z            
zstar];        



vecpara=[beta
delta
sigmava1
sigmaiy1
sigmadm1
sigmava2
sigmaiy2
sigmadm2
sigmavi2
sigmava3
sigmavi3
sigmaic
c1 
c2
x
xstar
pop
popstar
theta1
theta1star
omega1
omega1star
theta2
theta2star
omega2
omega2star
gamma2
gamma2star
theta3
theta3star
gamma3
gamma3star
kappac
kappah
sigmac
sigmah
sigmahh
psiw
psihab
pdotss
pdotssstar
chid1
chid1star 
chid2 
chid2star 
chix1 
chix1star
chix2 
chix2star 
chi3 
chi3star 
chib
epsilon1d 
epsilon1dstar 
epsilon2d 
epsilon2dstar 
epsilon1x 
epsilon1xstar 
epsilon2x 
epsilon2xstar 
epsilon3
epsilon3star 
epsilonk
epsilonw
sigmaz 
chik 
chiz 
thetar 
thetap1 
thetap2
thetap3
thetap4
thetay 
n
rho1
rho2
rho3
rho1star
rho2star
rho3star
sigma1
sigma2 
sigma3
sigma1star
sigma2star
sigma3star
y1
y2
y3
y1star
y2star
y3star
kappahstar];


clear r rstar w wstar mc1 mc1star pd1 pd1star px1 px1star pm1 pm1star p1 p1star mc2va mc2vastar mc2 mc2star pd2 pd2star px2 px2star pm2 pm2star p2 p2star mc3va mc3vastar mc3 mc3star p3 p3star rrf rrfstar Q wsquig wsquigstar Xi Xistar WTheta WThetastar xiw xiwstar Psi1d Upsilon1d xi1d Psi1x Upsilon1x xi1x Psi2d Upsilon2d xi2d Psi2x Upsilon2x xi2x Psi3 Upsilon3 xi3 Psi1dstar Upsilon1dstar xi1dstar Psi1xstar Upsilon1xstar xi1xstar Psi2dstar Upsilon2dstar xi2dstar Psi2xstar Upsilon2xstar xi2xstar Psi3star Upsilon3star xi3star c cstar k kstar z zstar ks ksstar I Istar y3 y3star VA3 y2bar VA3star y2barstar h3d k3d h3dstar k3dstar yD2bar yM2bar yD2barstar yM2barstar yX2bar yX2barstar y2 y2star VA2 y1bar VA2star y1barstar h2d k2d h2dstar k2dstar yD1bar yM1bar yD1barstar yM1barstar yX1bar yX1barstar y1 y1star h1d k1d h1dstar k1dstar hs hsstar b bf bfstar pi1 pi2 pi3 pi1star pi2star pi3star OG1 OG2 OG3 OG1star OG2star OG3star Lambda Lambdastar tfp1 tfp2 tfp3 tfp1star tfp2star tfp3star hssquig hssquig beta delta sigmava1 sigmaiy1 sigmadm1 sigmava2 sigmaiy2 sigmadm2 sigmavi2 sigmava3 sigmavi3 sigmaic c1 c2 xh xf xhstar xfstar theta1 theta1star omega1 omega1star theta2 theta2star omega2 omega2star gamma2 gamma2star theta3 theta3star gamma3 gamma3star kappac kappah sigmac sigmah sigmahh psiw psihab pdotss pdotssstar chid1 chid1star chid2 chid2star chix1 chix1star chix2 chix2star chi3 chi3star chib epsilon1d  epsilon1dstar  epsilon2d epsilon2dstar epsilon1x epsilon1xstar epsilon2x epsilon2xstar epsilon3 epsilon3star epsilonk epsilonw sigmaz chik chiz thetar thetap1 thetap2 thetap3 thetap4 thetay n rho1 rho2 rho3 rho1star rho2star rho3star sigma1 sigma2 sigma3 sigma1star sigma2star sigma3star a1 a1star a2 2star b2 b2star diff diffstar fx fxstar fy fystar hssquigstar i int istar kappahstar kgrid kgridstar kmax kmin mcdiff shared1 shared1star shared2 shared2star shareii2 shareii2star shareii3 shareii3star sharem1 sharem1star sharem2 sharem2star shareva2 shareva2star shareva3 shareva3star ans a2star sharek1 sharek2 sharek3 sharek1star sharek2star sharek3star shareva1h shareva2h shareva3h shareva1hstar shareva2hstar shareva3hstar shareva1k shareva2k shareva3k shareva1kstar shareva2kstar shareva3kstar adjust1 adjust2 adjust1star adjust2star d1 d1star d2 d2star d3 d3star diffk diffkstar dirrlam difflamstar j jstar lamwr lamwrstar wamwrgrid lamwrgridstar lamwrmax lamwrmin difflam lamwrgrid mcdiffgrid mcdiffmax mcdiffmin nstar pop popstar shareh1 shareh2 shareh3 shareh1star shareh2star shareh3star x xstar kappahstar

var bf bfstar c cstar h1d h1dstar h2d h2dstar h3d h3dstar hs hsstar I Istar k k1d k1dstar k2d k2dstar k3d k3dstar ks ksstar kstar Lambda Lambdastar mc1 mc1star mc2 mc2star mc2va mc2vastar mc3 mc3star mc3va mc3vastar OG1 OG1star OG2 OG2star OG3 OG3star p1 p1star p2 p2star p3 p3star pd1 pd1star pd2 pd2star pi1 pi1star pi2 pi2star pi3 pi3star pm1 pm1star pm2 pm2star Psi1d Psi1dstar Psi1x Psi1xstar Psi2d Psi2dstar Psi2x Psi2xstar Psi3 Psi3star px1 px1star px2 px2star Q r rrf rrfstar rstar tfp1 tfp1star tfp2 tfp2star tfp3 tfp3star Upsilon1d Upsilon1dstar Upsilon1x Upsilon1xstar Upsilon2d Upsilon2dstar Upsilon2x Upsilon2xstar Upsilon3 Upsilon3star VA2 VA2star VA3 VA3star w wsquig wsquigstar wstar WTheta WThetastar Xi xi1d xi1dstar xi1x xi1xstar xi2d xi2dstar xi2x xi2xstar xi3 xi3star Xistar xiw xiwstar y1 y1bar y1barstar y1star y2 y2bar y2barstar y2star y3 y3star yD1bar yD1barstar yD2bar yD2barstar yM1bar yM1barstar yM2bar yM2barstar yX1bar yX1barstar yX2bar yX2barstar z zstar;
varexo shtfp1 shtfp2 shtfp3 shtfp1star shtfp2star shtfp3star;
parameters beta delta sigmava1 sigmaiy1 sigmadm1 sigmava2 sigmaiy2 sigmadm2 sigmavi2 sigmava3 sigmavi3 sigmaic c1 c2 x xstar pop popstar theta1 theta1star omega1 omega1star theta2 theta2star omega2 omega2star gamma2 gamma2star theta3 theta3star gamma3 gamma3star kappac kappah sigmac sigmah sigmahh psiw psihab pdotss pdotssstar chid1 chid1star chid2 chid2star chix1 chix1star chix2 chix2star chi3 chi3star chib epsilon1d  epsilon1dstar  epsilon2d epsilon2dstar epsilon1x epsilon1xstar epsilon2x epsilon2xstar epsilon3 epsilon3star epsilonk epsilonw sigmaz chik chiz thetar thetap1 thetap2 thetap3 thetap4 thetay n rho1 rho2 rho3 rho1star rho2star rho3star sigma1 sigma2 sigma3 sigma1star sigma2star sigma3star y1ss y2ss y3ss y1ssstar y2ssstar y3ssstar kappahstar;
beta=vecpara(1,:);
delta=vecpara(2,:);
sigmava1=vecpara(3,:);
sigmaiy1=vecpara(4,:);
sigmadm1=vecpara(5,:);
sigmava2=vecpara(6,:);
sigmaiy2=vecpara(7,:);
sigmadm2=vecpara(8,:);
sigmavi2=vecpara(9,:);
sigmava3=vecpara(10,:);
sigmavi3=vecpara(11,:);
sigmaic=vecpara(12,:);
c1=vecpara(13,:);
c2=vecpara(14,:);
x=vecpara(15,:);
xstar=vecpara(16,:);
pop=vecpara(17,:);
popstar=vecpara(18,:);
theta1=vecpara(19,:);
theta1star=vecpara(20,:);
omega1=vecpara(21,:);
omega1star=vecpara(22,:);
theta2=vecpara(23,:);
theta2star=vecpara(24,:);
omega2=vecpara(25,:);
omega2star=vecpara(26,:);
gamma2=vecpara(27,:);
gamma2star=vecpara(28,:);
theta3=vecpara(29,:);
theta3star=vecpara(30,:);
gamma3=vecpara(31,:);
gamma3star=vecpara(32,:);
kappac=vecpara(33,:);
kappah=vecpara(34,:);
sigmac=vecpara(35,:);
sigmah=vecpara(36,:);
sigmahh=vecpara(37,:);
psiw=vecpara(38,:);
psihab=vecpara(39,:);
pdotss=vecpara(40,:);
pdotssstar=vecpara(41,:);
chid1=vecpara(42,:);
chid1star=vecpara(43,:);
chid2=vecpara(44,:);
chid2star=vecpara(45,:);
chix1=vecpara(46,:);
chix1star=vecpara(47,:);
chix2=vecpara(48,:);
chix2star=vecpara(49,:);
chi3=vecpara(50,:);
chi3star=vecpara(51,:);
chib=vecpara(52,:);
epsilon1d=vecpara(53,:);
epsilon1dstar=vecpara(54,:);
epsilon2d=vecpara(55,:);
epsilon2dstar=vecpara(56,:);
epsilon1x=vecpara(57,:);
epsilon1xstar=vecpara(58,:);
epsilon2x=vecpara(59,:);
epsilon2xstar=vecpara(60,:);
epsilon3=vecpara(61,:);
epsilon3star=vecpara(62,:);
epsilonk=vecpara(63,:);
epsilonw=vecpara(64,:);
sigmaz=vecpara(65,:);
chik=vecpara(66,:);
chiz=vecpara(67,:);
thetar=vecpara(68,:);
thetap1=vecpara(69,:);
thetap2=vecpara(70,:);
thetap3=vecpara(71,:);
thetap4=vecpara(72,:);
thetay=vecpara(73,:);
n=vecpara(74,:);
rho1=vecpara(75,:);
rho2=vecpara(76,:);
rho3=vecpara(77,:);
rho1star=vecpara(78,:);
rho2star=vecpara(79,:);
rho3star=vecpara(80,:);
sigma1=vecpara(81,:);
sigma2=vecpara(82,:);
sigma3=vecpara(83,:);
sigma1star=vecpara(84,:);
sigma2star=vecpara(85,:);
sigma3star=vecpara(86,:);
y1ss=vecpara(87,:);
y2ss=vecpara(88,:);
y3ss=vecpara(89,:);
y1ssstar=vecpara(90,:);
y2ssstar=vecpara(91,:);
y3ssstar=vecpara(92,:);
kappahstar=vecpara(93,:);

model;
ks=k1d+k2d+k3d;
ksstar=k1dstar+k2dstar+k3dstar;
w^(1-sigmahh)=(1-psiw)*(xiw(-1)*w(-1)/(1+pdotss))^(1-sigmahh)+psiw*wsquig^(1-sigmahh);
wstar^(1-sigmahh)=(1-psiw)*(xiwstar(-1)*wstar(-1)/(1+pdotssstar))^(1-sigmahh)+psiw*wsquigstar^(1-sigmahh);
mc1=(1/tfp1)*[(1-theta1)^sigmava1*r^(1-sigmava1)+theta1^sigmava1*w^(1-sigmava1)]^(1/(1-sigmava1));
mc1star=(1/tfp1star)*[(1-theta1star)^sigmava1*rstar^(1-sigmava1)+theta1star^sigmava1*wstar^(1-sigmava1)]^(1/(1-sigmava1));
pd1=Psi1d*mc1;
pd1star=Psi1dstar*mc1star;
px1=Psi1x*mc1;
px1star=Psi1xstar*mc1star;
pm1=(1+c1)*px1star*Q*(p3/p3star);
pm1star=(1+c1)*px1*(p3star/p3)/Q;
p1=((x*omega1)^sigmadm1*pd1^(1-sigmadm1)+(1-x*omega1)^sigmadm1*pm1^(1-sigmadm1))^(1/(1-sigmadm1));
p1star=((xstar*omega1star)^sigmadm1*pd1star^(1-sigmadm1)+(1-xstar*omega1star)^sigmadm1*pm1star^(1-sigmadm1))^(1/(1-sigmadm1));
mc2va=[(1-theta2)^sigmava2*r^(1-sigmava2)+theta2^sigmava2*w^(1-sigmava2)]^(1/(1-sigmava2));
mc2vastar=[(1-theta2star)^sigmava2*rstar^(1-sigmava2)+theta2star^sigmava2*wstar^(1-sigmava2)]^(1/(1-sigmava2));
mc2=[gamma2^sigmavi2*mc2va^(1-sigmavi2)+(1-gamma2)^sigmavi2*p1^(1-sigmavi2)]^(1/(1-sigmavi2));
mc2star=[gamma2star^sigmavi2*mc2vastar^(1-sigmavi2)+(1-gamma2star)^sigmavi2*p1star^(1-sigmavi2)]^(1/(1-sigmavi2));
pd2=Psi2d*mc2;
pd2star=Psi2dstar*mc2star;
px2=Psi2x*mc2;
px2star=Psi2xstar*mc2star;
pm2=(1+c2)*px2star*Q*(p3/p3star);
pm2star=(1+c2)*px2*(p3star/p3)/Q;
p2=((x*omega2)^sigmadm2*pd2^(1-sigmadm2)+(1-x*omega2)^sigmadm2*pm2^(1-sigmadm2))^(1/(1-sigmadm2));
p2star=((xstar*omega2star)^sigmadm2*pd2star^(1-sigmadm2)+(1-xstar*omega2star)^sigmadm2*pm2star^(1-sigmadm2))^(1/(1-sigmadm2));
mc3va=[(1-theta3)^sigmava3*r^(1-sigmava3)+theta3^sigmava3*w^(1-sigmava3)]^(1/(1-sigmava3));
mc3vastar=[(1-theta3star)^sigmava3*rstar^(1-sigmava3)+theta3star^sigmava3*wstar^(1-sigmava3)]^(1/(1-sigmava3));
mc3=[gamma3^sigmavi3*mc3va^(1-sigmavi3)+(1-gamma3)^sigmavi3*p2^(1-sigmavi3)]^(1/(1-sigmavi3));
mc3star=[gamma3star^sigmavi3*mc3vastar^(1-sigmavi3)+(1-gamma3star)^sigmavi3*p2star^(1-sigmavi3)]^(1/(1-sigmavi3));
p3=Psi3*mc3;
p3star=Psi3star*mc3star;
1+rrf=(1/beta)^(1-thetar)*(p3/p3(+1))*((1+rrf(-1))*(p3/p3(-1)))^(thetar)*((p3/p3(-1))^thetap1*(p3(-1)/p3(-2))^thetap2*(p3(-2)/p3(-3))^thetap3*(p3(-3)/p3(-4))^thetap4*OG1^thetay)^(1-thetar);
1+rrfstar=(1/beta)^(1-thetar)*(p3star/p3star(+1))*((1+rrfstar(-1))*(p3star/p3star(-1)))^(thetar)*((p3star/p3star(-1))^thetap1*(p3star(-1)/p3star(-2))^thetap2*(p3star(-2)/p3star(-3))^thetap3*(p3star(-3)/p3star(-4))^thetap4*OG1star^thetay)^(1-thetar);
wsquig^(1+sigmahh/sigmah)=w^(sigmahh/sigmah)*kappah^(-1/sigmah)*sigmahh/(sigmahh-1)*Xi/WTheta;
wsquigstar^(1+sigmahh/sigmah)=wstar^(sigmahh/sigmah)*kappahstar^(-1/sigmah)*sigmahh/(sigmahh-1)*Xistar/WThetastar;
Xi=hs^((sigmah+1)/sigmah)+beta*(1-psiw)*(xiw*w/(w(+1)*(1+pdotss)))^(-sigmahh*(sigmah+1)/sigmah)*Xi(+1);
Xistar=hsstar^((sigmah+1)/sigmah)+beta*(1-psiw)*(xiwstar*wstar/(wstar(+1)*(1+pdotssstar)))^(-sigmahh*(sigmah+1)/sigmah)*Xistar(+1);
WTheta=Lambda*hs/p3+beta*(1-psiw)/(1+pdotss)*xiw*(xiw*w/(w(+1)*(1+pdotss)))^(-sigmahh)*WTheta(+1);
WThetastar=Lambdastar*hsstar/p3star+beta*(1-psiw)/(1+pdotssstar)*xiwstar*(xiwstar*wstar/(wstar(+1)*(1+pdotssstar)))^(-sigmahh)*WThetastar(+1);
xiw=(1+pdotss)*(w/w(-1))^epsilonw;
xiwstar=(1+pdotssstar)*(wstar/wstar(-1))^epsilonw;
Psi1d=sigmaiy1/(sigmaiy1-1+chid1*Upsilon1d);
Upsilon1d=xi1d*(xi1d+1)-beta*(Lambda(+1)/Lambda)*(pd1(+1)/pd1)*(yD1bar(+1)/yD1bar)*xi1d(+1)*(xi1d(+1)+1);
xi1d=(pd1/pd1(-1))/((pd1(-1)/pd1(-2))^epsilon1d)-1;
Psi1x=sigmaiy1/(sigmaiy1-1+chix1*Upsilon1x);
Upsilon1x=xi1x*(xi1x+1)-beta*(Lambda(+1)/Lambda)*(px1(+1)/px1)*(yX1bar(+1)/yX1bar)*xi1x(+1)*(xi1x(+1)+1);
xi1x=(px1/px1(-1))/((px1(-1)/px1(-2))^epsilon1x)-1;
Psi2d=sigmaiy2/(sigmaiy2-1+chid2*Upsilon2d);
Upsilon2d=xi2d*(xi2d+1)-beta*(Lambda(+1)/Lambda)*(pd2(+1)/pd2)*(yD2bar(+1)/yD2bar)*xi2d(+1)*(xi2d(+1)+1);
xi2d=(pd2/pd2(-1))/((pd2(-1)/pd2(-2))^epsilon2d)-1;
Psi2x=sigmaiy2/(sigmaiy2-1+chix2*Upsilon2x);
Upsilon2x=xi2x*(xi2x+1)-beta*(Lambda(+1)/Lambda)*(px2(+1)/px2)*(yX2bar(+1)/yX2bar)*xi2x(+1)*(xi2x(+1)+1);
xi2x=(px2/px2(-1))/((px2(-1)/px2(-2))^epsilon2x)-1;
Psi3=sigmaic/(sigmaic-1+chi3*Upsilon3);
Upsilon3=xi3*(xi3+1)-beta*(Lambda(+1)/Lambda)*(p3(+1)/p3)*(y3(+1)/y3)*xi3(+1)*(xi3(+1)+1);
xi3=(p3/p3(-1))/((p3(-1)/p3(-2))^epsilon3)-1;
Psi1dstar=sigmaiy1/(sigmaiy1-1+chid1star*Upsilon1dstar);
Upsilon1dstar=xi1dstar*(xi1dstar+1)-beta*(Lambdastar(+1)/Lambdastar)*(pd1star(+1)/pd1star)*(yD1barstar(+1)/yD1barstar)*xi1dstar(+1)*(xi1dstar(+1)+1);
xi1dstar=(pd1star/pd1star(-1))/((pd1star(-1)/pd1star(-2))^epsilon1dstar)-1;
Psi1xstar=sigmaiy1/(sigmaiy1-1+chix1star*Upsilon1xstar);
Upsilon1xstar=xi1xstar*(xi1xstar+1)-beta*(Lambdastar(+1)/Lambdastar)*(px1star(+1)/px1star)*(yX1barstar(+1)/yX1barstar)*xi1xstar(+1)*(xi1xstar(+1)+1);
xi1xstar=(px1star/px1star(-1))/((px1star(-1)/px1star(-2))^epsilon1xstar)-1;
Psi2dstar=sigmaiy2/(sigmaiy2-1+chid2star*Upsilon2dstar);
Upsilon2dstar=xi2dstar*(xi2dstar+1)-beta*(Lambdastar(+1)/Lambdastar)*(pd2star(+1)/pd2star)*(yD2barstar(+1)/yD2barstar)*xi2dstar(+1)*(xi2dstar(+1)+1);
xi2dstar=(pd2star/pd2star(-1))/((pd2star(-1)/pd2star(-2))^epsilon2dstar)-1;
Psi2xstar=sigmaiy2/(sigmaiy2-1+chix2star*Upsilon2xstar);
Upsilon2xstar=xi2xstar*(xi2xstar+1)-beta*(Lambdastar(+1)/Lambdastar)*(px2star(+1)/px2star)*(yX2barstar(+1)/yX2barstar)*xi2xstar(+1)*(xi2xstar(+1)+1);
xi2xstar=(px2star/px2star(-1))/((px2star(-1)/px2star(-2))^epsilon2xstar)-1;
Psi3star=sigmaic/(sigmaic-1+chi3star*Upsilon3star);
Upsilon3star=xi3star*(xi3star+1)-beta*(Lambdastar(+1)/Lambdastar)*(p3star(+1)/p3star)*(y3star(+1)/y3star)*xi3star(+1)*(xi3star(+1)+1);
xi3star=(p3star/p3star(-1))/((p3star(-1)/p3star(-2))^epsilon3star)-1;
(kappac)^(1/sigmac)*(c/c(-1)^psihab)^(-1/sigmac)*(1/(c(-1)^psihab))=Lambda*p3;
(kappac)^(1/sigmac)*(cstar/cstar(-1)^psihab)^(-1/sigmac)*(1/(cstar(-1)^psihab))=Lambdastar*p3star;
beta*Lambda(+1)/Lambda*(r(+1)*z(+1)+p3(+1)*(1-delta-chiz*(z(+1)^(1+sigmaz)-1)/(1+sigmaz)+chik*(k(+1)-(k/k(-1))^epsilonk*k)/k*(k/k(-1))^epsilonk))=p3+p3*chik*(k-(k(-1)/k(-2))^epsilonk*k(-1))/k(-1);
beta*Lambdastar(+1)/Lambdastar*(rstar(+1)*zstar(+1)+p3star(+1)*(1-delta-chiz*(zstar(+1)^(1+sigmaz)-1)/(1+sigmaz)+chik*(kstar(+1)-(kstar/kstar(-1))^epsilonk*kstar)/kstar*(kstar/kstar(-1))^epsilonk))=p3star+p3star*chik*(kstar-(kstar(-1)/kstar(-2))^epsilonk*kstar(-1))/kstar(-1);
r/p3=chiz*z^sigmaz;
rstar/p3star=chiz*zstar^sigmaz;
ks=z*k(-1);
ksstar=zstar*kstar(-1);
I=k-(1-delta)*k(-1)+chiz*(z^(1+sigmaz)-1)*k(-1)/(1+sigmaz)+chik/2*(k-(k(-1)/k(-2))^epsilonk*k(-1))^2/k(-1);
Istar=kstar-(1-delta)*kstar(-1)+chiz*(zstar^(1+sigmaz)-1)*kstar(-1)/(1+sigmaz)+chik/2*(kstar-(kstar(-1)/kstar(-2))^epsilonk*kstar(-1))^2/kstar(-1);
y3=c+I;
y3star=cstar+Istar;
VA3=(1/tfp3)^(1-sigmavi3)*gamma3^sigmavi3*(mc3va/mc3)^(-sigmavi3)*y3;
y2bar=(1/tfp3)^(1-sigmavi3)*(1-gamma3)^sigmavi3*(p2/mc3)^(-sigmavi3)*y3;
VA3star=(1/tfp3star)^(1-sigmavi3)*gamma3star^sigmavi3*(mc3vastar/mc3star)^(-sigmavi3)*y3star;
y2barstar=(1/tfp3star)^(1-sigmavi3)*(1-gamma3star)^sigmavi3*(p2star/mc3star)^(-sigmavi3)*y3star;
h3d=theta3^sigmava3*(w/mc3va)^(-sigmava3)*VA3;
k3d=(1-theta3)^sigmava3*(r/mc3va)^(-sigmava3)*VA3;
h3dstar=theta3^sigmava3*(wstar/mc3vastar)^(-sigmava3)*VA3star;
k3dstar=(1-theta3)^sigmava3*(rstar/mc3vastar)^(-sigmava3)*VA3star;
yD2bar=(x*omega2)^sigmadm2*(pd2/p2)^(-sigmadm2)*y2bar;
yM2bar=(1-x*omega2)^sigmadm2*(pm2/p2)^(-sigmadm2)*y2bar;
yD2barstar=(xstar*omega2star)^sigmadm2*(pd2star/p2star)^(-sigmadm2)*y2barstar;
yM2barstar=(1-xstar*omega2star)^sigmadm2*(pm2star/p2star)^(-sigmadm2)*y2barstar;
yX2bar=(1+c2)*yM2barstar;
yX2barstar=(1+c2)*yM2bar;
y2=yD2bar+yX2bar;
y2star=yD2barstar+yX2barstar;
VA2=(1/tfp2)^(1-sigmavi2)*gamma2^sigmavi2*(mc2va/mc2)^(-sigmavi2)*y2;
y1bar=(1/tfp2)^(1-sigmavi2)*(1-gamma2)^sigmavi2*(p1/mc2)^(-sigmavi2)*y2;
VA2star=(1/tfp2star)^(1-sigmavi2)*gamma2star^sigmavi2*(mc2vastar/mc2star)^(-sigmavi2)*y2star;
y1barstar=(1/tfp2star)^(1-sigmavi2)*(1-gamma2star)^sigmavi2*(p1star/mc2star)^(-sigmavi2)*y2star;
h2d=theta2^sigmava2*(w/mc2va)^(-sigmava2)*VA2;
k2d=(1-theta2)^sigmava2*(r/mc2va)^(-sigmava2)*VA2;
h2dstar=theta2^sigmava2*(wstar/mc2vastar)^(-sigmava2)*VA2star;
k2dstar=(1-theta2)^sigmava2*(rstar/mc2vastar)^(-sigmava2)*VA2star;
yD1bar=(x*omega1)^sigmadm1*(pd1/p1)^(-sigmadm1)*y1bar;
yM1bar=(1-x*omega1)^sigmadm1*(pm1/p1)^(-sigmadm1)*y1bar;
yD1barstar=(xstar*omega1star)^sigmadm1*(pd1star/p1star)^(-sigmadm1)*y1barstar;
yM1barstar=(1-xstar*omega1star)^sigmadm1*(pm1star/p1star)^(-sigmadm1)*y1barstar;
yX1bar=(1+c1)*yM1barstar;
yX1barstar=(1+c1)*yM1bar;
y1=yD1bar+yX1bar;
y1star=yD1barstar+yX1barstar;
h1d=(1/tfp1)^(1-sigmava1)*theta1^sigmava1*(w/mc1)^(-sigmava1)*y1;
k1d=(1/tfp1)^(1-sigmava1)*(1-theta1)^sigmava1*(r/mc1)^(-sigmava1)*y1;
h1dstar=(1/tfp1star)^(1-sigmava1)*theta1^sigmava1*(wstar/mc1star)^(-sigmava1)*y1star;
k1dstar=(1/tfp1star)^(1-sigmava1)*(1-theta1)^sigmava1*(rstar/mc1star)^(-sigmava1)*y1star;
1/(1+rrfstar)*(p3star/p3star(+1))*(1+chib*bfstar/p3star)=beta*Lambdastar(+1)/Lambdastar;
1/(1+rrf)*(p3/p3(+1))=beta*Lambda(+1)/Lambda;
1/(1+rrfstar)*(Q(+1)/Q)*(1+chib*bf/(p3star*Q))=1/(1+rrf);
// 1/(1+rrfstar)*(Q(+1)/Q)*(1+chib*bfstar/p3star)=1/(1+rrf);
pi1=(Psi1x-1)*mc1*yX1bar+(Psi1d-1)*mc1*yD1bar;
pi2=(Psi2x-1)*mc2*yX2bar+(Psi2d-1)*mc2*yD2bar;
pi3=(Psi3-1)*mc3*y3;
pi1star=(Psi1xstar-1)*mc1star*yX1barstar+(Psi1dstar-1)*mc1star*yD1barstar;
pi2star=(Psi2xstar-1)*mc2star*yX2barstar+(Psi2dstar-1)*mc2star*yD2barstar;
pi3star=(Psi3star-1)*mc3star*y3star;
OG1=(h1d+h2d+h3d)/(pop*hs);
OG2=y3/(tfp3*y3ss);
OG3=(y1+y2+y3)/(tfp1*y1ss+tfp2*y2ss+tfp3*y3ss);
OG1star=(h1dstar+h2dstar+h3dstar)/(popstar*hsstar);
OG2star=y3star/(tfp3star*y3ssstar);
OG3star=(y1star+y2star+y3star)/(tfp1star*y1ssstar+tfp2star*y2ssstar+tfp3star*y3ssstar);
p3*c+p3*I+(p3/p3star)*bf/Q+p3*chib/2*(bf/(p3star*Q))^2=w*hs+r*ks+pi1+pi2+pi3+(1+rrfstar(-1))*(p3/p3star(-1))*bf(-1)/Q;
p3*(n/(1-n)*c+cstar/Q)+p3*(n/(1-n)*I+Istar/Q)+p3*chib/2*(n/(1-n)*(bf/(p3star*Q))^2)+p3*chib/2*(bfstar/p3star)^2/Q=(n/(1-n)*w*hs+(p3/p3star)*wstar*hsstar/Q)+(n/(1-n)*r*ks+(p3/p3star)*rstar*ksstar/Q)+(n/(1-n)*(pi1+pi2+pi3)+(p3/p3star)*(pi1star+pi2star+pi3star)/Q);
tfp1=exp(shtfp1)*tfp1(-1)^rho1;
tfp2=exp(shtfp2)*tfp2(-1)^rho2;
tfp3=exp(shtfp3)*tfp3(-1)^rho3;
tfp1star=exp(shtfp1star)*tfp1star(-1)^rho1star;
tfp2star=exp(shtfp2star)*tfp2star(-1)^rho2star;
tfp3star=exp(shtfp3star)*tfp3star(-1)^rho3star;
hs=Lambda^sigmah*w^sigmah*((sigmahh-1)/sigmahh)^sigmah*kappah;
hsstar=Lambdastar^sigmah*wstar^sigmah*((sigmahh-1)/sigmahh)^sigmah*kappahstar;
end;

initval;
bf=vecss(1,:); 
bfstar=vecss(2,:); 
c=vecss(3,:); 
cstar=vecss(4,:); 
h1d=vecss(5,:); 
h1dstar=vecss(6,:); 
h2d=vecss(7,:); 
h2dstar=vecss(8,:); 
h3d=vecss(9,:); 
h3dstar=vecss(10,:);  
hs=vecss(11,:);  
hsstar=vecss(12,:);
I=vecss(13,:); 
Istar=vecss(14,:);
k=vecss(15,:);
k1d=vecss(16,:);
k1dstar=vecss(17,:);
k2d=vecss(18,:);
k2dstar=vecss(19,:);
k3d=vecss(20,:);
k3dstar=vecss(21,:);
ks=vecss(22,:);
ksstar=vecss(23,:);
kstar=vecss(24,:);
Lambda=vecss(25,:);
Lambdastar=vecss(26,:);
mc1=vecss(27,:);
mc1star=vecss(28,:);
mc2=vecss(29,:);
mc2star=vecss(30,:);
mc2va=vecss(31,:);
mc2vastar=vecss(32,:);
mc3=vecss(33,:);
mc3star=vecss(34,:);
mc3va=vecss(35,:);
mc3vastar=vecss(36,:);
OG1=vecss(37,:);
OG1star=vecss(38,:);
OG2=vecss(39,:);
OG2star=vecss(40,:);
OG3=vecss(41,:);
OG3star=vecss(42,:);
p1=vecss(43,:);
p1star=vecss(44,:);
p2=vecss(45,:);
p2star=vecss(46,:);
p3=vecss(47,:);
p3star=vecss(48,:);
pd1=vecss(49,:);
pd1star=vecss(50,:);
pd2=vecss(51,:);
pd2star=vecss(52,:);
pi1=vecss(53,:);
pi1star=vecss(54,:);
pi2=vecss(55,:);
pi2star=vecss(56,:);
pi3=vecss(57,:);
pi3star=vecss(58,:);
pm1=vecss(59,:);
pm1star=vecss(60,:);
pm2=vecss(61,:);
pm2star=vecss(62,:);
Psi1d=vecss(63,:);
Psi1dstar=vecss(64,:);
Psi1x=vecss(65,:);
Psi1xstar=vecss(66,:);
Psi2d=vecss(67,:);
Psi2dstar=vecss(68,:);
Psi2x=vecss(69,:);
Psi2xstar=vecss(70,:);
Psi3=vecss(71,:);
Psi3star=vecss(72,:);
px1=vecss(73,:);
px1star=vecss(74,:);
px2=vecss(75,:);
px2star=vecss(76,:);
Q=vecss(77,:);
r=vecss(78,:);
rrf=vecss(79,:);
rrfstar=vecss(80,:);
rstar=vecss(81,:);
tfp1=vecss(82,:);
tfp1star=vecss(83,:);
tfp2=vecss(84,:);
tfp2star=vecss(85,:);
tfp3=vecss(86,:);
tfp3star=vecss(87,:);
Upsilon1d =vecss(88,:);
Upsilon1dstar=vecss(89,:);
Upsilon1x=vecss(90,:);
Upsilon1xstar=vecss(91,:);
Upsilon2d=vecss(92,:);
Upsilon2dstar=vecss(93,:);
Upsilon2x=vecss(94,:);
Upsilon2xstar=vecss(95,:);
Upsilon3=vecss(96,:);
Upsilon3star=vecss(97,:);
VA2=vecss(98,:);
VA2star=vecss(99,:);
VA3=vecss(100,:);
VA3star=vecss(101,:);
w=vecss(102,:);
wsquig=vecss(103,:);
wsquigstar=vecss(104,:);
wstar=vecss(105,:);
WTheta=vecss(106,:);
WThetastar=vecss(107,:);
Xi=vecss(108,:);
xi1d=vecss(109,:);
xi1dstar=vecss(110,:);
xi1x=vecss(111,:);
xi1xstar=vecss(112,:);
xi2d=vecss(113,:);
xi2dstar=vecss(114,:);
xi2x=vecss(115,:);
xi2xstar=vecss(116,:);
xi3=vecss(117,:);
xi3star=vecss(118,:);
Xistar=vecss(119,:);
xiw=vecss(120,:);
xiwstar=vecss(121,:);
y1=vecss(122,:);
y1bar=vecss(123,:);
y1barstar=vecss(124,:);
y1star=vecss(125,:);
y2=vecss(126,:);
y2bar=vecss(127,:);
y2barstar=vecss(128,:);
y2star=vecss(129,:);
y3=vecss(130,:);
y3star=vecss(131,:);
yD1bar=vecss(132,:);
yD1barstar=vecss(133,:);
yD2bar=vecss(134,:);
yD2barstar=vecss(135,:);
yM1bar=vecss(136,:);
yM1barstar=vecss(137,:);
yM2bar=vecss(138,:);
yM2barstar=vecss(139,:);
yX1bar=vecss(140,:);
yX1barstar=vecss(141,:);
yX2bar=vecss(142,:);
yX2barstar=vecss(143,:);
z=vecss(144,:);
zstar=vecss(145,:);
end;
steady;
// steady(solve_algo = 0 );

check;

shocks;
var shtfp1=sigma1^2;
var shtfp2=sigma2^2;
var shtfp3=sigma3^2;
var shtfp1star=sigma1star^2;
var shtfp2star=sigma2star^2;
var shtfp3star=sigma3star^2;
end;


// Perform long stochastic simulation
stoch_simul(order=1,irf=16) p3 p3star OG1 OG1star;