function   [smpl_coeff, smpl_st, smooth_t_keep, filter_t_keep, in_sample_results,hedge_insample]= ...
    Bay_RS_uni_all_OLS_fit_insample(y, x, nrep, warmup)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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 shift everything: intercept, coefficients, and variance
%% Have constraints: (1). Variance in regime 1 > Variance in regime 2
% Written By Weijun Xu, School of Banking & Finance, UNSW
% First edition, 01/08/08
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setting Priors & initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,k]=size(x);
% Transition probability prior
a11=10; b10=1;                    % tp11~Beta(a11, b10)
a00=10; 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
OLS_coeff=(x\y);
OLS_sigma=(y-x*OLS_coeff)'*(y-x*OLS_coeff)/(T-k);
OLS_coeff_stdev=OLS_sigma*inv(x'*x);
ub0=OLS_coeff; vb0=10*eye(k);           % Beta0~N(u1b, v1b)
ub1=OLS_coeff; vb1=10*eye(k);           % Beta1~N(u1b, v1b)

% Varianec low volatile regime
A0=10; B0=0.8*A0*OLS_sigma;        % Sigma0~IG(A0, B0), the higher the A0, the larger the volatile sample
A1=10; B1=1.2*A1*OLS_sigma;        % Sigma1~IG(A1, B1), Make sure Sigma1>Sigma0, since we assume regime 1 is high volatility regime

% Setting initial values
beta1_draw=multi_normrnd(ub1, vb1,1)';
beta0_draw=multi_normrnd(ub0, vb0,1)';
sigma1_draw=gamrnd(A1,1);  sigma1_draw=B1/sigma1_draw;
sigma0_draw=gamrnd(A0,1);  sigma0_draw=B0/sigma0_draw;
tp11_draw=0.9; tp00_draw=0.9;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Gibbs-Sampling
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for iter=1:nrep+warmup
    % step 1: draw latent state variables: st
    [st_draw, smooth_t, filter_t]=state_draw_all_uni(y,x,beta1_draw,beta0_draw, sigma1_draw, sigma0_draw, tp11_draw, tp00_draw);

    % step 2: draw the posterior coefficients: beta0 and beta1
    yt_1=y((st_draw==1),:);  yt_0=y((st_draw==0),:);  % sampling beta0 & beta1 we can seperate the data sets to sample them
    xt_1=x((st_draw==1),:);  xt_0=x((st_draw==0),:);
    % sample beta0
    Db0=(xt_0'*xt_0/sigma0_draw+inv(vb0))^-1;
    db0=xt_0'*yt_0/sigma0_draw+inv(vb0)*ub0;
    beta0_draw=multi_normrnd(Db0*db0, Db0,1)';
    % sample beta1
    Db1=(xt_1'*xt_1/sigma1_draw+inv(vb1))^-1;
    db1=xt_1'*yt_1/sigma1_draw+inv(vb1)*ub1;
    beta1_draw=multi_normrnd(Db1*db1, Db1,1)';

    % step 3: draw the posterior volatility variance: sigma0_draw & sigma1_draw, we have constraint that sigma1>sigma0
    num_1=sum(st_draw);
    shape1=A1+num_1/2;
    scale1=B1+sum((yt_1-xt_1*beta1_draw).^2)/2;
    sigma1_draw=gamrnd(shape1,1);
    sigma1_draw=scale1/sigma1_draw;
    
    num_0=T-sum(st_draw);
    shape0=A0+num_0/2;
    scale0=B0+sum((yt_0-xt_0*beta0_draw).^2)/2;
    sigma0_draw=gamrnd(shape0,1);
    sigma0_draw=scale0/sigma0_draw;
    % make sure sigma1_draw>sigma0_draw, since we assume state 1 has high volatility
    j=1;
    while sigma1_draw < sigma0_draw
        sigma0_draw=gamrnd(shape0,1);
        sigma0_draw=scale0/sigma0_draw;
        sigma1_draw=gamrnd(shape1,1);
        sigma1_draw=scale1/sigma1_draw;
        if j>=50
            sigma1_draw=gamrnd(A1,1);
            sigma1_draw=B1/sigma1_draw;
            break
        end
        j=j+1;
    end
    
    % step 4: 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)
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% One can have his own in-sample or out-sample functions to be drawn in each iteration
    %% the following two functions are used to calculate the hedge ratio in-sample and out-sample forecast

    % step 5: draw the in-sample hedge-ratio function
    if nargout > 4
        hedge1=beta1_draw(2); hedge0=beta0_draw(2);         % hedge ratio in regime 1 and regime 0
        hedge_insample=filter_t*[hedge1,hedge0]';           % in-sample hedge ration calculation
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if iter>warmup
        beta1_keep(iter-warmup,:)=beta1_draw';
        beta0_keep(iter-warmup,:)=beta0_draw';
        sigma1_keep(iter-warmup,:)=sigma1_draw;
        sigma0_keep(iter-warmup,:)=sigma0_draw;
        tp11_keep(iter-warmup,:)=tp11_draw;
        tp00_keep(iter-warmup,:)=tp00_draw;
        smpl_st(iter-warmup,:)=st_draw';
        smooth_t_keep(iter-warmup,:)=smooth_t(:,1)';
        filter_t_keep(iter-warmup,:)=filter_t(:,1)';
        if nargout > 4
            in_sample_results(iter-warmup,:)=hedge_insample';
            
        end
    end
    fprintf(' Iteration to %3g \n', iter)
end

smpl_coeff=[beta1_keep, beta0_keep, sigma1_keep, sigma0_keep, tp11_keep, tp00_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 filter probability calculation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [st_draw, smooth_t, filter_t]=state_draw_all_uni(y,x,beta1_draw,beta0_draw, sigma1_draw, sigma0_draw, tp11_draw, tp00_draw)

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
filter_t=zeros(T,2);       % which is the filter proability: filter_t(t,1)=Pr(St=1|Yt)
conditional_t=zeros(T,2);  % conditional_t=Pr(St=k|Yt-1)
smooth_t=zeros(T,2);       % smooth_t=Pr(St=k|Y_T)=Pr(St=k|Y_t)*Pr(St+1|St), which is the product of {filter probability}*{transtion probability}

filter_t(1,:)=[0.5, 0.5];  % assign the initial filter_t value
conditional_t(1,:)=[0.5, 0.5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate the filter probability, then calculate  posterior conditional probability: P(st|Y_T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)>0        % i.e. two likelihood are both zero, avoid 0/0 problem
        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

st_draw(end)=(U(end)<=filter_t(end,1));  % if random uniform number is less than filter_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
        smooth_t(i,:)=filter_t(i,:).*P(:,1)';                % -->[filter_t1*p11, filter_t0* p01]
    elseif st_draw(i+1)==0
        smooth_t(i,:)=filter_t(i,:).*P(:,2)';                % -->[filter_t1*p10, filter_t0* p00]
    end
    smooth_t(i,:)=smooth_t(i,:)./sum(smooth_t(i,:),2);   % Normalize smooth_t
    st_draw(i)=(U(i)<=smooth_t(i,1));
end

smooth_t(end,:)=filter_t(end,:);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Subfunction 2: calculate the state count number:  n11, n10, n01, n00; for inverse garma to generate sigmas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [n11, n10, n01, n00]=sigma_draw_number_uni(st_draw)

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

% Function End