function out = mcmc_func(x_,T_,S_,runs)

%% Unpack Chain Options

N            = S_.Options.EST.N;
burn_in_para = S_.Options.EST.burn_in;
burn_in      = ceil(burn_in_para*N) + 1;

numbersave  = S_.Options.EST.numbersave;         % Number of times chain is saved
display_num = S_.Options.EST.display_num ;
df          = S_.Options.EST.df ;                % Degrees of freedom for multivariate t proposal
kap_        = S_.Options.EST.kap_ ;

dim_x  = length(x_);
L      = S_.Options.EST.upper_xx;
pop    = 1:dim_x; 

param = param_func(x_);

param_names = S_.param_names ;

NumBreaks     = length(T_);
maxT          = S_.Options.EST.maxT ;
minT          = S_.Options.EST.minT ;
popT          = minT:maxT;
upperT        = S_.Options.EST.upperT ;         % maximum number of updated elements in T

if upperT > NumBreaks
    disp('Error: The number of dates breaks to be updated in the proposal cannot be higher than the number of breaks')
    return
end

c = S_.Options.EST.c ;

whensave = linspace(2,N,numbersave);

%% Posterior

S_.Options.EST.kron = 0;           % for mcmc

post_ = @(x_,T_) loglike(x_,T_,S_) + prior(x_);

%% SAMPLER

T_prop = zeros(N,NumBreaks);
T_prop(1,:) = T_;

Prop_xx = zeros(N,length(x_));
Prop_xx(1,:) = x_';

T_out  = zeros(N,NumBreaks);
xx_out = zeros(N,length(x_));
post_out = zeros(N,1);
post_out(1) = S_.Options.EST.post_val;

T_out(1,:) = T_prop(1,:);
xx_out(1,:) = x_';
xx_prop = x_';

kk_T = 0;
kk_xx = 0;
reject_rate = 0;
reject_rate_T = 0;

countersave = 1;

numerator2 = post_(xx_out(1,:),T_out(1,:));
denominator2 = numerator2;
acx = 1;

% Main Loop
for j=2:N

    % T-block
    T_prop(j,:)      = T_out(j-1,:);
    blocksizeT       = ceil(upperT*rand);
    blockT           = randsample(NumBreaks,blocksizeT);
    T_prop(j,blockT) = randsample(popT,blocksizeT);

    if acx == 1
        denominator = numerator2;
    else
        denominator = denominator2;
    end

    numerator        = post_(xx_out(j-1,:),T_prop(j,:));
    ratio            = exp(-numerator + denominator);
    alpha            = min(ratio,1);

    if isnan(ratio)
        alpha = 0;
    end

    if rand <= alpha
        T_out(j,:) = T_prop(j,:);
        kk_T       = kk_T + 1;
        acT        = 1;
        post_out(j,1) = -numerator;
    else
        T_out(j,:) = T_out(j-1,:);
        acT        = 0;
        post_out(j,1) = -denominator;
    end

    acc_T = 100*kk_T/j;

    % xx-block
    xx_prop = xx_out(j-1,:);
    blocksize = ceil(L*rand);
    block = randsample(pop,blocksize);
    xx_prop(1,block) = xx_out(j-1,block)+ (c(block,block)*trnd(df,blocksize,1))'; % multivariate-t with df degrees of freedom
    Prop_xx(j,:) = xx_prop(1,:);

    if min(xx_prop([param.order.h,param.order.rho_g:param.order.rho_b])) <= 0 || max(xx_prop([param.order.h,param.order.rho_g:param.order.rho_b])) >= 1 || xx_prop(param.order.psi_b)<=0
        alpha2 = 0;
        ratio2 = 0;
        reject_rate = reject_rate + 1;
    else
        if acT == 1
            denominator2 = numerator;
        else
            denominator2 = denominator;
        end

        numerator2    = post_(xx_prop,T_out(j,:));
        ratio2        = exp(-numerator2 + denominator2);
        alpha2        = min(ratio2,1);
    end

    if isnan(ratio2)
        alpha2 = 0;
    end

    if rand <= alpha2
        xx_out(j,:) = xx_prop;
        kk_xx       = kk_xx + 1;
        acx         = 1;
        post_out(j,1) = -numerator2;
    else
        xx_out(j,:) = xx_out(j-1,:);
        acx         = 0;
        post_out(j,1) = -denominator2;
    end

    acc_xx= 100*kk_xx/j;


    if mod(j/display_num,1)==0
        fprintf(sprintf('Chain %1.0f.... %2.2f per cent done\n      ***** accept rate: x, %2.3f and T, %2.3f \n',...
            runs,100*j/N,acc_xx,acc_T)) ;
        % disp(mean(xx_out(1:j,:)));
    end

    if mod(j/1000,1)==0
        header = {'Parameter','Mean','5','95'} ;
        Tabout1 = table(...
            [param_names(1:27)], ...
            round([(mean(xx_out(1:j,1:27)))'],4), ...
            round([(prctile(xx_out(1:j,1:27),5))'],4), ...
            round([(prctile(xx_out(1:j,1:27),95))'],4), ...
            'VariableNames',header)
        Tabout2 = table(...
            [param_names(28:end)], ...
            round([(mean(xx_out(1:j,28:end)))'],4), ...
            round([(prctile(xx_out(1:j,28:end),5))'],4), ...
            round([(prctile(xx_out(1:j,28:end),95))'],4), ...
            'VariableNames',header)
    end


    if not(isempty(whensave))
        if j ==floor(whensave(1,countersave))
            if j == 2
                from = 1;
            else
                from = floor(whensave(1,countersave-1))+1;
            end
            eval(['save mhall_' num2str(countersave),'_'])
            eval(['delete mhall_' num2str(countersave-1),'_.mat'])
            countersave = countersave+1;
        end
    end
end

eval(['delete mhall_' num2str(countersave),'_.mat'])

%% Outputs

out.xx_out   = xx_out ;
out.T_out    = T_out ;
out.post_out = post_out;
out.acc_xx   = acc_xx ;
out.acc_T    = acc_T ;
out.Prop_xx  = Prop_xx;
out.T_prop   = T_prop ; 


