clear all
a = 1;
b = 2;
rho = a/b;
% primo stato
s(1) = 1;
t(1) = 0;
% numero di transizioni di stato
% non c'e' tempo massimo
J = 10000;
for k = 2:J
  if (s(k-1) == 1)
% si e' nel primo stato
    w(k-1) = randexp(a);
    t(k) = t(k-1)+w(k-1);
    b0 = 0;
  else  
    w(k-1) = randexp(a+b);
    t(k) = t(k-1)+w(k-1);
    b0 = b;
  end
  s(k) = s(k-1)+2*randdisc([b0/(a+b0),a/(a+b0)])-3;
end
if (s(J) == 1)
  w(J) = randexp(a);
else  
  w(J) = randexp(a+b);
end
% stati visitati
nstati = max(s);  
iniziociclo = find(s == 1);
ncicli = length(iniziociclo)-1;
tautot = zeros(1,nstati);
tempociclo = 0;
for k = 1:ncicli
  ciclo = s(iniziociclo(k):iniziociclo(k+1)-1);
  tempi = w(iniziociclo(k):iniziociclo(k+1)-1);
  tempociclo = tempociclo+sum(tempi);
  tau = permanenza(nstati,ciclo,tempi);
  tautot = tautot+tau;
end
eta = (tautot/ncicli)/(tempociclo/ncicli);
norm(eta-(1-rho).*rho.^[0:nstati-1],1)
