function [parameters,LLF,smooth_t,filter_t]=MLE_RS_Multi_TVTP_fit(y, x, z)

% This function can be tested in Zhou(2008), Ang(2002)'s papers
% MLE Multivariate Regime-Switching Estimation with time varying transition probability
% This model allow to calculate the conditional mean, co-variance, co-skewness...
% We use short-term interest rate as prediction rate
% Model: y_1t=u_1(st)+b_1(st)*Tbill_t-1+e_1t(st)
%        y_2t=u_2(st)+b_2(st)*Tbill_t-1+e_2t(st) 
%       where H_it=Dt(st)*Rt(st)*Dt(st), with
%       Dt(st)=diag{sqrt[h1t(st),h2t(st)]}, Rt(st)=[1 rho(st); rho(st) 1]
% Pr(St=1|St-1=1,y_t-1)=exp(a1+b1*Tbill(t-1))/(1+exp(a1+b1*Tbill(t-1)))
% Pr(St=2|St-1=2,y_t-1)=exp(a2+b2*Tbill(t-1))/(1+exp(a2+b2*Tbill(t-1))) 
% st=1 or 0, where 1 represents high volatility regime with low intercept,
% and 0 low volatility regime with high intercept


[T,k]=size(x);
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 
zt=z;
for i=1:T                      % Xt: num_eq x k x T, three dimension matrix
    xt(:,:,i)=kron(eye(num_eq), x(i,:)); 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OLS_coeff=x\y;
OLS_H=(y-x*OLS_coeff)'*(y-x*OLS_coeff)/(T-k*num_eq); 
OLS_D=diag(diag(OLS_H)); 
OLS_R=inv(sqrt(OLS_D))*OLS_H*inv(sqrt(OLS_D));

% Two regime shifts coefficients
beta1_draw=[-OLS_coeff(1,:); OLS_coeff(2:end,:)]; % since regime 1 is high volatility with low intercept
beta1_draw=beta1_draw(:);
beta0_draw=OLS_coeff(:);
% two diagonal varaince matrix and correlation matrix
sigma1_D=1.1*OLS_D;  sigma1_D=nonzeros(sigma1_D);
rho1=OLS_R;          rho1=rho1(1,2);
sigma0_D=OLS_D;  sigma0_D=nonzeros(sigma0_D);
rho0=OLS_R;          rho0=0.9*rho0(1,2);
% transition probability parameters
tp_coeff=[2;1;-1];       
para_start=[beta1_draw; beta0_draw; sigma1_D; sigma0_D; rho1; rho0; tp_coeff];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Set the constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters=[b10;b11;b20;b21;sigma1_1;sigma1_2;sigma2_1;sigma2_2;rho1;rho2;a1;b1;a2;b2]
num_coeff=size(para_start,1);
LB=[-inf*ones(2*k*num_eq,1); 1e-6*(ones(num_cov*2,1)); -inf*ones(3,1)];
UB=[inf*ones(2*k*num_eq,1); inf*ones(2*num_eq,1); (1-1e-6)*ones(num_eq,1); inf*ones(3,1)]; % sigma>0, and 0<rho<1
Aineq=[1,zeros(1,3),-1,zeros(1,num_coeff-5);...            % regime 1's intercept lower than regime 2's intercept
    0,1,zeros(1,3),-1,zeros(1,num_coeff-6);...
    zeros(1,8), -1,0,1,zeros(1,num_coeff-11);...           % regime 1's sigma1_1 higher than regime 2's sigma2_1
    zeros(1,8),0,-1,0,1,zeros(1,num_coeff-12);...          % regime 1's sigma1_2 higher than regime 2's sigma2_2
    zeros(1,8),zeros(1,4),-1,1,zeros(1,num_coeff-14)];     % regime 1's rho1 higher than regime 2's rho2
Cineq=[-1e-2*ones(2,1);-1e-6*ones(3,1)];

options=optimset('fmincon');
options=optimset(options,'Display','iter','Diagnostics','off','LargeScale','off');

[parameters] = fmincon(@(parameters)MLE_RS_Multi_TVTP_likelihood(parameters, yt, xt,zt, T, num_eq, num_cov,k),...
   para_start,Aineq,Cineq,[],[],LB,UB,[],options);
%@(parameters)my_con(parameters,num_eq, num_cov,k),options);

%% calculate the standard errors
%[stder_general, stder_BHHH, stder_White, scores, hessian2]= std_error('MLE_RS_uni_likelihood2',parameters, y, x,sigma_type);

