clear;
clc;
close all;

addpath('data')
addpath('functions')
addpath('figures')
addpath('results')

%% |---------------- Simulations Settings

T           = [50,100,200,500];
DGPs        = 1:8;
MCs         = 500;

%% |---------------- Estimation settings

lags         = 4;
nhor         = 12;
H_max        = 5;            % Number of mixing components 
init_method  = 'cluster';   % cluster initialization method random|cluster|standard
iter         = 20;           % Print every "iter^th" iteration
nthin        = 100;
prior_mixing = 2;           % prior for mixing coefficients: 0: Normal Prior | 1: Tipping's ML prior | 2: Horseshoe Prior
prior_kernel = 0;           % prior for mixing coefficients: 0: Normal Prior | 1: Tipping's ML prior | 2: Horseshoe Prior
basis_exp    = 2;           % 0: no basis function expansion | 1: basis function expansion | 2: panel
max_iter     = 10000;
tol          = 1e-6;        % 1e-10;
R            = 1000;
quants       = [0.05,0.1,0.2,0.5,0.8,0.9,0.95]; 



% prior settings
prior.a_tau = 1;
prior.b_tau = 0.1;

prior.Ab_mix  = 1;
prior.tau_mix = 0.1;     
prior.nu_mix  = 0.1;  

sampler       = 2; 


for dgp = 1:length(DGPs)
    for tt = 1:length(T)
        clusters   = zeros(T(tt),MCs);
        parfor iter = 1:MCs
            if DGPs(dgp) == 1
                y = randn(T(tt),1);
            elseif DGPs(dgp) == 2
                mu    = 1;
                sig   = 1;
                alpha = -.5;
                nu    = 10;
                y     = rskt(T(tt),mu,sig,alpha,nu)';  
            elseif DGPs(dgp) == 3
                mu    = 1;
                sig   = 1;
                alpha = 1;
                nu    = 5;
                y     = rskt(T(tt),mu,sig,alpha,nu)';  
            elseif DGPs(dgp) == 4
                probs = rand(T(tt),1);
                y = (probs<=1/3).*(0.04*randn(T(tt),1)) + (probs>1/3).*(1+0.04*randn(T(tt),1));
            elseif DGPs(dgp) == 5
                probs = rand(T(tt),1);
                y = (probs<=1/6).*(0.2*randn(T(tt),1)) + (probs>1/6 & probs<=1/6+1/3).*(1+0.2*randn(T(tt),1))+ (probs>1/2).*(2+0.2*randn(T(tt),1));
            elseif DGPs(dgp) == 6
                probs = rand(T(tt),1);
                y = (probs<=3/4).*(-43/100+randn(T(tt),1)) + (probs>3/4).*(107/100+1/3*randn(T(tt),1));
            elseif DGPs(dgp) == 7
                probs = rand(T(tt),1);
                y = (probs<=1/10).*(randn(T(tt),1)) + (probs>1/10).*(1/10*randn(T(tt),1));
            elseif DGPs(dgp) == 8
                probs = rand(T(tt),1);
                y = (probs<=2/3).*(randn(T(tt),1)) + (probs>2/3).*(1/10*randn(T(tt),1));    
            end
    
            x = ones(T(tt),1);
            [out] = VB_BDR_pipe(y,x,x,H_max,init_method,prior_mixing,prior_kernel,max_iter,tol,prior);  
            clusters(:,iter)   = out.cluster ;
            
        end
        results.dgp(dgp).T(tt).cluster = clusters;
    end
end


for dgp = 1:length(DGPs)
    for tt = 1:length(T)
        for iter = 1:MCs
            counts(iter,dgp,tt) = length(unique(results.dgp(dgp).T(tt).cluster(:,iter)));
        end
    end
end


save("SIM2_H_select_"+ string(H_max) + "_priorM_"+ string(prior_mixing) + "_priorK_"+ string(prior_kernel) + "_sampler_" + string(sampler) + "_.mat");

nbH = [1 2 3 2 2 2];
hitrate = mean((counts(:,[1,4:8],:)-nbH)==0);
medest  = median(counts(:,[1,4:8],:));

metrics = [];
for i = 1:4
    metrics = [metrics; [hitrate(:,:,i); medest(:,:,i)]];
end

metn = ["Hitrate"; "Median"];
Metrics = [repmat(cellstr(metn),4,1) cellfun(@(x) sprintf('%.2f', x), num2cell(round(metrics,2)), 'UniformOutput', false)];

% Size info
[nRows, nCols, nLayer] = size(Metrics);

% New array: double the number of columns
B = cell(nRows, 2 * nCols, nLayer);

% Fill original columns into odd positions
B(:, 1:2:end, :) = Metrics;

% Fill symbol into even positions
B(:, 2:2:end,:) = repmat({'&'}, nRows, nCols, nLayer);

Metrics = B;
Metrics(:,end,:) = repmat({'"\\"'},nRows,nLayer);