% Script to perform Non-Parametric Hidden Markov Models Estimation.
%
% This routine takes a simulation file from the least squares estimation procedure
% "main_EstimCluster.m" and selects the right model by calibrating the penalty using
% the dimension jump heuristics.
%%

%% Simulation parameters
simulation_file = % data file name %

% Range of constant to check.
rho_min = 0;
rho_max = 5;
rho_step = 0.001;

%% Load data.
load(simulation_file);

% Get selected model for all constants in the specified range.
j = 0;
contrastvalues_redim = contrastvalue(:, min_m:max_m);
base_penalty = ((1:kmax)'*(min_m:max_m) + ((1:kmax).*(1:kmax))'*ones(1,max_m-min_m+1) - 1) * log(n)/n;     % penalty in (mk+k^2-1)log(n)/n
for rho = rho_min:rho_step:rho_max
	j = j+1;
	Penalized_contrast_values = contrastvalues_redim + rho * base_penalty;
	[valmin_m , indmin_m] = min(Penalized_contrast_values, [], 2);
	[~ , indmin_k] = min(valmin_m);
	K_selected(j) = indmin_k;
	M_selected(j) = indmin_m(indmin_k) + min_m - 1;
	M_selectedVect(j,:) = indmin_m + min_m - 1;
end
% Plot selected complexity.
plot(rho_min:rho_step:rho_max, M_selected .* K_selected + K_selected .* (K_selected - 1), 'k-');
ylim([0 (max_m*kmax+kmax*(kmax-1))]);
xlabel('$\rho$', 'Interpreter', 'Latex');
ylabel('Model complexity');

%% Get largest jump and the corresponding constant.
[~,indJump] = min(M_selected(2:end) - M_selected(1:end-1));
rho_jump = rho_min + rho_step*(indJump-1);
% Selected constant.
rho_selected = 2 * rho_jump;
% Get corresponding model.
Penalized_contrast_values = contrastvalues_redim + rho_selected * base_penalty;
[valmin_m , indmin_m] = min(Penalized_contrast_values, [], 2);
[~ , indmin_k] = min(valmin_m);
K_selected = indmin_k;
M_selected = indmin_m(indmin_k) + min_m - 1;

disp('Selected parameters:');
disp(['     M = ' num2str(M_selected)]);
disp(['     K = ' num2str(K_selected)]);
