function [parameters, LLF,smooth_t,filter_t, tp11, tp00]=MLE_RS_uni_TVTP_fit_all(y, x, z)

% This function allows each parameters to have regime shifts in uni-OLS
% model, with time varying transistion probability
% Model: Yt=Xt*B_st+ St*alpha_st + V(st)^1/2(ut), where V(st)=Sigma_low when
% st==0 and Sigma_high when 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)))
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Set initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,k]=size(x);
reg_coeff=(x\y);                                    % where y=x*coeff+e
coeff0=[reg_coeff(1,:); reg_coeff(2:end,:)];              % regime 0, low volatility, high intercept
coeff1=[-reg_coeff(1,:); reg_coeff(2:end,:)];             % regime 1, high volatility, low intercept
OLS_sigma=((y-x*reg_coeff)'*(y-x*reg_coeff)/(T-k)); 
sigma0=OLS_sigma;
sigma1=1.2*OLS_sigma;
tp_coeff=[4; 1; 1]; % Try to make Pr(St=1|St-1=1)~0.8 % [1.5 1.5 1] is OK for monthly      
para_start=[coeff0; coeff1; sigma0; sigma1; tp_coeff];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Set Constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num_coeff=size(para_start,1);
LB=[-inf*ones(2*k,1); 1e-6*ones(2,1); -inf*ones(3,1)];
UB=[inf*ones(2*k,1);inf*ones(2,1);inf*ones(3,1)];                                     
Aineq=[-0.8, 0, 1, 0, zeros(1,num_coeff-4);...                     % intercept1<intercept0
    zeros(1,4), 1, -1, zeros(1,num_coeff-6)];                    % Sigma1>Sigma0,
Cineq=[-1e-6; -1e-6];

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

[parameters] = fmincon(@(parameters)MLE_RS_uni_likelihood_all(parameters, y, x,z), ...
   para_start,Aineq,Cineq,[],[],LB,UB,[],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 inference probability for Pr(st=1|Yt)
[LLF,log_likelihood, filter_t, conditional_t, tp11, tp00, P]= MLE_RS_uni_likelihood_all(parameters, y, x, z);
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_uni_likelihood_all(parameters, y, x, z)
% Using MLE to estimate the univariate regime switching models
% sigma_type allows to choose variance type, '1' means varaince shifts, '0' means no variance shifts.
% conditional_t=Pr(St|Y_t-1)
% filter_t=Pr(St|Y_t)

[T,k]=size(x);
beta0_draw=parameters(1:k,:);        % coeff at regime 0
beta1_draw=parameters(k+1:2*k,:);    % coeff at regime 1
sigma0_draw=parameters(2*k+1:2*k+1,:);     % low volatility variance sigma0_draw
sigma1_draw=parameters(2*k+2:2*k+2,:);     % high volatility variance sigma0_draw
tp_coeff=parameters(2*k+3:end,:);     
tp11=zeros(T,1);
tp00=zeros(T,1);
% calculate 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)=1/sqrt(2*pi*sigma1_draw)*...
        exp(-0.5*(y(i,1)-x(i,:)*beta1_draw).^2/(sigma1_draw));       % when st=1
    likelihood(i,2)=1/sqrt(2*pi*sigma0_draw)*...
        exp(-0.5*(y(i,1)-x(i,:)*beta0_draw).^2/(sigma0_draw));       % when st=0
    if sum(likelihood(i,:)==0)>1        % i.e. two likelihood are both zero
        likelihood(i,:)=1e-30.*ones(2,1);
    end
   
    tp00(i)=exp(tp_coeff(1)+tp_coeff(3)*z(i-1))/(1+exp(tp_coeff(1)+tp_coeff(3)*z(i-1)));   
    tp11(i)=exp(tp_coeff(2)+tp_coeff(3)*z(i-1))/(1+exp(tp_coeff(2)+tp_coeff(3)*z(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);

        

