function [parameters, LLF, Ht, Rt, Qt, H, smooth_t, stderrors]=MLE_RS_Multi_DCC_fit(data)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MRS-DCC-Multivariate-GARCH, each single conditional varianec is measured
%% as GARCH(1,1), and the conditional correlation is also measured as
%% DCC(1,1), since under regime shifts framework, higher order has parameter estimation problem.
%% Note that, different from Single-state DCC, in MRS-DCC, one-step
%% estimation prcess is used to control all the variance and correlation shifts probability

% Written by Weijun Xu
% 01/08/2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Let's first set some initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,k]=size(data);

% Individual regime switching Univariate-GARCH parameters set
alpha=0.05;
beta=0.92;
omeg1=(1-(alpha+beta))*1.5*var(data);    %Use variance targeting form to its initial value at regime 1 (with higher variance)
omeg0=(1-(alpha+beta))*0.5*var(data);    %Use variance targeting form to its initial value at regime 0 (with lower variance)
para_start=[];
for i=1:k
    para_start=[para_start; omeg1(:,i);alpha;beta; omeg0(:,i);alpha;beta]; % regime 1 & rgime 0 parameters for each series
end

% DCC regime switching Correlation parameters set
dcc_a1=0.01;  dcc_a0=0.01; % in this case, initial value of the variance-targeting form: (1-dcc_a1-dcc_b1)>(1-dcc_b0-dcc_b0)
dcc_b1=0.8;  dcc_b0=0.97;
tp11=0.9;
tp00=0.9;
dcc_start=[dcc_a1,dcc_b1,dcc_a0,dcc_b0,tp11,tp00]'; % Dcc starting values under two regimes
para_start=[para_start; dcc_start];
num_para=length(para_start);

% Set parameters constraints: alpha(st)+beta(st)<1, dcc_a+dcc_b<1, and each parameter >0, except omeg
Aineq=[];
LB=[(1e-3)*ones(num_para,1)];
UB=(1-1e-3)*ones(num_para,1);
for i=1:k        
    Aineq=[Aineq; zeros(1,6*(i-1)), 0,1,1, zeros(1,num_para-3-6*(i-1));...            % Uni-GARCH regime 1
        zeros(1,6*(i-1)),zeros(1,3),0,1,1,zeros(1,num_para-6-6*(i-1))];               % Uni-GARCH regime 0
    UB([6*(i-1)+1,6*(i-1)+4])=inf;    % all parameters less than 1, while omeg has no uplimit
end
Aineq=[Aineq; zeros(1,6*k), 1,1, zeros(1, num_para-6*k-2); zeros(1,6*k+2),1,1,0,0];   %DCC correlation equation constraints set
Cineq=[(1-1e-3)*ones(2*k+2,1)];

options=optimset('fmincon');
options=optimset(options , 'Display','off','Diagnostics','off','LargeScale','off','MaxIter',80, 'MaxFunEvals', 1e+6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate the Regime-Switching DCC Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[parameters]=fmincon(@(parameters)MLE_RS_Multi_DCC_likelihood(parameters,data),...
    para_start,Aineq,Cineq,[],[],LB,UB,[],options); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate Kim's Smooth probability: Smooth_t(t)=filter_t(t).*{[Smooth_t(t+1)./Conditional_t(t+1)]*P'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[LLF,log_likelihood, Ht, Rt, Qt, H, filter_t, conditional_t]= MLE_RS_Multi_DCC_likelihood(parameters, data);
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 of DCC correlation equation
if nargout>7
    hessian=hessian_2sided('MLE_RS_Multi_DCC_likelihood',parameters, data);
    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_DCC_likelihood(hparameters, data);
        likelihoodsplus(:,i)=indivlike;
    end
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hminus(i);
        [HOLDER, indivlike] = MLE_RS_Multi_DCC_likelihood(hparameters, data);
        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




