function y=robusts(theta, seed);

format long;

% Initiation ... A, C, and S's
rho=.95;
mg=.35;
A=[0 0 0 rho mg*(1-rho); 
   1 0 0 0 0 ;
   0 1 0 0 0 ;
   0 0 1 0 0;
   0 0 0 0 1];
C=[mg/8*sqrt(1-rho^2);0;0;0;0];
Sg = [1 0 0 0 0];        
Sd = [0 0 0 0 0];       
Sb = [0 0 0 0 2.135];
Ss = [0 0 0 0 0];    

beta=.97;
T=50; % Length of simulation periods

% Initial condition
x0=[.35;.35;.35;.35;1];

m=size(C, 1); n=size(C, 2);
alpha = -(theta)^(-1);
V_lag=zeros(m,m);
tol1=1; tol2=1;
M1=(Sb-Sd-Ss)'*(Sb-Sd-Ss);
M2=(Sb-Sd-Sg)'*(Sg+Ss);
mu=.25; % Guess

while tol2 > (.1)^12;
   H=1/2/sqrt(2)*[-1 1+mu -1-mu mu; 1 mu-1 1-mu mu]*[Sg; Sd; Sb; Ss];
   
   while tol1 > (.1)^12;
      D_V = V_lag + alpha*V_lag*C*inv(eye(n) - alpha*C'*V_lag*C)*C'*V_lag;
      V=-H'*H + beta*A'*D_V*A;
      tol1=max(max(abs(V-V_lag)));
      V_lag=V;
   end;
   
   O=alpha*inv(eye(n)-alpha*C'*V*C)*C'*V;
   Ah = A + C*O;

   Q1=dlyap(sqrt(beta)*Ah', M1);
   q1=beta/(1-beta)*trace(C'*Q1*C);
   ax=x0'*Q1*x0+q1;
   Q2=dlyap(sqrt(beta)*Ah', M2);
   q2=beta/(1-beta)*trace(C'*Q2*C);
   bx=x0'*Q2*x0+q2;

   mu_new=mu^2 + bx/ax;
   tol2=abs(mu_new-mu);
   mu=mu_new;

   if mu <= 0 | mu >= .5;
      break; disp('The value of mu is contradictory');
   end;
end;

lambda = mu/(1-2*mu);

% Simulation
randn('state', seed);
W=randn(n,T); % Normal with variance 1
X=zeros(m, T+1); X(:,1)=x0;

% Preliminaries for computing the government debt
M3=[0 0.5 -0.5; 0.5 -1 0; -0.5 0 0];
bc=[1 -1-mu 1+mu -mu]*[Sg; Sd; Sb; Ss]; % b_t-c_t
l=[1 mu-1 1-mu mu]*[Sg; Sd; Sb; Ss]; %l_t
aa=[bc; l ; Sg];
MM=aa'*M3*aa;

% Solving for B_t
Q3=dlyap(sqrt(beta)*Ah', MM);
q3=beta/(1-beta)*trace(C'*Q3*C);

% Stacking
for i=1:T;
X(:,i+1)=A*X(:,i)+C*W(:,i); % Approximated law of motion
X1(:,i+1)=Ah*X(:,i); % Agents' conditional expectation
con(i,1)=.5*(Sb+Sd-Sg -mu*(Sb-Sd-Ss))*X(:,i); % Consumption
lab(i,1)=.5*(Sb-Sd+Sg -mu*(Sb-Sd-Ss))*X(:,i); % Labor
tax(i,1)=(2*mu*(Sb-Sd-Ss)*X(:,i))/((Sb-Sd+Sg +mu*(Sb-Sd-Ss))*X(:,i)); % Tax
gov(i,1)=Sg*X(:,i); % Government expenditure
b(i,1)=Sb*X(:,i);
p(i,1)=1/mu*(b(i,1)-con(i,1)); % Price
con1(i,1)=.5*(Sb+Sd-Sg -mu*(Sb-Sd-Ss))*X1(:,i+1); % Conditional expectation
b1(i,1)=Sb*X1(:,i+1); % Conditional expectation
rate(i,1)=(beta*(b1(i,1)-con1(i,1))/(b(i,1)-con(i,1)))^(-1) -1; % One-period risk-free net rate
BX(i,1)=(X(:,i)'*Q3*X(:,i)+q3)/(b(i,1)-con(i,1));%government Debt
end;

% Output
y=[con lab tax gov rate BX];
