function [parameters, LLF, smooth_t,filter_t, cov_matrix, stderrors]=MLE_RS_Multi_intercept_SUR_fit(y, x1, x2)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model: Y1t=X1t*B + St*Alpha + V(st)^1/2(ut), where V(st)=Sigma_low when st=0 and Sigma_high when st=1
%%        Y2t=X2t*B + St*Alpha + V(st)^1/2(ut), 
%% MLE approach to test MRS-SUR model, allow regime shifts in intercepts and cov-var matrix
%% Constraints: Alpha < 0 and varinace in regime 1 > regime 0
% Written By Weijun Xu, School of Banking & Finance, UNSW
% First edition, 01/08/08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T,k]=size(x1);
num_eq=size(y,2);
num_cov=num_eq*(num_eq+1)/2;   %num of covariance 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)=my_kron2(x1(i,:), x2(i,:)); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% simple way to get seemingly unrelated regression coefficients and cov
% matrix, though not precise, just set individual equation coefficients as
% initial values, since later on we still have correct MLE way to estimate
OLS_coeff1=x1\y(:,1);  % equation 1
OLS_coeff2=x2\y(:,2);  % equation 2
resid1=y(:,1)-x1*OLS_coeff1;
resid2=y(:,2)-x2*OLS_coeff2;
rho12=corr(resid1,resid2);
OLS_sigma1=resid1'*resid1/(T-k);
OLS_sigma2=resid2'*resid2/(T-k);
OLS_sigma=[OLS_sigma1, rho12*sqrt(OLS_sigma1*OLS_sigma2); rho12*sqrt(OLS_sigma1*OLS_sigma2), OLS_sigma2];

alpha_draw=[-1e-3*ones(num_eq,1)];             % regime shift intercetpts, High volatility has low intercept
beta_draw=[OLS_coeff1(:); OLS_coeff2(:);];

sigma0_chol=chol(OLS_sigma)';
sigma0_chol=0.8*nonzeros(sigma0_chol);         % stack the low-trangular matrix together in a column vector
sigma1_chol=1.2*chol(OLS_sigma)';
sigma1_chol=nonzeros(sigma1_chol);             % Regime 1 is high volatility regime
tp11=0.9;
tp00=0.9;
para_start=[alpha_draw; beta_draw; sigma0_chol; sigma1_chol; tp11; tp00];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the constraints:  (1). alpha<0 &  (2). tp11, tp00 ~(0,1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LB=[-inf*ones(num_eq,1);-inf*ones(k*num_eq,1);-inf*ones(2*num_cov,1);(1e-6)*ones(2,1)];
UB=[-1e-6*ones(num_eq,1);+inf*ones(k*num_eq,1); inf*ones(2*num_cov,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_intercept_likelihood(parameters, yt, xt, T, num_eq, num_cov,k),...
   para_start,[],[],[],[],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, filter_t, conditional_t]= MLE_RS_Multi_intercept_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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma0_chol=parameters(num_eq*(1+k)+1:num_eq*(1+k)+num_cov,:);                % low tranigular matrix of low volatility cov-var matrix
sigma1_chol=parameters(num_eq*(1+k)+num_cov+1:num_eq*(1+k)+2*num_cov,:);      % low tranigular matrix of high volatility cov-var matrix
% construct the cov-var matrix
sigma0_chol=ivech(sigma0_chol);
sigma1_chol=ivech(sigma1_chol);               % transfer the column vector to low-trangular matrix
sigma0_draw=sigma0_chol*sigma0_chol'; 
sigma1_draw=sigma1_chol*sigma1_chol';         % transfer low-trangular matrix to cov-variance matrix
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_intercept_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_intercept_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_intercept_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}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [C,Ceq]=my_con(parameters,num_eq, num_cov,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=[];

