gg1_v3

rng(42)

%% monte carlo 
N = 1000000;
t = exprnd(2, 1, N);
s = sum(exprnd(0.5, 2, N));  % Erlang-2 distributed r.v.

if (mean(t) <= mean(s))
    error('system is unstable');
end

wq = zeros(1, N + 1);
N_on_arriv = zeros(1, N);

for i = 1:N
    wq(i + 1) = max(0, wq(i) + s(i) - t(i));    
    if (wq(i + 1) ~= 0)
        prev_idx = i;
        tot_servt = s(prev_idx);
        while (prev_idx > 0 && wq(i + 1) > tot_servt)
            prev_idx = prev_idx - 1;
            tot_servt = tot_servt + s(prev_idx);
        end
        N_on_arriv(i + 1) = i - prev_idx  + 1;
    end
end

w = wq(1:end - 1) + s;