function [parameters, Ht_forecast, Ht, Rt_forecast, Rt]=MLE_RS_Multi_Scalar_BEKK_GARCH_forecast(data, forecast_step)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% when forecast step equals 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,k]=size(data);
num_cov=k*(k+1)/2;

[parameters, LLF, Ht, smooth_t]=MLE_RS_Multi_Scalar_BEKK_GARCH_fit(data);

% Calculate the Rt{1},{2},{3} represents: Rt_average, Rt_1(regime1), Rt_0(retime0)
for j=1:3
    for i=1:T
        Rt{j}(:,:,i)=inv(diag(sqrt(diag(Ht{j}(:,:,i)))))*Ht{j}(:,:,i)*inv(diag(sqrt(diag(Ht{j}(:,:,i)))));
    end
end
% Assign 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); % 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% when forecast step equals 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
conditional_t_forecast(1,:)=smooth_t(end,:)*P;  % Actually the last period of smooth_t is just the last period of filter_t

Ht_forecast{2}(:,:,1)=C1*C1'+A1*data(end,:)'*data(end,:)+B1*Ht{1}(:,:,end); % Regime 1's Last Ht term
Ht_forecast{3}(:,:,1)=C0*C0'+A0*data(end,:)'*data(end,:)+B0*Ht{1}(:,:,end); % Regime 0's Last Ht term
 
Rt_forecast{2}(:,:,1)=inv(diag(sqrt(diag(Ht_forecast{2}(:,:,1)))))*Ht_forecast{2}(:,:,1)*inv(diag(sqrt(diag(Ht_forecast{2}(:,:,1)))));
Rt_forecast{3}(:,:,1)=inv(diag(sqrt(diag(Ht_forecast{3}(:,:,1)))))*Ht_forecast{3}(:,:,1)*inv(diag(sqrt(diag(Ht_forecast{3}(:,:,1)))));

% Calculate the weighted average Ht_forecast
Ht_forecast{1}(:,:,1)=conditional_t_forecast(1,1)*Ht_forecast{2}(:,:,1)+conditional_t_forecast(1,2)*Ht_forecast{3}(:,:,1);  % Calculate Ht weighted average on Ht_1 and Ht_0
Rt_forecast{1}(:,:,1)=inv(diag(sqrt(diag(Ht_forecast{1}(:,:,1)))))*Ht_forecast{1}(:,:,1)*inv(diag(sqrt(diag(Ht_forecast{1}(:,:,1)))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% when forecast step is larger than 1
%% H(t+r)=Sum_0_r-2[C'C(A^2+B^2)^(i-2)]+ (A^2+B^2)^(r-1)o H(t+1), where o is the  hadamard  multiply symbol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if forecast_step>1
    tempHt_regime_1(:,:,1)=zeros(k,k);                   
    tempHt_regime_0(:,:,1)=zeros(k,k);               
    
    for i=2:forecast_step
        tempHt_regime_1(:,:,i)=tempHt_regime_1(:,:,i-1)+(C1*C1').*(A1^2+B1^2)^(i-2);        
        Ht_forecast{2}(:,:,i)=tempHt_regime_1(:,:,i)+((A1^2+B1^2).^(i-1)).*Ht_forecast{2}(:,:,1);
        Rt_forecast{2}(:,:,i)=inv(diag(sqrt(diag(Ht_forecast{2}(:,:,i)))))*Ht_forecast{2}(:,:,i)*inv(diag(sqrt(diag(Ht_forecast{2}(:,:,i)))));       
        
        tempHt_regime_0(:,:,i)=tempHt_regime_0(:,:,i-1)+(C0*C0').*(A0^2+B0^2)^(i-2);
        Ht_forecast{3}(:,:,i)=tempHt_regime_0(:,:,i)+((A0^2+B0^2).^(i-1)).*Ht_forecast{3}(:,:,1);
        Rt_forecast{3}(:,:,i)=inv(diag(sqrt(diag(Ht_forecast{3}(:,:,i)))))*Ht_forecast{3}(:,:,i)*inv(diag(sqrt(diag(Ht_forecast{3}(:,:,i)))));
        
        % Predict the multiple conditional probability
        conditional_t_forecast(i,:)=conditional_t_forecast(i-1,:)*P;  
        Ht_forecast{1}(:,:,i)=conditional_t_forecast(i,1)*Ht_forecast{2}(:,:,i)+conditional_t_forecast(i,2)*Ht_forecast{3}(:,:,i);  % Calculate Ht weighted average on Ht_1 and Ht_0
        Rt_forecast{1}(:,:,i)=inv(diag(sqrt(diag(Ht_forecast{1}(:,:,i)))))*Ht_forecast{1}(:,:,i)*inv(diag(sqrt(diag(Ht_forecast{1}(:,:,i)))));
    end
end



