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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% when forecast step equals 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,k]=size(data);

[parameters, LLF, Ht, Rt, Qt, H, smooth_t]=MLE_RS_Multi_DCC_fit(data);
% Note that Ht, Rt, Qt, H includes {1}{2}{3}sub_components, represents
% average, high volatility regime, low volatility regime.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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
% First let's get the standardized residuals
for i=1:T
    for j=1:k
        stdresid(i,j)=data(i,j)./sqrt(H{1}(i,j));
        stdresid_1(i,j)=data(i,j)./sqrt(H{2}(i,j));
        stdresid_0(i,j)=data(i,j)./sqrt(H{3}(i,j));
    end
end

% Calculate one-step forecast of each time series's conditional variance
for j=1:k
    % Calculate the conditional variance of each series under each reigme at time t
    H_forecast{2}(1,j)=garch{j}(1)+garch{j}(2)*data(end,j)^2+garch{j}(3)*H{1}(end,j);
    H_forecast{3}(1,j)=garch{j}(4)+garch{j}(5)*data(end,j)^2+garch{j}(6)*H{1}(end,j);
    % Collapse the conditional Variance of each series at time t for next time period use
    H_forecast{1}(1,j)=conditional_t_forecast(1,1)*H_forecast{2}(1,j)+conditional_t_forecast(1,2)*H_forecast{3}(1,j);                  % Calculate ht weighted average on ht_1 and ht_0
end

Qbar=corr(data);
Qt_forecast{2}(:,:,1)=Qbar*(1-dcc_a1-dcc_b1)+dcc_a1*(stdresid(end,:)'*stdresid(end,:))+ dcc_b1*Qt{1}(:,:,end);
Qt_forecast{3}(:,:,1)=Qbar*(1-dcc_a0-dcc_b0)+dcc_a0*(stdresid(end,:)'*stdresid(end,:))+ dcc_b0*Qt{1}(:,:,end);
Qt_forecast{1}(:,:,1)=conditional_t_forecast(1,1)*Qt_forecast{2}(:,:,1)+conditional_t_forecast(1,2)*Qt_forecast{3}(:,:,1);      % Calculate Qt weigQted average on Qt_1 and Qt_0

Rt_forecast{2}(:,:,1)=Qt_forecast{2}(:,:,1)./(sqrt(diag(Qt_forecast{2}(:,:,1)))*sqrt(diag(Qt_forecast{2}(:,:,1)))');
Rt_forecast{3}(:,:,1)=Qt_forecast{3}(:,:,1)./(sqrt(diag(Qt_forecast{3}(:,:,1)))*sqrt(diag(Qt_forecast{3}(:,:,1)))');

Ht_forecast{2}(:,:,1)=diag(H_forecast{2}(1,:))*Rt_forecast{2}(:,:,1)*diag(H_forecast{2}(1,:));
Ht_forecast{3}(:,:,1)=diag(H_forecast{3}(1,:))*Rt_forecast{3}(:,:,1)*diag(H_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
%% Q(t+r)=Sum_0_r-2[(1-A-B)*(A^2+B^2)^(i)]+ (A^2+B^2)^(r-1)o Q(t+1), where o is the  hadamard  multiply symbol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if forecast_step>1
    tempQt_regime_1(:,:,1)=zeros(k,k);
    tempQt_regime_0(:,:,1)=zeros(k,k);

    tempH_regime_1=zeros(1,k);
    tempH_regime_0=zeros(1,k);  
    
    for h=2:forecast_step
        conditional_t_forecast(h,:)=conditional_t_forecast(h-1,:)*P;
        for j=1:k
            % Calculate the conditional variance of each series under each reigme at time t
            tempH_regime_1(h,j)=tempH_regime_1(h-1,j)+garch{j}(1)*(garch{j}(2)+garch{j}(3))^(h-2);
            tempH_regime_0(h,j)=tempH_regime_0(h-1,j)+garch{j}(4)*(garch{j}(5)+garch{j}(6))^(h-2);

            H_forecast{2}(h,j)=tempH_regime_1(h,j)+(garch{j}(2)+garch{j}(3))^(h-1) * H_forecast{1}(1,j);  % Regime 1 variance
            H_forecast{3}(h,j)=tempH_regime_0(h,j)+(garch{j}(4)+garch{j}(5))^(h-1) * H_forecast{1}(1,j);  % Regime 0 variance
            % Collapse the conditional Variance of each series at time t for next time period use
            H_forecast{1}(h,j)=conditional_t_forecast(h,1)*H_forecast{2}(h,j)+conditional_t_forecast(h,2)*H_forecast{3}(h,j);                  % Calculate ht weighted average on ht_1 and ht_0
        end

        tempQt_regime_1(:,:,h)=tempQt_regime_1(:,:,h-1)+Qbar*(1-dcc_a1-dcc_b1).*(dcc_a1+dcc_a1)^(h-2);
        tempQt_regime_0(:,:,h)=tempQt_regime_0(:,:,h-1)+Qbar*(1-dcc_a0-dcc_b0).*(dcc_a0+dcc_a0)^(h-2);

        Qt_forecast{2}(:,:,h)=tempQt_regime_1(:,:,h)+ (dcc_a1+dcc_a1)^(h-1) .* Qt_forecast{1}(:,:,1);     % Regime 1 covariance
        Qt_forecast{3}(:,:,h)=tempQt_regime_0(:,:,h)+ (dcc_a0+dcc_a0)^(h-1) .* Qt_forecast{1}(:,:,1);     % Regime 1 covariance

        Rt_forecast{2}(:,:,h)=Qt_forecast{2}(:,:,h)./(sqrt(diag(Qt_forecast{2}(:,:,h)))*sqrt(diag(Qt_forecast{2}(:,:,h)))');
        Rt_forecast{3}(:,:,h)=Qt_forecast{3}(:,:,h)./(sqrt(diag(Qt_forecast{3}(:,:,h)))*sqrt(diag(Qt_forecast{3}(:,:,h)))');
            
        % Predict the multiple conditional probability
        Rt_forecast{1}(:,:,h)=conditional_t_forecast(h,1)*Rt_forecast{2}(:,:,h)+conditional_t_forecast(h,2)*Rt_forecast{3}(:,:,h);  % Calculate Ht weighted average on Ht_1 and Ht_0
        Ht_forecast{1}(:,:,h)=diag(H_forecast{1}(h,:))*Rt_forecast{1}(:,:,h)*diag(H_forecast{1}(h,:));
    end
end



