function [MLE_coeff, LLF,smooth_t,filter_t]=MLE_RS_uni_fit2(y, x, sigma_type)

% Model: Yt=Xt*B+ St*alpha_st + V(st)^1/2(ut), where V(st)=Sigma_low when
% st==0 and Sigma_high when st==1
[T,k]=size(x);
%% we need to set some smart starting values, without regime we have following paramteres
coeff=-(x\y);                                    % where y=x*coeff+e
sigma0=((y-x*coeff)'*(y-x*coeff)/(T-k)); 
sigma1=1.1*sigma0;
alpha_shift=(x\y);                   % make sure it is always >0
tp11=0.8;                                % make sure it is between 0~1               
tp00=0.8;
if sigma_type==1      % we have regime shifts on variance
    para_start=[alpha_shift; coeff; tp11; tp00; sigma0; sigma1]; 
    LB=[1e-6; -inf*ones(k,1); 1e-6*ones(2,1); 1e-6; 1e-6];
    UB=[inf; inf*ones(k,1);(1-1e-6)*ones(2,1); inf; inf];
    Aineq=[zeros(1, size(para_start,1)-2),1, -1];                   % make sure Sigma1>Sigma0
    Cineq=-1e-6;
else
    para_start=[alpha_shift; coeff; tp11; tp00; sigma0]; 
    LB=[ 1e-6; -inf*ones(k,1);1e-6*ones(2,1); 1e-6];
    UB=[ inf; inf*ones(k,1);(1-1e-6)*ones(2,1); inf];
    Aineq=[];Cineq=[];
end

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

[parameters,LLF,exitflag,output,grad,hessian] = fmincon(@(parameters)MLE_RS_uni_likelihood2(parameters, y, x, sigma_type), ...
   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);

%% Calculate the Filter probability: Pr(st=1|Y_t) & Conditional probability: Pr(St|Y_t-1)
[LLF,log_likelihood, filter_t, conditional_t]=MLE_RS_uni_likelihood2(parameters, y, x,sigma_type);
LLF=-LLF;

%% Re-organize the paramters
reg_coeff=parameters(1:k+1,:);
trans_p=parameters(k+2:k+3,:)
if sigma_type==1
    reg_sigma_low=parameters(end-1,:);
    reg_sigma_high=parameters(end,:);
    reg_sigma=[reg_sigma_low; reg_sigma_high];
else
    reg_sigma_low=parameters(end,:);
    reg_sigma=[reg_sigma_low];
end
MLE_coeff=[reg_coeff; reg_sigma; trans_p];

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sub function 1: {Log-likelihood function}
function   [LLF,log_likelihood, filter_t, conditional_t]= MLE_RS_uni_likelihood2(parameters, y, x, sigma_type)
% 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);
alpha_draw=parameters(1:1,:);        % regime shift intercetpt
beta_draw=parameters(2:k+1,:);
tp11=parameters(k+2:k+2,:);          % transistion probability tp11
tp00=parameters(k+3:k+3,:);          % transistion probability tp00
sigma0_draw=parameters(k+4:k+4,:);   % low volatility variance sigma_draw
if sigma_type==1
    sigma1_draw=parameters(k+5:k+5,:); % high volatility variance sigma1_draw
end

P=[tp11, 1-tp11; 1-tp00, tp00]; % which is the transistion matrix

% calculate the initial steady-state probability
P_transpose=P';
i_m=ones(2,1);
I_m=eye(2);
A=[I_m-P_transpose; i_m'];
temp_mat=((A'*A)^-1) *A';
filter_t(1,:)=temp_mat(:,end)';  % the last column of (A'A)^-1*A' is the initial steady-state probability
ystar=[y(1,:);y];                % extend the y and x 
xstar=[x(1,:);x];
conditional_t=zeros(T+1,2);

for i=2:T+1
    if sigma_type==1   %allow variance shifts     
        likelihood(i,1)=1/sqrt(2*pi*sigma1_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw-alpha_draw*1).^2/(sigma1_draw));       % when st=1
        likelihood(i,2)=1/sqrt(2*pi*sigma0_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw).^2/(sigma0_draw));                    % when st=0
    else               % No varinace shifts
        likelihood(i,1)=1/sqrt(2*pi*sigma0_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw-alpha_draw*1).^2/(sigma0_draw));       % when st=1
        likelihood(i,2)=1/sqrt(2*pi*sigma0_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw).^2/(sigma0_draw));                    % when st=0
    end    
    if sum(likelihood(i,:)==0)>1        % i.e. two likelihood are both zero
        likelihood(i,:)=1e-30.*ones(2,1);
    end
    conditional_t(i,:)=filter_t(i-1,:)*P;    
    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);
filter_t=filter_t(2:end,:);
conditional_t=conditional_t(2:end,:);




        

