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

% Model: Yt=Xt*B+ St*alpha_st + V(st)^1/2(ut), where V(st)=Sigma(1+omeg)
[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
sigma=(y-x*coeff)'*(y-x*coeff)/(T-k); 
alpha_shift=0.8*(x\y);                   % make sure it is always >0
tp11=0.8;                            % we use tp11=(exp(tp11)/(1+exp(tp11)) to make sure it is between 0~1               
tp00=0.8;
if sigma_type==1      % we have regime shifts on variance
    omeg=0.5;           % we are going to use omeg=exp(log_omeg) to make sure always omeg>0 
    para_start=[alpha_shift; coeff; tp11; tp00; sigma; omeg]; 
    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];    
else
    para_start=[alpha_shift; coeff; tp11; tp00; sigma]; 
    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];
end

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

[parameters,LLF,exitflag,output,grad,hessian] = fmincon(@(parameters)MLE_RS_intercept_uni_likelihood(y, x, parameters,sigma_type), ...
   para_start,[],[],[],[],LB,UB,[],options);

%% calculate the standard errors

%% get the inference probability for Pr(st=1|Yt)
[LLF,log_likelihood, filter_t]=MLE_RS_intercept_uni_likelihood(y, x, parameters,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=reg_sigma_low*(1+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];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Log-likelihood function 
function   [LLF,log_likelihood, filter_t]= MLE_RS_intercept_uni_likelihood(y, x, parameters,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.

[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
sigma_draw=parameters(k+4:k+4,:);   % low volatility variance sigma_draw

if sigma_type==1
    omeg_draw=parameters(k+5:k+5,:);   % high volatility variance proportion, where simga_high=sigma(1+omeg*1)
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];

test=0;
for i=2:T+1
    if sigma_type==1        
        likelihood(i,1)=1/sqrt(2*pi*sigma_draw*(1+omeg_draw *1))*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw-alpha_draw*1).^2/(sigma_draw*(1+omeg_draw *1)));  % when st=1
        likelihood(i,2)=1/sqrt(2*pi*sigma_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw).^2/(sigma_draw));                                % when st=0
    else
        likelihood(i,1)=1/sqrt(2*pi*sigma_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw-alpha_draw*1).^2/(sigma_draw));                   % when st=1
        likelihood(i,2)=1/sqrt(2*pi*sigma_draw)*...
            exp(-0.5*(ystar(i,1)-xstar(i,:)*beta_draw).^2/(sigma_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
    
    joint_likelihood(i,1)=sum((filter_t(i-1,:)*P).*likelihood(i,:));
    filter_t(i,:)=((filter_t(i-1,:)*P).*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,:);

