function   [smpl_coeff, smpl_tp, smpl_st,smpl_p_st, smpl_ip_t, smpl_marginal_t ]= ...
    Bay_RS_intercept_uni_fit(y, x, nrep, warmup, sigma_type)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model: Yt=Xt*B(St) +  V(St)^1/2(ut), where V(st)=Sigma_low when st=0 and Sigma_high when st=1
%% Bayesian Gibbs-Sampling for MRS-OLS model, allow regime shifts on: intercept and varaiance
%% Have constraints: (1). intercept in regime 1 < intercept in regime 2
%%                   (2). Variance in regime 1 > Variance in regime 2
% estimate general regression for univaraite analysis under two regimes cases
% All priors are given below, one can change some initial values, while
% since I use Gibbs-sampling, the results are less sensitive to the intial values
% Note: for the scale y and x, when they represents the return (which is usually the case), it is bettet to use 
% return itself not the return percentage//

% Written by Weijun Xu, 18/07/2008
% School of Banking and Finance, UNSW
% First edition

[T,k]=size(x);
%% Setting the priors: Assume there is only two latent states
%%%%% Transition probability prior %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a11=5; b10=1;                    % tp11~Beta(a11, b10)                   
a00=5; b01=1;                    % tp00~Beta(a00, b01), tp10->St=0|St-1=1
% the tranition matrix is like:          St=1            St=0
%                               St-1=1   tp11            (1-tp11)
%       `                           St-1=0   (1-tp00)         tp00
%%%%%%% Coefficients beta %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ub=0*ones(k,1); vb=5*eye(k);  % Beta1~N(u1b, v1b)
%%%%%%% Coefficients regime-shift alpha %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ualpha=0.5; valpha=5;    
%%%%%%% Varianec low volatile regime        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=2; B=1;              % sigma~IG(A, B)
%%%%%%% Varianec high volatile regime, =sigma(1+omeg*st)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=1; b=3;              % 1+omeg~IG(a, b), 1+omeg>1

%% Setting initial values
beta_draw=multi_normrnd(ub, vb,1)';
alpha_draw=multi_normrnd(ualpha, valpha,1)';
%sigma_draw=gamrnd(A,1;    omeg_draw=gamrnd(a,1); 
%sigma_draw=B/sigma_draw;  omeg_draw=b/omeg_draw;                           
sigma_uncond=cov(y);
sigma_draw=sigma_uncond;
omeg_draw=3;
tp11_draw=0.9; tp00_draw=0.9;

%% Gibbs-Sampling
for iter=1:nrep+warmup   
    % step 1: draw latent state variables: st  
    if sigma_type==1
        sigma=[sigma_draw; omeg_draw];
    else 
        sigma=sigma_draw;
    end    
    [st_draw, p_st, ip_t, marginal_t]=state_draw_uni(y,x,beta_draw,alpha_draw, sigma, tp11_draw, tp00_draw, sigma_type);
    % step 2: draw the posterior coefficients: beta
    % we need do some data transformation to make the beta is linear
    for i=1:T
        if sigma_type==1
            ystar(i,1)=(1+omeg_draw*st_draw(i))^-0.5 * y(i,1);
            xstar(i,:)=(1+omeg_draw*st_draw(i))^-0.5 .* x(i,:);
            st_star(i,1)=(1+omeg_draw*st_draw(i))^-0.5 * st_draw(i,1);
        else
            ystar(i,1)= y(i,1);
            xstar(i,:)=x(i,:);
            st_star(i,1)=st_draw(i,1);
        end

    end
    y_new=ystar-alpha_draw.*st_star;   % since we are sampling beta, we treat (ystar-alpha * st_draw) as new 'y'
    x_new=xstar;
    Db=(x_new'*x_new/sigma_draw+inv(vb))^-1;
    db=x_new'*y_new/sigma_draw+inv(vb)*ub;
    beta_draw=multi_normrnd(Db*db, Db,1)';
    clear y_new x_new

    % step 3: draw the posterior regime shift intercept, alpha
    % Notice that, we assume regime 1 is high level, high volatile, thus
    % alpha should be greater than 0
    
    y_new=ystar-xstar*beta_draw;   % since we are sampling RS-intercept now, we treat (ystar-xstar*beta) as new 'y'
    x_new=st_star;                 % and we treat st_star as new 'x' now
    Dalpha=(x_new'*x_new/sigma_draw+inv(valpha))^-1;
    dalpha=x_new'*y_new/sigma_draw+inv(valpha)*ualpha;
    alpha_draw=multi_normrnd(Dalpha*dalpha, Dalpha,1)';
    % since we need alpha>0, we need to use truncated sample
    % alpha_draw=mu+sigma*(IV~CDF(U*(P2-P1)+P1)
    % where P1=CDF[(ll-mu)/sigma]; P2=CDF[(ul-mu)/sigma];
    ll=1e-6; ul=1e+6; % for alpha>0
    u=unifrnd(0,1);
    p1=normcdf((ll-Dalpha*dalpha)/sqrt(Dalpha));
    p2=normcdf((ul-Dalpha*dalpha)/sqrt(Dalpha));
    alpha_draw=Dalpha*dalpha+sqrt(Dalpha)* norminv(u*(p2-p1)+p1,0,1);
    
    clear y_new x_new
   
    % step 4: draw the posterior low volatility variance: sigma      
    shape=A+T/2;
    scale=B+ sum((ystar-xstar*beta_draw-alpha_draw.*st_star).^2)/2;
    sigma_draw=gamrnd(shape,1);
    sigma_draw=scale/sigma_draw;
    
    % step 5: draw the posterior high volatilty: omeg, where
    % sigma_high=sigma(1+omeg*1)
    % agian, we need to do some data transformation
    % ybar=xbar*beta+alpha*st_bar +(1+omeg)^.5 *ut, we treat
    % omegbar~IG(shape, scale), and then calculate omeg_draw
    if sigma_type==1
        obs_1=sum(st_draw,1);
        ybar=y((st_draw==1),:)./sqrt(sigma_draw);
        xbar=x((st_draw==1),:)./sqrt(sigma_draw);
        st_bar=st_draw((st_draw==1),:)./sqrt(sigma_draw);
        shape2=a+obs_1/2;
        scale2=b+ sum((ybar-xbar*beta_draw- alpha_draw.*st_bar).^2)/2;
        omeg2=gamrnd(shape2,1);
        omeg2=scale2/omeg2;
        % make sure omeg2>1, since we assume state 1 has high volatility
        while omeg2<1
            omeg2=gamrnd(shape2,1);
            omeg2=scale2/omeg2;
        end
        omeg_draw=omeg2-1;    % since we need omeg_draw not omeg2
    end

    % step 6: draw the transistion probability: tp11,and tp00
    [n11, n10, n01, n00]=sigma_draw_number_uni(st_draw);  % n10 is number of states that St=0|St-1=1
    tp11_draw=betarnd(a11+n11, b10+n10);   % tp11~Beta(a11+n11, b10+n10), where n10=number of times that St-1=1|St=0
    tp00_draw=betarnd(a00+n00, b01+n01);   % tp00~Beta(a00+n00, b01+n01)
        
    if iter>warmup
        beta_keep(iter-warmup,:)=beta_draw';
        alpha_keep(iter-warmup,:)=alpha_draw';
        sigma_keep(iter-warmup,:)=sigma_draw;
        if sigma_type==1
            sigma_high_keep(iter-warmup,:)=sigma_draw*(1+omeg_draw*1); % high volatility variance
        end
        tp11_keep(iter-warmup,:)=tp11_draw;
        tp00_keep(iter-warmup,:)=tp00_draw;
        st_keep(iter-warmup,:)=st_draw';
        smpl_p_st(iter-warmup,:)=p_st';
        smpl_ip_t(iter-warmup,:)=ip_t';
        smpl_marginal_t(iter-warmup,:)=marginal_t';
        
    end
    fprintf(' Iteration to %3g \n', iter)
end
if sigma_type==1
    smpl_coeff=[alpha_keep, beta_keep, sigma_keep, sigma_high_keep];
else
    smpl_coeff=[alpha_keep, beta_keep, sigma_keep];
end
smpl_tp=[tp11_keep, tp00_keep];
smpl_st=st_keep;
end



%% Subfunction 1: Generate latent state variable, the approach used here
%% follows Chib (1996), which is quite like Original Markov-Switching
%% Hamilton (2005)'s inference probability calculation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [st_draw, p_st, ip_t, marginal_t]=state_draw_uni(y,x,beta_draw,alpha_draw, sigma, tp11_draw, tp00_draw, sigma_type)

% Usage: to generate the latent state variable for uni-variate regression/
% this function is only for two state case
if sigma_type==1
    sigma_draw=sigma(1,:);
    omeg_draw=sigma(2,:);
else
    sigma_draw=sigma(1,:);
end

T=size(y,1);
st_draw=zeros(T,1);
U=unifrnd(0,1,T,1);

% first let's define some initial settings
P=[tp11_draw, 1-tp11_draw; 1-tp00_draw, tp00_draw]; % which is the transistion matrix
ip_t=zeros(T,2);       % which is the inference proability: ip_t(t,1)=Pr(St=1|Yt)
marginal_t=zeros(T,2); % marginal_t=Pr(St=k|Yt-1)    
p_st=zeros(T,2);
ip_t(1,:)=[0.5, 0.5];  % assign the initial ip_t value
marginal_t(1,:)=[0.5, 0.5];

for i=2:T
    if sigma_type==1        
        likelihood(i,1)=exp(-0.5*(y(i,1)-x(i,:)*beta_draw-alpha_draw*1).^2/(sigma_draw*(1+omeg_draw *1)));  % when st=1
        likelihood(i,2)=exp(-0.5*(y(i,1)-x(i,:)*beta_draw-alpha_draw*0).^2/(sigma_draw*(1+omeg_draw *0)));  % when st=0
    else
        likelihood(i,1)=exp(-0.5*(y(i,1)-x(i,:)*beta_draw-alpha_draw*1).^2/(sigma_draw));  % when st=1
        likelihood(i,2)=exp(-0.5*(y(i,1)-x(i,:)*beta_draw-alpha_draw*0).^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);
        %likelihood(i,:)=likelihood(i-1,:);
    end
    marginal_t(i,:)=ip_t(i-1,:)*P;
    ip_t(i,:)=marginal_t(i,:).*likelihood(i,:);
    ip_t(i,:)=ip_t(i,:)./sum(ip_t(i,:),2);  % Normalized ip_t
end
st_draw(end)=(U(end)<=ip_t(end,1));  % if random uniform number is less than ip_t, we set it as 1, otherwise set it as 0
% caculate the posterior conditional probability: P(st|Y_T)
for i=T-1:-1:1
    if st_draw(i+1)==1 
        p_st(i,:)=ip_t(i,:).*P(:,1)';   % -->[ip_1*p11, ip_0* p01]
        p_st(i,:)=p_st(i,:)./sum(p_st(i,:),2);   % Normalize p_st1
    elseif st_draw(i+1)==0
        p_st(i,:)=ip_t(i,:).*P(:,2)';   % -->[ip_1*p10, ip_0* p00]
        p_st(i,:)=p_st(i,:)./sum(p_st(i,:),2);
    end            
    st_draw(i)=(U(i)<=p_st(i,1));
end

marginal_t=marginal_t(:,1);       % marginal probabiilty:  P(St=1|Y_t-1
ip_t=ip_t(:,1);                   % inference probability: P(St=1|Y_t) 
p_st(end,:)=ip_t(end,1);
p_st=p_st(:,1);                   % posterior density:     P(St=1|Y_T, S_T+1)
end
    

%% Subfunction 2: calculate the some count number for inverse garma to
%% generate sigmas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [n11, n10, n01, n00]=sigma_draw_number_uni(st_draw)
% to calculate the count number n11, n10, n01, n00;
% generally, one count will be missing in the total sum of counts
T=size(st_draw,1);
n11=0; n10=0; n01=0; n00=0;
for i=2:T
    if  st_draw(i)==1
        if st_draw(i-1)==1
            n11=n11+1;     % n10 is number of times that St=0|St-1=1
        else
            n01=n01+1;
        end
    elseif  st_draw(i)==0
        if st_draw(i-1)==0
            n00=n00+1;
        else
            n10=n10+1;
        end
    end
end
end

   





















