function [LLF,log_likelihood, Ht, Ht_1, Ht_0, filter_t, conditional_t]= MLE_RS_Multi_BEKK_GARCH_likelihood(parameters, data, T, k, num_cov)

% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
A1=parameters(2*num_cov+1:2*num_cov+k,:);
A0=parameters(2*num_cov+k+1:2*num_cov+2*k,:);
B1=parameters(2*num_cov+2*k+1:2*num_cov+3*k,:);
B0=parameters(2*num_cov+3*k+1:2*num_cov+4*k,:);
tp11=parameters(2*num_cov+4*k+1:2*num_cov+4*k+1,:);         % transistion probability tp11
tp00=parameters(2*num_cov+4*k+2:2*num_cov+4*k+2,:);         % 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=zeros(k,k,T+1);   Ht_1=zeros(k,k,T+1);  Ht_0=zeros(k,k,T+1);
uncond_cov=cov(data);
Ht(:,:,1)=uncond_cov; Ht_1(:,:,1)=uncond_cov;  Ht_0(:,:,1)=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_1(:,:,i)=C1*C1'+ diag(A1)*data(i-1,:)'*data(i-1,:)*diag(A1)+diag(B1)*Ht(:,:,i-1)*diag(B1);
    Ht_0(:,:,i)=C0*C0'+ diag(A0)*data(i-1,:)'*data(i-1,:)*diag(A0)+diag(B0)*Ht(:,:,i-1)*diag(B0);   
    
    % Calculate the likelihood under each regime at time t
    likelihood(i,1)=(2*pi)^(-0.5*k) *(det(Ht_1(:,:,i))).^(-0.5)*...
        exp(-0.5*data(i,:)*inv(Ht_1(:,:,i))*data(i,:)');                     % when st=1, we use Ht_1
    likelihood(i,2)=(2*pi)^(-0.5*k) *(det(Ht_0(:,:,i))).^(-0.5)*...
        exp(-0.5*data(i,:)*inv(Ht_0(:,:,i))*data(i,:)');                     % when st=0, we use Ht_0
    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(:,:,i)=conditional_t(i,1)*Ht_1(:,:,i-1)+conditional_t(i,2)*Ht_0(:,:,i-1);  % Calculate Ht weighted average on Ht_1 and Ht_0
end

log_likelihood=log(joint_likelihood(2:end,:));
LLF=-sum(log_likelihood);
Ht=Ht(:,:,2:end);
Ht_1=Ht_1(:,:,2:end);
Ht_0=Ht_0(:,:,2:end);

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