function [parameters, LLF, Ht, Ht_1, Ht_0, smooth_t, filter_t, stderrors]=MLE_RS_Multi_BEKK_GARCH_fit(data)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model: Yt=Xt*B + St*A(st)+ V(st)^1/2(ut), where V(st)=Sigma_low when st=0 and Sigma_high when st=1
%% Allow regime shifts in intercepts and cov-var matrix
% where Yt=[y1,..,ym]'and Xt=[1,x1,...,xn], St=I(m) when St=1 otherwise St=0
% Used to test the Regime-Switching Mulitrivate case, through MLE approach
% Written By Weijun Xu, School of Banking & Finance, UNSW
% First edition, 01/08/08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T,k]=size(data);
num_cov=k*(k+1)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% We need to Set some smart initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%cd ('C:\Research_Weijun_Xu\Research\Hedging_MRS\matlab\Multi_GARCH')
%[parameters_single] = diagonal_bekk_mvgarch(data,1,1); % we draw A, B from single-state Diagonal-BEKK regression

cd ('C:\Research_Weijun_Xu\Research\Hedging_MRS\matlab\Regime_Switching')
fprintf('\n %%%%%%%%%%%%%%%%%%%% Estimating Markov-Regime-Switching on Multivariate GARCH model %%%%%%%%%%%%%%%%%%%%%% \n')

%A=parameters_single(num_cov+1:num_cov+k);
%B=parameters_single(num_cov+k+1:num_cov+2*k);
%A1_draw=A; A0_draw=0.8*A;
%B1_draw=B; B0_draw=0.8*B;
A1_draw=0.2*ones(k,1);  A0_draw=0.4*ones(k,1);
B1_draw=0.98*ones(k,1); B0_draw=0.9*ones(k,1);

