
function [LLF,log_likelihood, ht, filter_t,conditional_t]=MLE_RS_Uni_GARCH_likelihood(parameters, error)
%% likelihood function of 'MLE_RS_Uni_GARCH_fit' function.

T=length(error);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assign the parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta1_draw=parameters(1:3,:);              % coeff at regime 1
beta0_draw=parameters(4:2*3,:);            % coeff at regime 0
tp11=parameters(7:7,:);                    % transistion probability tp11
tp00=parameters(8:8,:);                    % transistion probability tp00

P=[tp11, 1-tp11; 1-tp00, tp00];            % transistion matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the initial settings: initial steady-probability, filter,
%% conditional probability
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% first observations are set to be the unconditional variance

ht=zeros(T+1,1);   ht_1=zeros(T+1,1);  ht_0=zeros(T+1,1);
uncond_cov=cov(error);
ht(1)=uncond_cov; ht_1(1)=uncond_cov;  ht_0(1)=uncond_cov; 
error=[sqrt(cov(error));error];

% 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 variance under each regime at time t
    ht_1(i)=beta1_draw(1)+beta1_draw(2)*error(i-1,:)^2+beta1_draw(3)*ht(i-1);
    ht_0(i)=beta0_draw(1)+beta0_draw(2)*error(i-1,:)^2+beta0_draw(3)*ht(i-1);
    
    % Calculate the likelihood under each regime at time t
    likelihood(i,1)=1/sqrt((2*pi*ht_1(i)))*exp((-0.5)*error(i)^2 /ht_1(i));       % when st=1, we use ht_1
    likelihood(i,2)=1/sqrt((2*pi*ht_0(i)))*exp((-0.5)*error(i)^2 /ht_0(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)+conditional_t(i,2)*ht_0(i);                  % 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,:);
ht=[ht,ht_1,ht_0];      % combine ht, ht_1, ht_0 in the same matrix

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


