function [parameters, LLF, smooth_t,filter_t, cov_matrix, stderrors]=MLE_RS_Multi_all_fit(y, x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model: Yt=Xt*B(st)+ V(st)^1/2(ut), where V(st)=Sigma_low when st=0 and
%% Sigma_high when st=1, we assume constant correlation matrix!!

% Used to test the Regime-Switching Mulitrivate case, through MLE approach
% Can be used to test Zhou(2008)'s paper
% Written By Weijun Xu, School of Banking & Finance, UNSW
% First edition, 01/08/08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T,k]=size(x);
num_eq=size(y,2);
num_cov=num_eq*(num_eq+1)/2;   %num of covariance matrix parameters need to be estimated

yt=y';                         % Yt: num_eq x T           
for i=1:T                      % Xt: num_eq x k x T, three-D matrix
    xt(:,:,i)=kron(eye(num_eq), x(i,:)); 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OLS_coeff=x\y;
OLS_sigma=(y-x*OLS_coeff)'*(y-x*OLS_coeff)/(T-k);

beta1_draw=[-OLS_coeff(1,:);OLS_coeff(2:end,:)]; % regime shift intercetpts, High volatility has low intercept
beta1_draw=beta1_draw(:);
beta0_draw=OLS_coeff(:);

sigma1_chol=1.2*chol(OLS_sigma)';
sigma1_chol=nonzeros(sigma1_chol);               % Regime 1 is high volatility regime
sigma0_chol=chol(OLS_sigma)';
sigma0_chol=0.9*nonzeros(sigma0_chol);           % stack the low-trangular matrix together in a column vector

tp11=0.9;
tp00=0.9;
para_start=[beta1_draw; beta0_draw; sigma1_chol; sigma0_chol; tp11; tp00];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the constraints:  intercept in regime 1 < intercept in regime 2  & tp11, tp00 ~(0,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num_para=length(para_start);
Aineq=zeros(num_eq,num_para);
for i=1:num_eq     % intercept in regime 1 < intercept in regime 2
    Aineq(i,:)=[zeros(1, (i-1)*k), 1, zeros(1, num_eq*k-1), -1, zeros(1, num_para-(i-1)*k-num_eq*k-1)];
end
Cineq=-1e-6*ones(num_eq,1);

LB=[-inf*ones(num_para-2,1);(1e-6)*ones(2,1)];   % Constraint tp11, tp00 ~(0,1)
UB=[inf*ones(num_para-2,1); (1-1e-6)*ones(2,1)]; 

options=optimset('fmincon');
options=optimset(options,'Display','iter','Diagnostics','off','LargeScale','off','DiffMaxChange', 1e-2);

[parameters] = fmincon(@(parameters)MLE_RS_Multi_all_likelihood(parameters, yt, xt, T, num_eq, num_cov,k),...
   para_start,[],[],[],[],LB,UB,@(parameters)my_con(parameters, num_eq, num_cov,k),options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate Kim's Smooth probability: Smooth_t(t)=filter_t(t).*{[Smooth_t(t+1)./Conditional_t(t+1)]*P'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[LLF,log_likelihood, filter_t, conditional_t]= MLE_RS_Multi_all_likelihood(parameters, yt, xt, T, num_eq, num_cov, 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Transfer the low-tranigular paramters to covariance-matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma1_chol=parameters(2*num_eq*k+1:2*num_eq*k+num_cov,:);                % low tranigular matrix of low volatility cov-var matrix
sigma0_chol=parameters(2*num_eq*k+num_cov+1:2*num_eq*k+2*num_cov,:);      % low tranigular matrix of high volatility cov-var matrix

% construct the cov-var matrix
sigma1_chol=ivech(sigma1_chol);               % transfer the column vector to low-trangular matrix
sigma0_chol=ivech(sigma0_chol);
sigma1_draw=sigma1_chol*sigma1_chol';         % transfer low-trangular matrix to cov-variance matrix
sigma0_draw=sigma0_chol*sigma0_chol'; 

cov_matrix=zeros(num_eq,num_eq,3);
cov_matrix(:,:,1)=sigma0_draw;                % low-volatility-cov-matrix
cov_matrix(:,:,2)=sigma1_draw;                % high-volatility-cov-matrix
cov_matrix(:,:,3)=OLS_sigma;                  % single regime cov-matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the standard errors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Hessian
if nargout>5
    fprintf('\n Calculating the standard errors.... \n')
    hessian=hessian_2sided('MLE_RS_Multi_all_likelihood',parameters, yt, xt, T, num_eq, num_cov,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_all_likelihood(hparameters, yt, xt, T, num_eq, num_cov,k);
        likelihoodsplus(:,i)=indivlike;
    end
    for i=1:length(parameters)
        hparameters=parameters;
        hparameters(i)=hminus(i);
        [HOLDER, indivlike] =MLE_RS_Multi_all_likelihood(hparameters, yt, xt, T, num_eq, num_cov,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} Varince in regime 1 > variance in regime 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [C,Ceq]=my_con(parameters, num_eq, num_cov,k)

% the constraint is that the diagonal of covariance matrix of regime 1 is greater
% than diagonal of covariance matrix regime 2
% i.e. volatility in regime 1 is larger than regime 2

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

% construct the cov-var matrix
sigma1_chol=ivech(sigma1_chol);               % transfer the column vector to low-trangular matrix
sigma0_chol=ivech(sigma0_chol);
sigma1_draw=sigma1_chol*sigma1_chol';         % transfer low-trangular matrix to cov-variance matrix
sigma0_draw=sigma0_chol*sigma0_chol'; 

sigma1=diag(sigma1_draw);  
sigma0=diag(sigma0_draw); 
C=sigma0-sigma1;                              % each element of varaince of reigme 2 is less than regime 1
Ceq=[];


