function   [LLF,log_likelihood, filter_t, conditional_t]= MLE_RS_uni_all_likelihood(parameters, y, x)
% Using MLE to estimate the univariate regime switching models
% conditional_t=Pr(St|Y_t-1)
% filter_t=Pr(St|Y_t)

[T,k]=size(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assign the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta1_draw=parameters(1:k,:);              % coeff at regime 1
beta0_draw=parameters(k+1:2*k,:);          % coeff at regime 0
sigma1_draw=parameters(2*k+1:2*k+1,:);     % high volatility variance sigma1_draw
sigma0_draw=parameters(2*k+2:2*k+2,:);     % low volatility variance sigma0_draw
tp11=parameters(2*k+3:2*k+3,:);            % transistion probability tp11
tp00=parameters(2*k+4:2*k+4,:);            % transistion probability tp00

P=[tp11, 1-tp11; 1-tp00, tp00];            % transistion matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the initial settings: initial steady-probability, filter,
%% conditional probability
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P_transpose=P';
i_m=ones(2,1);
I_m=eye(2);
A=[I_m-P_transpose; i_m'];
temp_mat=((A'*A)^-1) *A';
filter_t(1,:)=temp_mat(:,end)';    % the last column of (A'A)^-1*A' is the initial steady-state probability

ystar=[mean(y);y];                  % extend the y and x
xstar=[mean(x);x];
conditional_t=zeros(T+1,2);

for i=2:T+1
    likelihood(i,1)=1/sqrt(2*pi*sigma1_draw)*...
        exp(-0.5*(ystar(i,1)-xstar(i,:)*beta1_draw).^2/(sigma1_draw));       % when st=1
    likelihood(i,2)=1/sqrt(2*pi*sigma0_draw)*...
        exp(-0.5*(ystar(i,1)-xstar(i,:)*beta0_draw).^2/(sigma0_draw));       % when st=0
    if sum(likelihood(i,:)==0)>0        % i.e. two likelihood are both zero
        likelihood(i,:)=1e-30.*ones(2,1);
    end
    
    conditional_t(i,:)=filter_t(i-1,:)*P;
    joint_likelihood(i,1)=sum(conditional_t(i,:).*likelihood(i,:));
    filter_t(i,:)=(conditional_t(i,:).*likelihood(i,:))./joint_likelihood(i,1);   % Normalized filter_t
end

log_likelihood=log(joint_likelihood(2:end,:));
LLF=-sum(log_likelihood);
filter_t=filter_t(2:end,:);
conditional_t=conditional_t(2:end,:);
