function [LLF,log_likelihood, Ht, Rt,Q, H, filter_t, conditional_t]= MLE_RS_Multi_DCC_likelihood(parameters, data)

%% likelihood function of 'MLE_RS_Multi_DCC_fit' function.
%% Unlike single-state DCC, in the MRS-DCC framework, one-step estimation is used

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assign the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,k]=size(data);
for i=1:k
    garch{i}=parameters((i-1)*6+1:i*6,:);         % Each series's GARCH parameters including two regime coefficients
end
dcc_a1=parameters(6*k+1,:);                       % DCC correlation parameters
dcc_b1=parameters(6*k+2,:);
dcc_a0=parameters(6*k+3,:);
dcc_b0=parameters(6*k+4,:);
tp11=parameters(6*k+5,:);                         % transistion probability tp11
tp00=parameters(6*k+6,:);                         % transistion probability tp00

P=[tp11, 1-tp11; 1-tp00, tp00];                   % This is the transistion matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Let's set some reasonable initial value at time 0, to make sure whole conditional varaiance be positive definite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extend individual variance series
ht=zeros(T+1,k);   ht_1=zeros(T+1,k);  ht_0=zeros(T+1,k);
ht_bar=var(data);
ht(1,:)=ht_bar; ht_1(1,:)=1.2*ht_bar;  ht_0(1,:)=0.8*ht_bar;     % regime 1 is high volatility regime, regime 0 is low volatility regime
data=[sqrt(diag(cov(data)))';data];
stdresid=zeros(T+1,k);
stdresid(1,:)=data(1,:)./sqrt(ht_bar);
    
% Extend DCC correlation series
Qt=zeros(k,k,T+1);   Qt_1=zeros(k,k,T+1);  Qt_0=zeros(k,k,T+1);
Qbar=corr(data);
Qt(:,:,1)=Qbar; Qt_1(:,:,1)=Qbar;  Qt_0(:,:,1)=Qbar; 

% 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;  
    for j=1:k        
        % Calculate the conditional variance of each series under each reigme at time t
        ht_1(i,j)=garch{j}(1)+garch{j}(2)*data(i-1,j)^2+garch{j}(3)*ht(i-1,j);
        ht_0(i,j)=garch{j}(4)+garch{j}(5)*data(i-1,j)^2+garch{j}(6)*ht(i-1,j);
        % Collapse the conditional Variance of each series at time t for next time period use
        ht(i,j)=conditional_t(i,1)*ht_1(i,j)+conditional_t(i,2)*ht_0(i,j);                  % Calculate ht weighted average on ht_1 and ht_0        
        % Calculate the standardized residuals, under each regime 
        stdresid(i,j)=data(i,j)./sqrt(ht(i,j)); 
        stdresid_1(i,j)=data(i,j)./sqrt(ht_1(i,j)); 
        stdresid_0(i,j)=data(i,j)./sqrt(ht_0(i,j)); 
    end
    
    % Calculate the conditional covariance matrix under each regime using Variance Targeting skill at time t
    Qt_1(:,:,i)=Qbar*(1-dcc_a1-dcc_b1)+dcc_a1*stdresid(i-1,:)'*stdresid(i-1,:)+dcc_b1*Qt(:,:,i-1);
    Qt_0(:,:,i)=Qbar*(1-dcc_a0-dcc_b0)+dcc_a0*stdresid(i-1,:)'*stdresid(i-1,:)+dcc_b0*Qt(:,:,i-1);        
    Rt_1(:,:,i)=Qt_1(:,:,i)./(sqrt(diag(Qt_1(:,:,i)))*sqrt(diag(Qt_1(:,:,i)))');    
    Rt_0(:,:,i)=Qt_0(:,:,i)./(sqrt(diag(Qt_0(:,:,i)))*sqrt(diag(Qt_0(:,:,i)))');

    % Calculate the likelihood under each regime at time t    
    % One step DCC Log_Likelihood=-.5*{k*log(2pi)+2*log(|Dt|)+log(|Rt|)+et'*inv(Rt)*et}    
    Dt(:,:,i)=diag(sqrt(ht(i,:)));     % Dt=diag{sqrt(ht)}
    Dt_1(:,:,i)=diag(sqrt(ht_1(i,:))); 
    Dt_0(:,:,i)=diag(sqrt(ht_0(i,:))); 
        
    likelihood(i,1)=(2*pi)^(-0.5*k) *(det(Dt_1(:,:,i))).^(-0.5*2)*(det(Rt_1(:,:,i))).^(-0.5)*...
        exp(-0.5*stdresid_1(i,:)*inv(Rt_1(:,:,i))*stdresid_1(i,:)');                     % when st=1, we use Rt_1, Dt_1, stdresid_1
    likelihood(i,2)=(2*pi)^(-0.5*k) *(det(Dt_0(:,:,i))).^(-0.5*2)*(det(Rt_0(:,:,i))).^(-0.5)*...
        exp(-0.5*stdresid_0(i,:)*inv(Rt_0(:,:,i))*stdresid_0(i,:)');                     % when st=0, we use Rt_0, Dt_0, stdresid_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 weigQted 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 Qt at time t for next time period use
    Qt(:,:,i)=conditional_t(i,1)*Qt_1(:,:,i)+conditional_t(i,2)*Qt_0(:,:,i);      % Calculate Qt weigQted average on Qt_1 and Qt_0

    % Calculate the conditional correlation matrix Rt, using filter state probability
    Rt_average(:,:,i)=filter_t(i,1)*Rt_1(:,:,i)+filter_t(i,2)*Rt_0(:,:,i);        % Calculate Rt weigQted average on Rt_1 and Rt_0    
      
    % Calculate the conditional correlation matrix Rt, using filter state probability
    Ht_average(:,:,i)=Dt(:,:,i)*Rt_average(:,:,i)*Dt(:,:,i);
    Ht_1(:,:,i)=Dt_1(:,:,i)*Rt_1(:,:,i)*Dt_1(:,:,i);
    Ht_0(:,:,i)=Dt_0(:,:,i)*Rt_0(:,:,i)*Dt_0(:,:,i);   
end

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

Rt{1}=Rt_average(:,:,2:end); % weighted average regime 1 & 0
Rt{2}=Rt_1(:,:,2:end);       % High volatility regime 1
Rt{3}=Rt_0(:,:,2:end);       % Low volatility regime 0

Ht{1}=Ht_average(:,:,2:end); % weighted average regime 1 & 0
Ht{2}=Ht_1(:,:,2:end);       % High volatility regime 1
Ht{3}=Ht_0(:,:,2:end);       % Low volatility regime 0

H{1}=ht(2:end,:);         % weighted average regime 1 & 0
H{2}=ht_1(2:end,:);       % High volatility regime 1
H{3}=ht_0(2:end,:);       % Low volatility regime 0

Q{1}=Qt(:,:,2:end);         % weighted average regime 1 & 0
Q{2}=Qt_1(:,:,2:end);       % High volatility regime 1
Q{3}=Qt_0(:,:,2:end);       % Low volatility regime 0

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