%% get the Filter probability for Pr(st=1|Yt) & conditioal probability for Pr(st=1|Yt-1)
[LLF,log_likelihood, filter_t, conditional_t, tp11, tp00, P]= MLE_RS_Multi_TVTP_likelihood(parameters, yt, xt,zt, T, num_eq, num_cov,k);
LLF=-LLF;

%% Calculate the Kim's smoothed probability
% Smooth_t(t)=filter_t(t).*{[Smooth_t(t+1)./Conditional_t(t+1)]*P'};
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(:,:,i)');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sub function 1: {Log-likelihood function}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function   [LLF,log_likelihood, filter_t, conditional_t, tp11, tp00, P]= MLE_RS_Multi_TVTP_likelihood(parameters, yt, xt,zt, T, num_eq, num_cov,k)
% Using MLE to estimate the Multivariate regime switching models
% sigma_type allows to choose variance type, '1' means varaince shifts, '0' means no variance shifts.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% % assign the parameter coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters=[b10;b11;b20;b21;sigma1_1;sigma1_2;sigma2_1;sigma2_2;rho1;rho2;a1;a2;b]

beta1_draw=parameters(1:k*num_eq,:);
beta0_draw=parameters(k*num_eq+1:2*k*num_eq,:);
sigma1_D=parameters(2*k*num_eq+1:(2*k+1)*num_eq,:); 
sigma0_D=parameters((2*k+1)*num_eq+1:(2*k+2)*num_eq,:);
rho1=parameters((2*k+2)*num_eq+1:(2*k+2)*num_eq+1,:);
rho0=parameters((2*k+2)*num_eq+2:(2*k+2)*num_eq+2,:);
tp_coeff=parameters((2*k+2)*num_eq+2+1:(2*k+2)*num_eq+2+3,:);


sigma1_draw=diag(sqrt(sigma1_D))*[1 rho1;rho1 1]*diag(sqrt(sigma1_D)); % high volatility regime cov-var matrix
sigma0_draw=diag(sqrt(sigma0_D))*[1 rho0;rho0 1]*diag(sqrt(sigma0_D)); % low volatility regime cov-var matrix

% Set the initial steady-state probability
filter_t(1,:)=[0.5, 0.5];  
conditional_t(1,:)=[0.5, 0.5]; 
tp11(1)=0.8; tp00(1)=0.8;
P(:,:,1)=[tp11(1), 1-tp11(1); 1-tp00(1), tp00(1)];

for i=2:T
    likelihood(i,1)=(2*pi)^(-1/2*num_eq)*(det(sigma1_draw)).^(-1/2)*...
        exp(-1/2*(yt(:,i)-xt(:,:,i)*beta1_draw)'*inv(sigma1_draw)*...
        (yt(:,i)-xt(:,:,i)*beta1_draw));     % when st=1, we use sigma1_draw and beta1_draw
    likelihood(i,2)=(2*pi)^(-1/2*num_eq)*(det(sigma0_draw)).^(-1/2)*...
        exp(-1/2*(yt(:,i)-xt(:,:,i)*beta0_draw)'*inv(sigma0_draw)*...
        (yt(:,i)-xt(:,:,i)*beta0_draw));     % when st=0, we use sigma0_draw, and beta0_draw
    
    if sum(likelihood(i,:)==0)>1        % i.e. two likelihood are both zero, avoid 0/0 problem
        likelihood(i,:)=1e-30.*ones(2,1);
    end
    % time varying transtion probability
    
    tp00(i)=exp(tp_coeff(1)+tp_coeff(3)*zt(i-1))/(1+exp(tp_coeff(1)+tp_coeff(3)*zt(i-1)));
    tp11(i)=exp(tp_coeff(2)+tp_coeff(3)*zt(i-1))/(1+exp(tp_coeff(2)+tp_coeff(3)*zt(i-1)));
    P(:,:,i)=[tp11(i), 1-tp11(i); 1-tp00(i), tp00(i)];               % This is the transistion matrix

    conditional_t(i,:)=filter_t(i-1,:)*P(:,:,i);
    joint_likelihood(i,1)=sum(conditional_t(i,:).*likelihood(i,:));
    filter_t(i,:)=(conditional_t(i,:).*likelihood(i,:))./joint_likelihood(i,1);   % Normalized filter_t
end
log_likelihood=log(joint_likelihood(2:end,:));
LLF=-sum(log_likelihood);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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=[];




