function y=robustm(theta, seed);

format long g;

% Initiation
beta=.97;
T=50; % Length of simulation periods
P = [0.8 0.2   0; 
       0 0.5 0.5; 
       0   0   1]; % Transition matrix
n=size(P,1); 

% Rows correspond to g, d, b, s; colums to states
state=[ .5  .5 .25;  % Gov expenditure
         0   0   0;  % Endowment
       2.2 2.2 2.2;  % Pref shock
         0   0   0]; % Coupon payment
V_lag=zeros(n,1);
U=zeros(n,1);
tol1=1; tol2=1;
mu=.25; % Guess

while tol1 > (.1)^12;
   H=.5/sqrt(2)*[-1 1+mu -1-mu mu; 1 mu-1 1-mu mu];
   for k=1:n;
      U(k,1)=-state(:,k)'*H'*H*state(:,k);
   end;
         
   while tol2 > (.1)^12;
      V_lm=V_lag/(-theta); 
      V_lm2=exp(V_lm);
      V=U -beta*theta*log(P*V_lm2);
      tol2=max(abs(V-V_lag));
      V_lag=V;
   end;
   
   w=exp(-V/theta);
   Pw=zeros(n);
   for i=1:n;
      for j=1:n;       
         Pw(i,j)= w(j,1)*P(i,j)/(P(i,:)*w);
      end;
   end;
   
   % Let's assume that the initial state is #1
   h1=zeros(n,1); h2=zeros(n,1);
   for k=1:n;   
      h1(k,1)=([0 -1 1 -1]*state(:,k))^2;
      h2(k,1)=[1 -1 1 0]*state(:,k)*[1 0 0 1]*state(:,k);
   end;
   invP=inv(eye(n)-beta*Pw);
   ax=.5*[1 0 0]*invP*h1; % Note that initial state is #1
   bx=.5*[1 0 0]*invP*h2;
   
   mu_new=mu^2 +bx/ax;
   tol1=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 
rand('state', seed);
W=rand(T,1); % Uniform distribution
S=zeros(T+1, 1); S(1,1)=1;

% Preliminaries for computing government Debt
for i=1:3
  cc(i,1)=.5*[-1 1+mu 1-mu mu]*state(:,i);
  ll(i,1)=.5*[1 -1+mu 1-mu mu]*state(:,i);  
end;

% The following is applicable only to this case. Imitating the transition matrix.
for i=1:T;
   if S(i,1)==1 & W(i,1) < P(1,1);
      S(i+1)=1;
   elseif S(i,1)==1 & W(i,1) >=P(1,1);
      S(i+1)=2;
   elseif S(i,1)==2 & W(i,1)  < P(2,2);
      S(i+1)=2;
   elseif S(i,1)==2 & W(i,1) >=P(2,2);
      S(i+1)=3;
   elseif S(i,1)==3;
      S(i+1)=3;
   end;
con(i,1)=.5*[-1 1+mu 1-mu mu]*state(:,S(i,1));
lab(i,1)=.5*[1 -1+mu 1-mu mu]*state(:,S(i,1));
tax(i,1)=(2*mu*[0 -1 1 -1]*state(:,S(i,1)))/([1 -1-mu 1+mu -mu]*state(:,S(i,1)));
gov(i,1)=state(1, S(i,1));

% Govt debt
B(i,1)=invP(S(i,1),:)*((state(3,:)'-cc).*(ll-state(1,:)')-ll.^2)/(state(3,S(i,1))-con(i,1));

% One-period risk-free interest rate
bc=state(3,S(i,1))-con(i,1); % b_t-c_t
Rate(i,1)=(beta*(Pw(S(i,1),:)*(state(3,:)'-cc))/bc)^(-1)-1;
end;

y=[con lab tax gov Rate B];
