function p=probs(y,x,alpha,theta)
% calculates p(y(3:T)|y(1), y(2),x,alpha)

    T=length(y);
    p=1;
    for t=3:T
        aa=exp(alpha+y(t-1)*theta(1)+y(t-2)*theta(2)+x(t,:)*theta(3:length(theta)));
        if y(t)==1
            p=p*aa/(1+aa);
        else
            p=p/(1+aa);
        end
    end
end