tp11=0.9;
tp00=0.9;
para_start=[A1_draw; A0_draw; B1_draw; B0_draw; tp11; tp00];
num_para=length(para_start);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the constraints: (1). tp11, tp00 ~(0,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LB=[-inf*ones(num_para-2,1); (1e-6)*ones(2,1)];
UB=[(1-1e-6)*ones(num_para,1)]; 
options=optimset('fmincon');
options=optimset(options,'Display','iter','Diagnostics','off','LargeScale','off','DiffMaxChange', 1e-2);

[parameters] = fmincon(@(parameters) MLE_RS_Multi_BEKK_GARCH_likelihood(parameters, data, T, k),...
   para_start,[],[],[],[],LB,UB,[],options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the Dynamic Conditioal Ht and 
%% Kim's Smooth probability: Smooth_t(t)=filter_t(t).*{[Smooth_t(t+1)./Conditional_t(t+1)]*P'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[LLF,log_likelihood, Ht, Ht_1, Ht_0, filter_t, conditional_t]= MLE_RS_Multi_BEKK_GARCH_likelihood(parameters, data, T, k);
LLF=-LLF;

P=[parameters(end-1), 1-parameters(end-1); 1-parameters(end), parameters(end)];
smooth_t=zeros(T,2);
smooth_t(end,:)=filter_t(end,:);
for i=T-1:-1:1
    smooth_t(i,:)=filter_t(i,:).*((smooth_t(i+1,:)./conditional_t(i+1,:))*P');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the standard errors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Hessian
if nargout>5
    fprintf('\n Calculating the standard errors.... \n')
    hessian=hessian_2sided('MLE_RS_Multi_BEKK_GARCH_likelihood',parameters, data, T, k);
    h=max(abs(parameters/2),1e-2)*eps^(1/3);
    hplus=parameters+h;
    hminus=parameters-h;
    likelihoodsplus=zeros(T,length(parameters));
    likelihoodsminus=zeros(T,length(parameters));
    % Calculate the Gradients
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hplus(i);
        [HOLDER, indivlike] = MLE_RS_Multi_BEKK_GARCH_likelihood(hparameters, data, T, k);
        likelihoodsplus(:,i)=indivlike;
    end
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hminus(i);
        [HOLDER, indivlike] = MLE_RS_Multi_BEKK_GARCH_likelihood(hparameters, data, T, k);
        likelihoodsminus(:,i)=indivlike;
    end
    scores=(likelihoodsplus-likelihoodsminus)./(2*repmat(h',T,1));
    B=scores'*scores;
    stderrors=inv(hessian)*B*inv(hessian);
    stderrors=sqrt(diag(stderrors));
end

% Function End



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sub function 2: {Non-Linear Constraints}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [C,Ceq]=my_con(parameters,k)

% the constraint is that the diagonal of covariance matrix 1 is greater
% than covariance matrix 0



sigma0_chol_low=parameters(num_eq*(1+k)+3:num_eq*(1+k)+num_cov+2,:);            % low tranigular matrix of low volatility cov-var matrix
sigma1_chol_low=parameters(num_eq*(1+k)+num_cov+3:num_eq*(1+k)+2*num_cov+2,:);  % low tranigular matrix of high volatility cov-var matrix

sigma0_chol_low=ivech(sigma0_chol_low);
sigma1_chol_low=ivech(sigma1_chol_low);       % transfer the column vector to low-trangular matrix
sigma0_draw=sigma0_chol_low*sigma0_chol_low'; 
sigma1_draw=sigma1_chol_low*sigma1_chol_low'; % transfer low-trangular matrix to cov-variance matrix
sigma0=diag(sigma0_draw); 
sigma1=diag(sigma1_draw);  
C=sigma0-sigma1;
Ceq=[];








function [LLF,log_likelihood, Ht, Ht_1, Ht_0, filter_t, conditional_t]= MLE_RS_Multi_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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A1=parameters(1:k,:);
A0=parameters(k+1:2*k,:);
B1=parameters(2*k+1:3*k,:);
B0=parameters(3*k+1:4*k,:);
tp11=parameters(4*k+1:4*k+1,:);         % transistion probability tp11
tp00=parameters(4*k+2: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)=(uncond_cov-diag(A1)*uncond_cov*diag(A1)-diag(B1)*uncond_cov*diag(B1))+...
        diag(A1)*data(i-1,:)'*data(i-1,:)*diag(A1)+diag(B1)*Ht(:,:,i-1)*diag(B1);
    Ht_0(:,:,i)=(uncond_cov-diag(A0)*uncond_cov*diag(A0)-diag(B0)*uncond_cov*diag(B0))+...
        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,:);





function [parameters, LLF, Ht, Ht_1, Ht_0, smooth_t, filter_t, stderrors]=MLE_RS_Multi_BEKK_GARCH_fit(data)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model: Yt=Xt*B + St*A(st)+ V(st)^1/2(ut), where V(st)=Sigma_low when st=0 and Sigma_high when st=1
%% Allow regime shifts in intercepts and cov-var matrix
% where Yt=[y1,..,ym]'and Xt=[1,x1,...,xn], St=I(m) when St=1 otherwise St=0
% Used to test the Regime-Switching Mulitrivate case, through MLE approach
% Written By Weijun Xu, School of Banking & Finance, UNSW
% First edition, 01/08/08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T,k]=size(data);
num_cov=k*(k+1)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% We need to Set some smart initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd ('C:\Research_Weijun_Xu\Research\Hedging_MRS\matlab\Multi_GARCH')
[parameters_single] = scalar_bekk_mvgarch(data,1,1); % we draw A, B from single-state Diagonal-BEKK regression

cd ('C:\Research_Weijun_Xu\Research\Hedging_MRS\matlab\Regime_Switching')
fprintf('\n %%%%%%%%%%%%%%%%%%%% Estimating Markov-Regime-Switching on Multivariate GARCH model %%%%%%%%%%%%%%%%%%%%%% \n')

C=parameters_single(1:num_cov);
A=parameters_single(num_cov+1:num_cov+1);
B=parameters_single(num_cov+2:num_cov+2);
C1_draw=1.2*C; C0_draw=C; 
A1_draw=ones(k,1)*A; A0_draw=ones(k,1)*A;
B1_draw=ones(k,1)*B; B0_draw=ones(k,1)*B;
%A1_draw=0.2*ones(k,1);  A0_draw=0.4*ones(k,1);
%B1_draw=0.98*ones(k,1); B0_draw=0.9*ones(k,1);

tp11=0.8;
tp00=0.8;
para_start=[C1_draw; C0_draw; A1_draw; A0_draw; B1_draw; B0_draw; tp11; tp00];
num_para=length(para_start);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the constraints: (1). tp11, tp00 ~(0,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LB=[-inf*ones(2*num_cov,1);(1e-6)*ones(num_para-2,1)];
UB=[inf*ones(2*num_cov,1); (1-1e-6)*ones(num_para-2*num_cov,1)]; 
options=optimset('fmincon');
options=optimset(options,'Display','iter','Diagnostics','off','LargeScale','off','DiffMaxChange', 1e-2);

[parameters] = fmincon(@(parameters)MLE_RS_Multi_BEKK_GARCH_likelihood(parameters, data, T, k, num_cov),...
   para_start,[],[],[],[],LB,UB,[],options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the Dynamic Conditioal Ht and 
%% Kim's Smooth probability: Smooth_t(t)=filter_t(t).*{[Smooth_t(t+1)./Conditional_t(t+1)]*P'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[LLF,log_likelihood, Ht, Ht_1, Ht_0, filter_t, conditional_t]= MLE_RS_Multi_BEKK_GARCH_likelihood(parameters, data, T, k, num_cov);
LLF=-LLF;

P=[parameters(end-1), 1-parameters(end-1); 1-parameters(end), parameters(end)];
smooth_t=zeros(T,2);
smooth_t(end,:)=filter_t(end,:);
for i=T-1:-1:1
    smooth_t(i,:)=filter_t(i,:).*((smooth_t(i+1,:)./conditional_t(i+1,:))*P');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the standard errors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Hessian
if nargout>7
    fprintf('\n Calculating the standard errors.... \n')
    hessian=hessian_2sided('MLE_RS_Multi_BEKK_GARCH_likelihood',parameters, data, T, k, num_cov);
    h=max(abs(parameters/2),1e-2)*eps^(1/3);
    hplus=parameters+h;
    hminus=parameters-h;
    likelihoodsplus=zeros(T,length(parameters));
    likelihoodsminus=zeros(T,length(parameters));
    % Calculate the Gradients
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hplus(i);
        [HOLDER, indivlike] = MLE_RS_Multi_BEKK_GARCH_likelihood(hparameters, data, T, k, num_cov);
        likelihoodsplus(:,i)=indivlike;
    end
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hminus(i);
        [HOLDER, indivlike] = MLE_RS_Multi_BEKK_GARCH_likelihood(hparameters, data, T, k, num_cov);
        likelihoodsminus(:,i)=indivlike;
    end
    scores=(likelihoodsplus-likelihoodsminus)./(2*repmat(h',T,1));
    B=scores'*scores;
    stderrors=inv(hessian)*B*inv(hessian);
    stderrors=sqrt(diag(stderrors));
end

% Function End



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sub function 2: {Non-Linear Constraints}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [C,Ceq]=my_con(parameters,k)

% the constraint is that the diagonal of covariance matrix 1 is greater
% than covariance matrix 0



sigma0_chol_low=parameters(num_eq*(1+k)+3:num_eq*(1+k)+num_cov+2,:);            % low tranigular matrix of low volatility cov-var matrix
sigma1_chol_low=parameters(num_eq*(1+k)+num_cov+3:num_eq*(1+k)+2*num_cov+2,:);  % low tranigular matrix of high volatility cov-var matrix

sigma0_chol_low=ivech(sigma0_chol_low);
sigma1_chol_low=ivech(sigma1_chol_low);       % transfer the column vector to low-trangular matrix
sigma0_draw=sigma0_chol_low*sigma0_chol_low'; 
sigma1_draw=sigma1_chol_low*sigma1_chol_low'; % transfer low-trangular matrix to cov-variance matrix
sigma0=diag(sigma0_draw); 
sigma1=diag(sigma1_draw);  
C=sigma0-sigma1;
Ceq=[];




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)=1.2*uncond_cov;  Ht_0(:,:,1)=0.8*uncond_cov; 
%data=[sqrt(diag(cov(data)))';data];
data=[zeros(1,k);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,:);


