function [LLF,log_likelihood, filter_t, conditional_t]= MLE_RS_Multi_all_likelihood(parameters, yt, xt, T, num_eq, num_cov, k)

%% Design for 'MLE_RS_Multi_all_fit' & 'MLE_RS_Multi_all_SUR_fit' function use.
%% Using MLE to estimate the Multivariate regime switching models, allowing all coefficients and cov-var matrix to regime shift
% Conditional_t=Pr(St|Y_t-1)
% filter_t=Pr(St|Y_t)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assign the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta1_draw=parameters(1:num_eq*k,:);                                      % Coefficients in regime 1
beta0_draw=parameters(num_eq*k+1:2*num_eq*k,:);                           % Coefficients in regime 2
sigma1_chol=parameters(2*num_eq*k+1:2*num_eq*k+num_cov,:);                % low tranigular matrix of high volatility cov-var matrix
sigma0_chol=parameters(2*num_eq*k+num_cov+1:2*num_eq*k+2*num_cov,:);      % low tranigular matrix of low volatility cov-var matrix
tp11=parameters(2*num_eq*k+2*num_cov+1:2*num_eq*k+2*num_cov+1,:);         % transistion probability tp11
tp00=parameters(2*num_eq*k+2*num_cov+2:2*num_eq*k+2*num_cov+2,:);         % transistion probability tp00

% construct the cov-var matrix
sigma1_chol=ivech(sigma1_chol);               % transfer the column vector to low-trangular matrix
sigma0_chol=ivech(sigma0_chol);
sigma1_draw=sigma1_chol*sigma1_chol';         % transfer low-trangular matrix to cov-variance matrix
sigma0_draw=sigma0_chol*sigma0_chol'; 
P=[tp11, 1-tp11; 1-tp00, tp00];               % This is the 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
conditional_t=zeros(T+1,2);

ystar=[mean(yt,2), yt];                  % extend the yt and xt
xstar=zeros(num_eq, k*num_eq, T+1);
xstar(:,:,1)=mean(xt,3);
xstar(:,:,2:end)=xt;

for i=2:T+1
    likelihood(i,1)=(2*pi)^(-0.5*num_eq) *(det(sigma1_draw)).^(-0.5)*...
        exp(-0.5*(ystar(:,i)-xstar(:,:,i)*beta1_draw)'*inv(sigma1_draw)*...
        (ystar(:,i)-xstar(:,:,i)*beta1_draw));      % when st=1, we use sigma1_draw, and beta1_draw
    likelihood(i,2)=(2*pi)^(-0.5*num_eq)*(det(sigma0_draw)).^(-0.5)*...
        exp(-0.5*(ystar(:,i)-xstar(:,:,i)*beta0_draw)'*inv(sigma0_draw)*...
        (ystar(:,i)-xstar(:,:,i)*beta0_draw));      % when st=0, we use sigma0_draw, and beta0_draw
    
    if sum(likelihood(i,:)==0)>0        % i.e. two likelihood are both zero, avoid 0/0 problem
        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,:);
