function [tau,R,y,h,c,pai,g,z,TM]=simu(state0,MU)
%[tau,R,y,h,c,pai,g,z,TM]=simu(state0,MU) computes Ramsey allocations for a given state of the economy (state0) and a given markup (MU). It also provides the processes for the exogenous gov't purchases (g) and technology (z) shocks and the state transition matrix of these shocks, TM. Note that for tau, R, y, h, c, g, and z, the stochastic process is defined by 4 numbers, one corresponding to each possible state of the economy. The process pai is given by a 4x4 matrix because both the previous and current states are relevant in determining this variable. 
%Inputs: state0 and MU (respectively, the state of the economy and the markup)
%Output: Ramsy equilibrium processes for: tau (labor income tax rate),R (nominal interest rate), y (output),h (hours), c (consumption), pai (inflation). The exogenous processes g (vog't purchases) and z (technology shock), and the state transition matrixTM. 
%Calls: L.m, param.m, findvh (via fsolve) 
%(c) Stephanie Schmitt-Grohe and Martin Uribe

LAMBDA = L(state0,MU);

[BETA,THETA,D,ETA,B,A,gl,gh,phig,zl,zh,phiz,ALFA,TM]=param(MU);

% STATES
%(gh,zh)=1
%(gh,zl)=2;
%(gl,zh)=3;
%(gl,zl)=4;

g=[gh gh gl gl];
z=[zh zl zh zl];

g0 = g(state0);
z0 = z(state0);

OPTIONS = OPTIMSET('tolx',1e-15,'tolfun',1e-15);

Vs = sqrt(B/A);
Vmax = sqrt((1+B)/A);

P0=0; %t>0$
for j=1:4
vh = fsolve('findvh',[2.7 .2],OPTIONS,LAMBDA,g(j),z(j),MU,P0);
v(j) = Vs + (Vmax-Vs) * abs(vh(1)) / (1+abs(vh(1))); %Vs + abs(vh(1)-Vs);
h(j) = abs(vh(2)) / (1+abs(vh(2)));
s(j)=A*v(j)+B/v(j)-2*sqrt(A*B);
c(j) = (z(j)*h(j)-g(j)) / (1+ALFA*s(j));
uc(j)=1/c(j);
uh(j)=-THETA/(1-h(j));
sp(j)=A-B/v(j)^2;
R(j)=1/(1-v(j)^2*sp(j)); %nominal interest rate
gama(j)=1+s(j)+v(j)*sp(j);
phi(j) = (1+(1-1/R(j))/v(j))/gama(j) +ALFA*s(j)/gama(j);
w(j) = (1+ETA)/ETA * z(j); %wage rate
tau(j) = 1 + uh(j)/uc(j) * gama(j) / w(j); %labor tax rate
y(j) = z(j)*h(j); %output
ps(j) = tau(j)*w(j)*h(j)-g(j); %primary surplus
lhs(j,1) = uc(j)*c(j)*phi(j) + uh(j)*h(j) + uc(j)*z(j)*h(j)/gama(j)/ETA;
end

P0=1; %t=0$
vh0 = fsolve('findvh',[2.7 .2],OPTIONS,LAMBDA,g0,z0,MU,P0);
v0 = Vs + (Vmax-Vs) * abs(vh0(1)) / (1+abs(vh0(1))); %Vs + abs(vh(1)-Vs);
h0 = abs(vh0(2)) / (1+abs(vh0(2)));
s0=A*v0+B/v0-2*sqrt(A*B);
c0 = (z0*h0-g0) / (1+ALFA*s0);
uc0=1/c0;
ucc0=-1/c0^2;
uch0=0;
uhc0=0;
uh0=-THETA/(1-h0);
uhh0 = -THETA/(1-h0)^2;
sp0=A-B/v0^2;
spp0=2*B/v0^3;
R0=1/(1-v0^2*sp0);
gama0=1+s0+v0*sp0;
gamap0=2*A;
phi0 = (1+(1-1/R0)/v0)/gama0 + ALFA * s0/gama0;
y0=z0*h0;%Output in period 0
w0 = (1+ETA)/ETA * z0; %wage rate
tau0 = 1 + uh0/uc0 * gama0 / w0; %labor tax rate

for pr=1:length(g)
ELHS(pr) = lhs(pr) + BETA * TM(pr,:) * inv(eye(size(TM))-BETA*TM) * lhs;
d(pr) = gama(pr)/uc(pr)*ELHS(pr);
end

for pr=1:length(g)
for pa=1:length(g)
pai(pa,pr) = R(pa) / d(pr) * (d(pa) + g(pa) - tau(pa)*w(pa)*h(pa) - c(pa)/v(pa) * (1-1/R(pa))); 
end
end