function [LLF,log_likelihood, Ht, filter_t, conditional_t]= MLE_RS_Multi_Scalar_BEKK_GARCH_likelihood(parameters, data, T, k)

% This likelihood function can be used for MRS-VAR and MRS-SUR models
% Using MLE to estimate the Multivariate regime switching models, allowing
% intercept and cov-var matrix to regime shift
% conditional_t=Pr(St|Y_t-1)
% filter_t=Pr(St|Y_t)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assign the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num_cov=k*(k+1)/2; 
C1=parameters(1:num_cov,:);             C1=ivech(C1); % transfer column vector to low trangular matrix
C0=parameters(num_cov+1:2*num_cov,:);   C0=ivech(C0); % transfer column vector to low trangular matrix
A1=parameters(2*num_cov+1,:);
A0=parameters(2*num_cov+2,:);
B1=parameters(2*num_cov+3,:);
B0=parameters(2*num_cov+4,:);

tp11=parameters(2*num_cov+5,:);               % transistion probability tp11
tp00=parameters(2*num_cov+6,:);               % transistion probability tp00
P=[tp11, 1-tp11; 1-tp00, tp00];               % This is the transistion matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Let's first do some data extention, and set some initial values
%% filter_t=Pr(St|y_t); conditional_t=Pr(St|y_t-1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first observations are set to be the unconditional covariance
Ht{1}=zeros(k,k,T+1);   Ht{2}=zeros(k,k,T+1);  Ht{3}=zeros(k,k,T+1); % Ht{1}{2}{3} represents Average Ht, Ht(regime 1) and Ht(regime 0)
uncond_cov=cov(data);
Ht{1}(:,:,1)=uncond_cov; Ht{2}(:,:,1)=1.5*uncond_cov;  Ht{3}(:,:,1)=0.5*uncond_cov; 
data=[sqrt(diag(cov(data)))';data];
% initial steady-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);

for i=2:T+1
    conditional_t(i,:)=filter_t(i-1,:)*P;  
    % Calculate the conditional covariance matrix under each regime using Variance Targeting skill at time t
    Ht{2}(:,:,i)=C1*C1'+A1*data(i-1,:)'*data(i-1,:)+B1*Ht{1}(:,:,i-1); % Regime 1 Ht
    Ht{3}(:,:,i)=C0*C0'+A0*data(i-1,:)'*data(i-1,:)+B0*Ht{1}(:,:,i-1); % Regime 0 Ht
    
    % Calculate the likelihood under each regime at time t
    likelihood(i,1)=(2*pi)^(-0.5*k) *(det(Ht{2}(:,:,i))).^(-0.5)*...
        exp(-0.5*data(i,:)*inv(Ht{2}(:,:,i))*data(i,:)');                     % when st=1, we use Ht{2}
    likelihood(i,2)=(2*pi)^(-0.5*k) *(det(Ht{3}(:,:,i))).^(-0.5)*...
        exp(-0.5*data(i,:)*inv(Ht{3}(:,:,i))*data(i,:)');                     % when st=0, we use Ht{3}
 
    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
    
    % Calculate the joint likelihood by weighted average of likelihood under each regime 
    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
        
    % Collapse the conditional covariance matrix at time t for next time period use
    Ht{1}(:,:,i)=conditional_t(i,1)*Ht{2}(:,:,i)+conditional_t(i,2)*Ht{3}(:,:,i);  % Calculate Ht weighted average on Ht_1 and Ht_0
end

log_likelihood=log(joint_likelihood(2:end,:));
LLF=-sum(log_likelihood);
Ht{1}=Ht{1}(:,:,2:end);
Ht{2}=Ht{2}(:,:,2:end);
Ht{3}=Ht{3}(:,:,2:end);

filter_t=filter_t(2:end,:);
conditional_t=conditional_t(2:end,:);

