%
%    The author of this script is Danilo Aringhieri.
%
clear all
License=char({'','This software is released by Danilo Aringhieri, who is its sole author. ',...
' Everyone is warned against making use of the specific content of this script, ',...
' wholly or partially, for any scientific work ',...
' without including Danilo Aringhieri as a coauthor. '});
disp(License)
disp('    Press the "Enter" button on the keyboard to START...')
response=input(' ','s');
while isempty(response)==0
disp('    Please, press the "Enter" button to START!')
response=input(' ','s');
end
Intro1=char({'', ' A spreadsheet, able to read and save xls-type or xlsx-type files, ',...
' must be installed on your computer and the variable values, to be imported, ',...
' must be tabled by columns and rows. ',...
' Opening columns on the left may be filled with some lists of strings ',...
' (firms'' names in alphabetic characters, alphanumeric identity codes or addresses ) ',...
' as well as opening rows on the top may be filled with the dates of observations ',...
' or any other details written in string characters, without compromising the import. ',...
'                        DO YOU WANT TO CONTINUE ? ',...
' (Answer "Y" for ''Yes'' or "N" for ''Not'', then press [Enter]) '});
disp(Intro1)
response=input('  ','s');
while strcmp(response,{'Y','y','N','n'})==[0,0,0,0]
disp('   (You must answer exclusively "Y" or "N" !) ')
response=input('  ','s');
end
if sum(strcmp(response,{'N','n'}))==1
disp(char({' ',' OK: you are going to return directly ',...
' to the prompt state in the Command Window! ',' '}))
clear Intro1 response
return
elseif sum(strcmp(response,{'Y','y'}))==1
clear Intro1
end
Input0=char({'',' (You can STOP at any time the program ',...
'  typing  "Ctrl" + "C"  on the keyboard !) '});
disp(Input0)
Input1=char({'','  To import automatically your data, ',...
'       WRITE THE NAME OF YOUR FILE INCLUDING, ',...
'       its extension (".xls" or ".xlsx"), ',...
'       then press [Enter]. ',...
'    (any other command will provide an error message as a result; ',...
'     if it should appear a warning message, please wait few seconds. '});
disp(Input1)
filename=input('   ','s');
[AuxData1,Titles]=xlsread(filename);
dim_AuxData1=size(AuxData1);
for k=1:dim_AuxData1(2)
nan_col(k)=(sum(isnan(AuxData1(:,k))))./dim_AuxData1(1);
end
% If the k-th column in AuxData1 is made only by 'NaN' elements,
% it is nan_col(k)=1 otherwise, if there exist some finite numbers
% in the same column, it holds nan_col(k)<1.
not_nan_col=find(nan_col-1);
% The elements of the row vector 'not_nan_col' identify the positions,
% in the matrix AuxData1, of all the columns not completely filled
% with 'NaN-type' elements.
dim_notnan_col=size(not_nan_col);
AuxData2=AuxData1(:,[not_nan_col]);
% All the columns in the original datasheet (and thus in the matrix AuxData1, too),
% wholly filled with 'NaN-type' elementshas been now cancelled from AuxData1:
% the result is the matrix Auxdata2.
NumNotData1=char({'',' This tool tolerates that, in the original datasheet, ',...
' some columns are filled with numbers which ARE NOT NUMERICAL DATA, ',...
' like when the statistical units are indentified (also) by ordinals. ',...
'       ALL these kind of numerical non-data columns ',...
'       HAVE TO STAY AT THE LEFT of the variable values !! ',...
' Do some numerical non-data columns exist in the original datasheet? ',...
'     ( Answer "Y" if they exist or "N" otherwise, then press [Enter] ) '});
disp(NumNotData1)
response=input('   ','s');
while strcmp(response,{'Y','y','N','n'})==[0,0,0,0]
disp('   (You must answer exclusively "Y" or "N" !) ')
response=input(' ','s');
end
if sum(strcmp(response,{'N','n'}))==1
All_Data=AuxData2;
elseif sum(strcmp(response,{'Y','y'}))==1
NumNotData2=char({' How many are the numerical non-data columns in the original sheet ? '});
disp(NumNotData2)
nondata_col=input('   ');
if (isscalar(nondata_col)==0 || isreal(nondata_col)==0 || isnan(nondata_col)==1 || isinf(nondata_col)==1 ||...
        nondata_col<0 || floor(nondata_col)~=nondata_col)
while (isscalar(nondata_col)==0 || isreal(nondata_col)==0 || isnan(nondata_col)==1 || isinf(nondata_col)==1 ||...
        nondata_col<0 || floor(nondata_col)~=nondata_col)
if isscalar(nondata_col)==0
    disp(' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ')
    disp('   BUT  you must insert exclusively a positive integer ! ) ')
    nondata_col=input('   ');
elseif (isscalar(nondata_col)==1 && (isnan(nondata_col)==1 || isinf(nondata_col)==1 || isreal(nondata_col)==0))
    disp(' Your input is a NaN value or infinite or a complex number: WHY? ')
    nondata_col=input('   ');
elseif nondata_col<0
    disp(' ( You must insert exclusively a positive number ! ) ')
    nondata_col=input('   ');
elseif (nondata_col>0 && floor(nondata_col)~=nondata_col)
    disp(' ( You must insert exclusively a positive integer ! ) ')
    nondata_col=input('   ');
end
end
elseif (isscalar(nondata_col)==1 && isreal(nondata_col)==1 && isnan(nondata_col)==0 && isinf(nondata_col)==0 ...
    && nondata_col>0 && floor(nondata_col)==nondata_col)
end
All_Data=AuxData2(:,nondata_col+1:dim_notnan_col(2));
end
% Now all the columns wholly filled with numerical values, which are not values
% for the variable, are cancelled from Auxdata2: the result is the matrix
% All_Data, only containing either the available variable's values
% for every statistical unit or a 'NaN' element, where and where the variable is unknown.
% The number of columns in All_Data is the number of all observation periods
% recorded in the sheet of the variable.
clear NumNotData1 NumNotData2
dim_All_Data=size(All_Data);
% Now the tool is going to extract, from the matrix 'All_Data', the subsamples
% of the incumbents, that is all the statistical units in the datasheet expressing
% a value for everyone of the observation periods registered, for the variable of interest.
for k=1:dim_All_Data(1)
incumbent_row(k)=prod(isfinite(All_Data(k,:)));
end
% The k-th individual in the spreadsheet is an incumbent if and only if
% isfinite(All_Data(k,:) is a vector whose elements are all equal to 1
% so the component 'incumbent_row(k)= prod(isfinite(All_Data(k,:)))'
% is 1 for an incumbent but it is zero for any agent lacking of a datum for some periods.
incumbent_posit=find(incumbent_row);
% The elements of the row vector 'incumbent_posit' identify the positions,
% in the matrix All_Data, of all agents (one agent corresponds to one and only one row)
% which are incumbents. The element inside All_Data which are not incumbent
% are going to be excluded from the workspace and the final result
% is 'incumbData', the cleaned matrix of the only incumbent agents.
incumbData=All_Data([incumbent_posit],:);
dim_incumb_pos=size(incumbent_posit);
num_of_incumb=dim_incumb_pos(2);
share_of_incumb=num_of_incumb/dim_All_Data(1);
dim_incumb=size(incumbData);
disp(char({' ','   The number of all individuals ',...
'   in the original spreadsheet '}))
fprintf('   is %-1.0d \n',dim_All_Data(1))
disp('   but number of the extracted incumbents  ')
fprintf('   is %-1.0d statistical units. \n',num_of_incumb)
MaxValues=max(incumbData);
MinValues=min(incumbData);
disp(char({' ','  What is the dimensional unit of your data? ',...
'  (Use the expression you like and press [Enter]) '}))
units=input('  ','s');
disp(char({' ','  Insert the time unit in which your observations are marked, ',...
'  then press [Enter] '}))
time_scan=input('  ','s');
Input2_1=char({'',' It is assumed observation dates are in chronological order: ',...
' so they will be implicitly marked, according to their relative position ',...
' in the datasheet, by progressive positive integers assigned to the auxiliary variable ',...
' "t_start" of the time when a transition starts. ',...
'',' All the numerical data from your file ',...
' will be recorded into the matrix   "incumbData"   in the workspace. ',...
' The opening rows and columns in a string format will be copied ',...
' into the array "Titles" which will be cancelled soon after. ',...
' The ''NaN'' elements inside "incumbData" correspond, inside the spreadsheet, ',...
' to empty cells or to cells where no values for the variable are available. ',...
'         (press [Enter] to go on...) '});
disp(Input2_1)
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to go on...)')
response=input(' ','s');
end
Input2_2=char({'',' Three vectors will be created: ',...
' "dim_incumb", "MaxValues" and  "MinValues" . ',...
'',' "dim_incumb" consists of the dimensions of incumbData:',...
'  dim_incumb(1) is the number of incumbents enumerated inside the datasheet. ',...
'  dim_incumb(2) is the number of periods when the variable''s values have been measured. ',...
'','  " MaxValues "  and  " MinValues " are row vectors: ',...
'  " MaxValues(m) "  and  " MinValues(m) "  are respectively ',...
'  the biggest datum and the smallest one registered in the m-th column of "incumbData"  ',...
'  referring to the m-th observation period. ',...
'         (press [Enter] to see) ','',});
disp(Input2_2)
clear resp
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to see)')
response=input(' ','s');
end
disp('  According to your sample, the vector  " MinValues "  is: ')
if dim_incumb(2)==2
fprintf('     [%-1.3f, ',MinValues(1));...
fprintf('%-1.3f] \n',MinValues(2))
else
fprintf('     [%-1.3f, ',MinValues(1));...
fprintf('%-1.3f, ',MinValues(2:dim_incumb(2)-1));...
fprintf('%-1.3f] \n',MinValues(dim_incumb(2)))
end
disp('  and the vector  " MaxValues "  is: ')
if dim_incumb(2)==2
fprintf('     [%-1.3f, ',MaxValues(1));...
fprintf('%-1.3f] \n',MaxValues(2))
else
fprintf('     [%-1.3f, ',MaxValues(1));...
fprintf('%-1.3f, ',MaxValues(2:dim_incumb(2)-1));...
fprintf('%-1.3f] \n',MaxValues(dim_incumb(2)))
end
clear Titles
lowest_datum=min(MinValues');
highest_datum=max(MaxValues');
fprintf('   (all the values are ''%s'') \n',units)
disp(' ')
disp(' In your whole spreadsheet, absolute minimum and maximum values are: ')
fprintf('   Minimal datum = %-1.3f \n',lowest_datum)
fprintf('   Maximal datum = %-1.3f \n',highest_datum)
fprintf('                  (%s) \n',units)
disp(char({' Knowing such extremes is helpful to decide the partition ',...
' on which the states, available for the transitions of agents, ',...
' will be set. ',...
'                 (Press [Enter] to continue...) '}))
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to continue...)')
response=input(' ','s');
end
clear AuxData1 AuxData2 dim_AuxData1 dim_incumb_pos
disp(char({' Are you interested in knowing also ',...
' the set of quantiles of a certain degree, ',...
' along with the median and the mean, of your imported data? ',...
' You might adopt the quantiles to create the partition ',...
' on which the states for the transition matrix are based. ',...
' ("Y" for ''Yes'' or "N" for ''Not'', then [Enter]) '}))
response=input('   ','s');
while strcmp(response,{'Y','y','N','n'})==[0,0,0,0]
disp('   (You must answer exclusively "Y" or "N" !) ')
response=input(' ','s');
end
if sum(strcmp(response,{'Y','y'}))==1
disp(char({' Insert the degree of quantiles you are interested in ',...
' (e.g. ''4'' for the quartiles, ''5'' for the quintiles, ',...
' ''10'' for the deciles, ''100'' for centiles, and so on). '}));
qntl_degree=input('   ');
if (isscalar(qntl_degree)==0 || isreal(qntl_degree)==0 || isnan(qntl_degree)==1 || isinf(qntl_degree)==1 ||...
        qntl_degree<=0 || floor(qntl_degree)~=qntl_degree)
while (isscalar(qntl_degree)==0 || isreal(qntl_degree)==0 || isnan(qntl_degree)==1 || isinf(qntl_degree)==1 ||...
        qntl_degree<=0 || floor(qntl_degree)~=qntl_degree)
if isscalar(qntl_degree)==0
    disp(' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ')
    disp('   BUT  you must insert exclusively a positive integer ! ) ')
    qntl_degree=input('   ');
elseif (isscalar(qntl_degree)==1 && (isnan(qntl_degree)==1 || isinf(qntl_degree)==1 || isreal(qntl_degree)==0))
    disp(' Your input is a NaN value or infinite or a complex number: WHY? ')
    qntl_degree=input('   ');
elseif qntl_degree<=0
    disp(' ( You must insert exclusively a positive number ! ) ')
    qntl_degree=input('   ');
elseif (qntl_degree>0 && floor(qntl_degree)~=qntl_degree)
    disp(' ( You must insert exclusively a positive integer ! ) ')
    qntl_degree=input('   ');
end
end
elseif (isscalar(qntl_degree)==1 && isreal(qntl_degree)==1 && isnan(qntl_degree)==0 && isinf(qntl_degree)==0 ...
        && qntl_degree>0 && floor(qntl_degree)==qntl_degree)
end
q=zeros(1,qntl_degree);
for k=1:qntl_degree
q(k)=k/qntl_degree;
end
Quantiles=(quantile(incumbData,q,1))';
Means=mean(incumbData);
Medians=median(incumbData);
Modes=mode(incumbData);
disp(char({'',' The matrix of quantiles you have chosen is named " Quantiles " ',...
' and it is shown below: every observation period in the datasheet corresponds orderly ',...
' to a row in "Quantiles" and its columns from the second to the last one ',...
' corresponds to the quantile ordering for such a date. ','',' Quantiles = '}))
disp(Quantiles)
disp(char({' The means, medians and modes of the data for each period, ',...
' instead, are ordered in three row vectors, ',...
' "Means", "Medians" and "Modes" respectively: ',...
'',' Means = '}))
disp(Means)
disp(' Medians = ')
disp(Medians)
disp(' Modes = ')
disp(Modes)
disp('    Press [Enter] to continue...  ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to continue...)')
response=input(' ','s');
end
elseif sum(strcmp(response,{'N','n'}))==1
end
Input3=char({'',' ALL AGENTS are going to be GROUPED ',...
' IN "n" CLASSES OF NOT-PREDETERMINED WIDTH, ',...
' splitting the whole spectrum of the variable in the following manner: ',...
'','[0;sup(1)), [sup(1);sup(2)), [sup(2);sup(3)),.., [sup(n-1);infinity) ',...
' where sup(1) is bigger than the minimal datum in the sheet',...
' and sup(n-1) is smallest than the maximal datum. ',...
'',' Each class represents one of the states among which a transition happens ',...
' and will be identified by a progressive positive integer. ',...
'',' INSERT THE NUMBER OF CLASSES ',...
' and press  [Enter]: '});
disp(Input3)
classes_num=input('   ');
while (isscalar(classes_num)==0 || isreal(classes_num)==0 || isnan(classes_num)==1 || isinf(classes_num)==1 ...
        || classes_num<=0 || floor(classes_num)~=classes_num)
if isscalar(classes_num)==0
disp(char({' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ',...
'   BUT  you must insert exclusively a positive integer ! ) '}))
classes_num=input('   ');
elseif (isscalar(classes_num)==1 && (isnan(classes_num)==1 || isinf(classes_num)==1 || isreal(classes_num)==0))
disp(' Your input is a NaN value or a not finite one or it is not a real number: WHY? ')
classes_num=input('   ');
elseif classes_num<=0
disp(' ( You must insert exclusively a positive number ! ) ')
classes_num=input('   ');
elseif (classes_num>0 && floor(classes_num)~=classes_num)
disp(' ( You must insert exclusively a positive integer ! ) ')
classes_num=input('   ');
end
end
disp(' ')
fprintf('   Now write the sup for everyone of the first %-1.0d classes, \n',classes_num-1)
disp(char({'   which are the inner points of your partition ',...
'   of the variable spectrum [ Minimal datum, Maximal  datum ]. '}))
for k=1:classes_num-1
    fprintf('Insert the top (expressed as number of ''%s'') of the %-1.0d° class:\n',units,k)
    sup_k=input(' ');
    if k==1
        if not(isscalar(sup_k)==1 && isreal(sup_k)==1 && isnan(sup_k)==0 && isinf(sup_k)==0 &&...
                sup_k>0 && sup_k>lowest_datum && sup_k<highest_datum)
            while not(isscalar(sup_k)==1 && isreal(sup_k)==1 && isnan(sup_k)==0 && isinf(sup_k)==0 &&...
                    sup_k>0 && sup_k>lowest_datum && sup_k<highest_datum)
                if ( isscalar(sup_k)==0)
                    disp(' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ')
                    disp('   BUT  you must insert exclusively a positive integer ! ) ')
                    sup_k=input('   ');
                elseif (isscalar(sup_k)==1 && (isnan(sup_k)==1 || isinf(sup_k)==1 || isreal(sup_k)==0))
                    disp(' Your input is a NaN value or infinite or a complex number: WHY? ')
                    sup_k=input('   ');
                elseif sup_k<=0
                    disp(' ( You cannot insert a negative number or zero ! ) ')
                    sup_k=input('   ');
                elseif sup_k<=lowest_datum
                    disp(char({' Insert a value bigger than Minimal datum (but smaller than Maximal datum): '}))
                    sup_k=input(' ');
                elseif sup_k>=highest_datum
                    disp(char({' Insert a value smaller than Maximal datum (and bigger than Minimal datum): '}))
                    sup_k=input(' ');
                end
            end
        else
        end
        sup(1)=sup_k;
    elseif k>1
        if not(isscalar(sup_k)==1 && isreal(sup_k)==1 && isnan(sup_k)==0 && isinf(sup_k)==0 &&...
                sup_k>0 && sup_k>lowest_datum && sup_k<highest_datum && sup(k-1)<sup_k)
            while not(isscalar(sup_k)==1 && isreal(sup_k)==1 && isnan(sup_k)==0 && isinf(sup_k)==0 &&...
                    sup_k>0 && sup_k>lowest_datum && sup_k<highest_datum && sup(k-1)<sup_k)
                if ( isscalar(sup_k)==0)
                    disp(' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ')
                    disp('   BUT  you must insert exclusively a positive integer ! ) ')
                    sup_k=input('   ');
                elseif (isscalar(sup_k)==1 && (isnan(sup_k)==1 || isinf(sup_k)==1 || isreal(sup_k)==0))
                    disp(' Your input is a NaN value or infinite or a complex number: WHY? ')
                    sup_k=input('   ');
                elseif sup_k<=0
                    disp(' ( You cannot insert a negative number or zero ! ) ')
                    sup_k=input('   ');
                elseif sup_k<=lowest_datum
                    disp(char({' Insert a value bigger than Minimal datum (but smaller than Maximal datum): '}))
                    sup_k=input(' ');
                elseif sup_k>=highest_datum
                    disp(char({' Insert a value smaller than Maximal datum (and bigger than Minimal datum): '}))
                    sup_k=input(' ');
                elseif  sup_k<=sup(k-1)
                    disp(char({' Insert a value greater than the previous one '}))
                    sup_k=input(' ');
                end
            end
        else
        end
        sup(k)=sup_k;
    end
end
clear sup_k
partition=[0,sup,Inf];
Input4=char({'',' The extremes of the partition you have chosen ',...
' are arranged orderly as components inside the vector "partition". ',''});
disp(Input4)
clear Input1 Input2_1 Input2_2 Input3 Input4
disp('   ( Press [Enter] to continue... ) ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to continue...)')
response=input(' ','s');
end
%
%   From this point, forward to the creation of the array "mov_k_tstart",
%   we are referrring to the statistical inference about first order
%   Markov chains presented in Sections 2.1 and 2.2 from the 1956 paper 
%   by Anderson & Goodman, "Statistical inference about Markov Chains".
%
%   n_k_tstart     is the number of agents belonging to the k-th class
%                  at the starting period (t_start), ready to make a one-period transition.
%
%   sh_k_tstart is the percentage of units belonging to the k-th class at t_start.
%                  For such array, as well as for "n_k_tstart", t_start is an integer  
%                  varying from 1 to dim_incumb(2), differently from the other arrays
%                  which count agents or transitions
%
%   num_k_h_tstart is the number of agents making a transition
%                  from the state 'k' occupied at time t_start
%                  to the state 'h' occupied at time t_start+1
%
%   mov_k_tstart   is the percentage of moving agents, which were inside the k-th class
%                  at the starting period t_start, observed inside any other class, 
%                  different from the k-th one, at the successive period t_start+1
%
%   nstar_k_tult   is the cumulated number of all one-period transitions
%                  starting from the k-th class and occurring between t_start=1
%                  and t_start=t_ult where (t_ult>=1 and t_ult<=dim_incumb(2)-1
%
%   n_k_h_tult     is the cumulated number of one-period transitions occurring
%                  from state 'k' at time t_start=1 to state 'h'
%                  at every subsequent period t_start=t_ult<=dim_incumb(2)-1
%
%   p_stat_k_h     is the maximum likelihood estimation of the first order
%                  stationary transition probability from the k-th class
%                  to the h-th one, computed between t_start=1 and t_start=dim_incumb(2)-1
%
%   p_nonstat_k_h_tstart  is the maximum likelihood estimation of the first order,
%                         not-stationary transition probability from the k-th state
%                         to the h-th one, computed between t_start and t_start+1
%
% PAY ATTENTION: differently from the other counting arrays,
% "n_k_tstart" is defined for everyone of the observation dates 
% in the spreadsheet, from t_start=1 to t_start=dim_incumb(2),
% instead of till to dim_incumb(2)-1 only.
%
n_k_tstart=zeros(classes_num,dim_incumb(2));
num_k_h_tstart=zeros(classes_num,classes_num,dim_incumb(2)-1);
for t_start=1:dim_incumb(2)-1
for k=1:classes_num
for h=1:classes_num
num_k_h_tstart(k,h,t_start)=sum((incumbData(:,t_start)>=partition(k) & incumbData(:,t_start)<partition(k+1))&...
    (incumbData(:,t_start+1)>=partition(h) & incumbData(:,t_start+1)<partition(h+1)));
end
end
end
for t_start=1:dim_incumb(2)
for k=1:classes_num
n_k_tstart(k,t_start)=sum(incumbData(:,t_start)>=partition(k) & incumbData(:,t_start)<partition(k+1));
sh_k_tstart(k,t_start)=n_k_tstart(k,t_start)./dim_incumb(1);
end
end
p_nonstat_k_h_tstart=zeros(classes_num,classes_num,dim_incumb(2)-1);
for t_start=1:dim_incumb(2)-1
for k=1:classes_num
for h=1:classes_num
if n_k_tstart(k,t_start)==0
    continue
end
p_nonstat_k_h_tstart(k,h,t_start)=num_k_h_tstart(k,h,t_start)./n_k_tstart(k,t_start);
end
end
end
n_k_h_tult=zeros(classes_num,classes_num,dim_incumb(2)-1);
nstar_k_tult=zeros(classes_num,dim_incumb(2)-1);
for k=1:classes_num
for t_ult=1:dim_incumb(2)-1
for h=1:classes_num
n_k_h_tult(k,h,t_ult)=sum(num_k_h_tstart(k,h,1:t_ult));
end
end
end
for k=1:classes_num
for t_ult=1:dim_incumb(2)-1
nstar_k_tult(k,t_ult)=sum(n_k_h_tult(k,:,t_ult));
end
end
p_stat_k_h=zeros(classes_num,classes_num);
for k=1:classes_num
for h=1:classes_num
if nstar_k_tult(k,dim_incumb(2)-1)==0
continue
end
p_stat_k_h(k,h)=n_k_h_tult(k,h,dim_incumb(2)-1)./nstar_k_tult(k,dim_incumb(2)-1);
end
end
p_up_nonstat_k_h_tstart=zeros(classes_num,classes_num,dim_incumb(2)-1);
p_dn_nonstat_k_h_tstart=zeros(classes_num,classes_num,dim_incumb(2)-1);
for t_start=1:dim_incumb(2)-1
p_up_nonstat_k_h_tstart(:,:,t_start)=triu(p_nonstat_k_h_tstart(:,:,t_start),1);
p_dn_nonstat_k_h_tstart(:,:,t_start)=tril(p_nonstat_k_h_tstart(:,:,t_start),-1);
end
p_up_stat_k_h=triu(p_stat_k_h,1);
p_dn_stat_k_h=tril(p_stat_k_h,-1);
%
% Now we are going to reckon the products of both the stationary transition matrix
% and the not-stationary ones between the first period of observations in the datasheet
% and everyone of the subsequent periods until the last one, in order to calculate
% also the aggregate values of the mobility indices of our intererest.
%
%    power_p_stat_k_h_tult    is the compound stationary transition matrix
%                             from the first period t_start=1 to each of the following 
%                             periods t_ult=1,2,.., dim_incumb(2)-1
%
%    aggr_p_nonstat_k_h_tult  is the compound not-stationary transition matrix
%                             from the first period t_start=1 to everyone  
%                             of the following t_ult=2,3,..,dim_incumb(2)-1.
%                             
for t_ult=1:dim_incumb(2)-1
power_p_stat_k_h_tult(:,:,t_ult)=(p_stat_k_h)^(t_ult);   
end
for t_start=1:dim_incumb(2)-1
prdct(:,:,t_start)=eye(classes_num,classes_num);
end
aggr_p_nonstat_k_h_tult=zeros(classes_num,classes_num,dim_incumb(2)-1);
for t_ult=1:dim_incumb(2)-1
if t_ult==1 
aggr_p_nonstat_k_h_tult(:,:,t_ult)=p_nonstat_k_h_tstart(:,:,1);
continue
end
for t_start=1:t_ult
prdct(:,:,t_start+1)=prdct(:,:,t_start)*p_nonstat_k_h_tstart(:,:,t_start);
end
aggr_p_nonstat_k_h_tult(:,:,t_ult)=prdct(:,:,t_start+1);
end
clear prdct
%
% Chi-Square Test and Likelihood Ratio Criterion to test the null hypothesis
% that the one-period transition probabilities from a first order Markov chain 
% are constant respect to the time (Anderson & Goodman 1956, Sect. 3.2).
%
disp(char({'',' In this part we are going to check up the 1st null hypothesis ',...
' that the one-period transition probabilities from a first order Markov chain, ',... 
' based on the imported data, are constant respect to the time, ',...
' applying the Chi-Square Test and the Maximum Likelihood Ratio Criterion (MLRC) ',...
' as proposed in Section 3.2 of the 1956 paper by Anderson & Goodman. ',''}))
r_X_timedepend=zeros(classes_num,classes_num,dim_incumb(2)-1);
r_MLR_timedepend=zeros(classes_num,classes_num,dim_incumb(2)-1);
for k=1:classes_num
for h=1:classes_num
for t_start=1:dim_incumb(2)-1
if p_stat_k_h(k,h)==0
   continue 
end
r_X_timedepend(k,h,t_start)=((p_nonstat_k_h_tstart(k,h,t_start)-p_stat_k_h(k,h)).^2)./p_stat_k_h(k,h);
end
end
end
for k=1:classes_num
for h=1:classes_num
for t_start=1:dim_incumb(2)-1
if p_stat_k_h(k,h)==0
r_MLR_timedepend(k,h,t_start)=1;  
else
r_MLR_timedepend(k,h,t_start)=(p_stat_k_h(k,h)./p_nonstat_k_h_tstart(k,h,t_start)).^num_k_h_tstart(k,h,t_start);
end
end
end
end
% "df" is the number of degree of freedom in both the Chi-Square Test
% and the Likelihood Ratio Criterion.
a_k=zeros(1,classes_num);
b_k=zeros(1,classes_num);
for k=1:classes_num
a_k(k)=length(find(p_stat_k_h(k,:)));
%  By construction, "n_k_tstart" is the number of individuals staying in 'k'
%  for each observation date; the number of observation dates is "dim_incumb(2)"
%  BUT we are interested in the number of transition periods
%  with a positive number of individuals concluding a transition. 
%  From the last observation date in the database, no transitions start !!.
%  So, for the degrees of freedom, we are not interested in "n_k_tstart(:,dim_incumb(2))."
%  We have to count all the positive values of "n_k_tstart"
%  only from the first date to the penultimate one, marked as "dim_incumb(2)-1",
%  because the number of periods in the datasheet is just "dim_incumb(2)-1".
b_k(k)=length(find(n_k_tstart(k,1:dim_incumb(2)-1)));
end
df_timedepend_k(:)=(a_k(:)-1).*(b_k(:)-1);
df_timedepend_tot=(a_k-1)*(b_k-1)';
%
% "chisq_timedepend_oss_k" is a column vector, relative to the specific imported data matrix,
% collecting, for each starting state 'k', the observed value of the Chi-Square statistics
% representing the null hypothesis respect to the alternative.
% "chisq_timedepend_oss" is the scalar of the total observed statistics 
% for the null hypothesis against the alternative, without distinguishing
% the different starting states.
% Both of them admit as asymptotic distributon the Chi-Square one. 
%
chisq_timedepend_oss_k=sum(n_k_tstart(:,1:dim_incumb(2)-1).*squeeze(sum(r_X_timedepend,2)),2);
chisq_timedepend_oss=sum(chisq_timedepend_oss_k);
%
% "p_X2_timedepend_k" and "p_X2_timedepend_tot" are the probabilities that the value of a chi-square variable 
% may be greater than the observed chi-square statistics, for the hypothesis of 1st order Markovian
% stationary transition probabilities, respectively in the cases of distinguished starting states 'k'
% and without any distinction.
%
for k=1:classes_num
p_X2_timedepend_k(k)=chi2cdf(chisq_timedepend_oss_k(k),df_timedepend_k(k),'upper');
end
disp(char({' According to the Chi Square Test, the probability that, ',...
' depending on every starting state ''k'', Markovian first-order ',...
' transition probabilities, based on the imported data, are constant in time, ',...
' is the k-th component in the vector : ','',' p_X2_timedepend_k(k) = '}))
disp(p_X2_timedepend_k)
p_X2_timedepend_tot=chi2cdf(chisq_timedepend_oss,df_timedepend_tot,'upper');
disp(char({' and the probability that the whole 1st order Markov chain, ',...
' based on the data, is stationary, is: ','',' p_X2_timedepend_tot = '}))
disp(p_X2_timedepend_tot)
if chisq_timedepend_oss<=df_timedepend_tot
disp(char({' The value of the Chi^2 statistics computed onto the data for the tested sample ',...
' isn''t greater than number of degree of freedom "df_timedepend_tot", so we can conclude ',...
' that the 1st null hypothesis is confirmed. '}))
elseif chisq_timedepend_oss>df_timedepend_tot
if p_X2_timedepend_tot>=5/100
disp(char({' The value of the Chi^2 statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_timedepend_tot", but the probability ',...
' "p_X2_timedepend_tot" for such an event isn''t lower than 5%, ',...
' so we can accept the 1st null hypothesis at the 5% significance level. '}))    
elseif p_X2_timedepend_tot<5/100
disp(char({' The value of the Chi^2 statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_timedepend_tot" and the probability ',...
' "p_X2_timedepend_tot" for such an event is lower than 5%, ',...
' so we must reject the 1st null hypothesis at the 5% significance level. '}))
end
end
disp('   ( Press [Enter] to continue... ) ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to continue...)')
response=input(' ','s');
end
%
% "MLRcrit_timedepend_k" is a column vector, relative to the imported data, collecting, 
% for each starting state 'k', the opposite of the double logarithm
% of the observed value of the Maximum Likelihood Ratio 
% for the null hypothesis against its alternative.
% "MLRcrit_timedepend_tot" is the scalar of the total observed statistics
% of the Maximum Likelihood Ratio Criterion (MLRC)
% for the null hypothesis against the alternative, without distinguishing
% the different starting states.
% Both of them admit, as an asymptotic distributon, the Chi-Square one. 
%
lambda_timedepend_oss_k=prod(squeeze(prod(r_MLR_timedepend,2)),2);
MLRcrit_timedepend_k=-2*log(lambda_timedepend_oss_k);
MLRcrit_timedepend_tot=sum(MLRcrit_timedepend_k);
disp(char({' According to the Maximum Likelihood Ratio Criterion, ',...
' the probability, depending on each starting state ''k'', that Markovian first-order ',...
' transition probabilities, from the imported data, are time-independent, ',...
' is the k-th component of the vector : ','',' p_MLR_timedepend_k(k) = '}))
%
% The degrees of freedom for the Maximum Likelihood Ratio Criterion 
% are the same, "nu_X_k" and "nu_X_tot", as for the previous Chi-Square Test.
%
% "p_MLR_timedepend_k" and "p_MLR_timedepend_tot" are the probabilities that the value of a chi-square variable 
% may be greater than the observed MLRC statistics, for the hypothesis of 1st order Markovian
% stationary transition probabilities, respectively by distinguishing, or not, 
% the starting states 'k'.
%
for k=1:classes_num
p_MLR_timedepend_k(k)=chi2cdf(MLRcrit_timedepend_k(k),df_timedepend_k(k),'upper');
end
disp(p_MLR_timedepend_k)
p_MLR_timedepend_tot=chi2cdf(MLRcrit_timedepend_tot,df_timedepend_tot,'upper');
disp(char({' and the probability that the whole 1st order Markov chain, ',...
' built up on the data, is stationary, is: ','',' p_MLR_timedepend_tot = '}))
disp(p_MLR_timedepend_tot)
if MLRcrit_timedepend_tot<=df_timedepend_tot
disp(char({' The value of the MLRC statistics computed onto the data for the tested sample ',...
' isn''t greater than number of degree of freedom "df_timedepend_tot", so we can conclude ',...
' that the 1st null hypothesis is confirmed by the MLR Criterion. ',''}))
elseif MLRcrit_timedepend_tot>df_timedepend_tot
if p_MLR_timedepend_tot>=5/100
disp(char({' The value of the MLRC statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_timedepend_tot", but the probability ',...
' "p_MLR_timedepend_tot" for such an event isn''t lower than 5%, ',...
' so the 1st null hypothesis is accepted by the MLR Criterion ',...
' at the 5% significance level. ',''}))    
elseif p_MLR_timedepend_tot<5/100
disp(char({' The value of the MLRC statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_timedepend_tot" and the probability ',...
' "p_MLR_timedepend_tot" for such an event is lower than 5%, ',...
' so the 1st null hypothesis is rejected by the MLR Criterion ',...
' at the 5% significance level. ',''}))
end
end
disp('   ( Press [Enter] to continue... ) ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to continue...)')
response=input(' ','s');
end
%
% Chi-Square Test and Likelihood Ratio Criterion to test the null hypothesis
% that the Markov chain from the data is a stationary first-order chain 
% respect to the alternative that it is a stationary second-order one 
% (Anderson & Goodman 1956, Sect. 3.3).
%
disp(char({'',' We want also to verify the 2nd null hypothesis ',...
' that the Markov chain from the data is a stationary 1st-order chain, ',... 
' against the alternative that it is a stationary 2nd-order one, ',...
' by the relative Chi-Square Test and the Maximum Likelihood Ratio Criterion (MLRC) ',...
' as proposed in Section 3.3 of the 1956 paper by Anderson & Goodman. ',''}))
num_k_h_l_tstart=zeros(classes_num,classes_num,classes_num,dim_incumb(2)-2);
%
%   num_k_h_l_tstart is the number of individuals making a transition
%                    from the state 'k' occupied at time t_start
%                    to the intermediate state 'h' occupied at time t_start+1
%                    and finally to the state 'l' occupied at time t_start+2,
%                    varying t_start between 1 and dim_incumb(2)-1.
%                                       
for t_start=1:dim_incumb(2)-2
for k=1:classes_num
for h=1:classes_num
for l=1:classes_num
num_k_h_l_tstart(k,h,l,t_start)=sum((incumbData(:,t_start)>=partition(k) & incumbData(:,t_start)<partition(k+1))&...
    (incumbData(:,t_start+1)>=partition(h) & incumbData(:,t_start+1)<partition(h+1))...
    &(incumbData(:,t_start+2)>=partition(l) & incumbData(:,t_start+2)<partition(l+1)));
end
end
end
end
%
%   p_stat_k_h_l    is the maximum likelihood estimation, for a Markovian 2nd order
%                   stationary chain, of the probability to find an agent in the k-th class
%                   at time t_start, in the intermediate h-th one at the date t_start+1
%                   and finally in the l-th class at t_start+2, varying t_start              
%                   between t_start=1 and t_start=dim_incumb(2)-2.
%                   
p_stat_k_h_l=zeros(classes_num,classes_num,classes_num);
num_k_h_l=squeeze(sum(num_k_h_l_tstart,4));
%
%   num_k_h_l    is the number of one-period transitions from the state 'k' at the date t_start                   
%                to the intermediate class 'h' at t_start+1 and to the final one 'l'              
%                at time t_start+2, varying t_start between 1 and dim_incumb(2)-2.
%
%   nstarM2ord_k_h  is the number of agents making a transition from the k-th class
%                   to the h-th one, between t_start and t_start+1, whatever the final state
%                   occupied at t_start+2, varying t_start between 1 and dim_incumb(2)-2.
%
nstarM2ord_k_h=squeeze(sum(num_k_h_tstart(:,:,1:(dim_incumb(2)-2)),3));
% A first alternative definition for "nstarM2ord_k_h" is:
% nstarM2ord_k_h(:,:)=(n_k_h_tult(:,:,dim_incumb(2)-1)-num_k_h_tstart(:,:,dim_incumb(2)-1));
for k=1:classes_num
for h=1:classes_num
for l=1:classes_num
if nstarM2ord_k_h(k,h)==0
   continue
end
p_stat_k_h_l(k,h,l)=num_k_h_l(k,h,l)./nstarM2ord_k_h(k,h);
end
end
end
for k=1:classes_num
for h=1:classes_num
% A second alternative definition for "nstarM2ord_k_h" is:   
% nstarM2ord_k_h(k,h)=n_k_h_tult(k,h,dim_incumb(2)-1)-num_k_h_tstart(k,h,dim_incumb(2)-1);
sum2_M2ord_k_h(k,h)=n_k_h_tult(k,h,dim_incumb(2)-1)-num_k_h_tstart(k,h,1);
end
end
% An alternative definition for "sum2_M2ord_k_h" is:
% sum2_M2ord_k_h(:,:)=sum(num_k_h_tstart(:,:,2:dim_incumb(2)-1),3);
%
%   p_statM2ord_k_h   is the maximum likelihood estimation, for a Markovian 2nd order
%                     stationary chain, of the probability to find an agent in the k-th class
%                     at time t_start+1 and in the h-th class at t_start+2, 
%                     whatever the starting state at the date t_start,
%                     varying t_start between 1 and dim_incumb(2)-2.
%          
p_statM2ord_k_h=zeros(classes_num,classes_num);
for k=1:classes_num
for h=1:classes_num
if nstar_k_tult(k,dim_incumb(2)-1)-nstar_k_tult(k,1)==0
   continue 
end
p_statM2ord_k_h(k,h)=sum2_M2ord_k_h(k,h)/(nstar_k_tult(k,dim_incumb(2)-1)-nstar_k_tult(k,1));
% An alternative definition for "p_statM2ord_k_h" is:
% p_statM2ord_k_h(k,h)=sum2_M2ord_k_h(k,h)/(nstar_k_tult(k,dim_incumb(2)-1)-n_k_tstart(k,1));
end
end
r_X_1stVS2nd_k_h_l=zeros(classes_num,classes_num,classes_num);
r_MLR_1stVS2nd_k_h_l=zeros(classes_num,classes_num,classes_num);
for k=1:classes_num
for h=1:classes_num
for l=1:classes_num
if p_statM2ord_k_h(h,l)==0
   continue
end
r_X_1stVS2nd_k_h_l(k,h,l)=((p_stat_k_h_l(k,h,l)-p_statM2ord_k_h(h,l)).^2)./p_statM2ord_k_h(h,l);
end
end
end
for k=1:classes_num
for h=1:classes_num
for l=1:classes_num
if p_statM2ord_k_h(h,l)==0
r_MLR_1stVS2nd_k_h_l(k,h,l)=1;   
else
r_MLR_1stVS2nd_k_h_l(k,h,l)=(p_statM2ord_k_h(h,l)./p_stat_k_h_l(k,h,l)).^num_k_h_l(k,h,l);
end
end
end
end
% "df" is the number of degree of freedom in both the Chi-Square Test
% and the Likelihood Ratio Criterion.
c_k=zeros(1,classes_num);
d_h=zeros(1,classes_num);
for k=1:classes_num
c_k(k)=length(find(p_statM2ord_k_h(k,:)));
end
for h=1:classes_num
d_h(h)=length(find(nstarM2ord_k_h(:,h)));
end
df_1stVS2nd_k=(c_k-1).*(d_h-1);
df_1stVS2nd_tot=(c_k-1)*(d_h-1)';
%
% "chisq_1stVS2nd_oss_h" is a column vector, relative to the imported data matrix,
% collecting, for each intermediate state 'k', the observed value of the Chi-Square statistics
% representing the null hypothesis respect to the alternative.
% "chisq_1stVS2nd_oss" is the scalar of the total observed statistics 
% for the null hypothesis against the alternative, without distinguishing
% the different intermediate states.
% Both of them admit as asymptotic distributon the Chi-Square one. 
%
chisq_1stVS2nd_oss_h=squeeze(sum(nstarM2ord_k_h.*sum(r_X_1stVS2nd_k_h_l,3),1));
chisq_1stVS2nd_oss=sum(chisq_1stVS2nd_oss_h);
%
% "p_X2_1stVS2nd_h" and "p_X2_1stVS2nd_tot" are the probabilities that the value of a chi-square variable 
% may be greater than the observed chi-square statistics, for the hypothesis of 1st order 
% Markovian stationary chain, respectively in the cases of distinguished intermediate classes 
% 'h' and without any distinction.
%
for h=1:classes_num
p_X2_1stVS2nd_h(h)=chi2cdf(chisq_1stVS2nd_oss_h(h),df_1stVS2nd_k(h),'upper');
end
p_X2_1stVS2nd_tot=chi2cdf(chisq_1stVS2nd_oss,df_1stVS2nd_tot,'upper');
disp(char({' According to the Chi Square Test, the probability that, ',...
' depending on every intermediate state ''h'', the 2nd-order stationary Markov chain, ',...
' based on the imported data, reduces itself to a stationary 1st order one starting from ''h'', ',...
' is the h-th component in the vector : ','',' p_X2_1stVS2nd_h(h) = '}))
disp(p_X2_1stVS2nd_h)
disp(char({' and the probability that the whole stationary Markov chain, ',...
' based on the data, is a 1st order chain but not a 2nd order one, is: ','',...
' p_X2_1stVS2nd_tot = '}))
disp(p_X2_1stVS2nd_tot)
if chisq_1stVS2nd_oss<=df_1stVS2nd_tot
disp(char({' The value of the Chi^2 statistics computed onto the data for the tested sample ',...
' isn''t greater than number of degree of freedom "df_1stVS2nd_tot", so we can conclude ',...
' that the 2nd null hypothesis is confirmed. '}))
elseif chisq_1stVS2nd_oss>df_1stVS2nd_tot
if p_X2_1stVS2nd_tot>=5/100
disp(char({' The value of the Chi^2 statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_1stVS2nd_tot", but the probability ',...
' "p_X2_1stVS2nd_tot" for such an event isn''t lower than 5%, ',...
' so the 2nd null hypothesis is accepted at the 5% significance level. '}))    
elseif p_X2_1stVS2nd_tot<5/100
disp(char({' The value of the Chi^2 statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_1stVS2nd_tot" and the probability ',...
' "p_X2_1stVS2nd_tot" for such an event is lower than 5%, ',...
' so the 2nd null hypothesis is rejected at the 5% significance level. '}))
end
end
disp('   ( Press [Enter] to continue... ) ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to continue...)')
response=input(' ','s');
end
%
% "MLRcrit_1stVS2nd_h" is a column vector, based on the imported data, collecting, 
% for each intermediate state 'h', the opposite of the double logarithm
% of the observed value of the Maximum Likelihood Ratio 
% for the null hypothesis against its alternative.
% "MLRcrit_1stVS2nd_tot" is the scalar of the total observed statistics
% of the Maximum Likelihood Ratio Criterion (MLRC)
% for the null hypothesis against the alternative, without distinguishing
% the different intermediate classes.
% Both of them admit, as an asymptotic distributon, the Chi-Square one. 
%
lambda_1stVS2nd_oss_h=prod(squeeze(prod(r_MLR_1stVS2nd_k_h_l,3)),1);
MLRcrit_1stVS2nd_h=-2*log(lambda_1stVS2nd_oss_h);
MLRcrit_1stVS2nd_tot=sum(MLRcrit_1stVS2nd_h);
%
% The degrees of freedom for the Maximum Likelihood Ratio Criterion 
% are the same, "df_1stVS2nd_k" and "df_1stVS2nd_tot", as for the previous Chi-Square Test.
%
% "p_MLR_1stVS2nd_h" and "p_MLR_1stVS2nd_tot" are the probabilities that the value of a chi-square variable 
% may be greater than the observed MLRC statistics, for the hypothesis of a 1st order Markovian
% stationary chain, respectively by distinguishing, or not, the intermediate states 'h'.
%
for h=1:classes_num
p_MLR_1stVS2nd_h(h)=chi2cdf(MLRcrit_1stVS2nd_h(h),df_1stVS2nd_k(h),'upper');
end
p_MLR_1stVS2nd_tot=chi2cdf(MLRcrit_1stVS2nd_tot,df_1stVS2nd_tot,'upper');
disp(char({' According to the Maximum Likelihood Ratio Criterion, ',...
' the probability that, depending on every intermediate state ''h'' ',...
' the 2nd-order stationary Markov chain, based on the imported data, ',...
' reduces itself to a stationary 1st order one starting from ''h'', ',...
' is the h-th component in the vector : ','',' p_MLR_1stVS2nd_h(h) = '}))
disp(p_MLR_1stVS2nd_h)
disp(char({' and the probability that the whole stationary Markov chain, ',...
' based on the data, is a 1st order chain but not a 2nd order one, is: ','',...
' p_MLR_1stVS2nd_tot = '}))
disp(p_MLR_1stVS2nd_tot)
if MLRcrit_1stVS2nd_tot<=df_1stVS2nd_tot
disp(char({' The value of the MLRC statistics computed onto the data for the tested sample ',...
' isn''t greater than number of degree of freedom "df_1stVS2nd_tot", so we can conclude ',...
' that the 2nd null hypothesis is confirmed by the MLR Criterion. '}))
elseif MLRcrit_1stVS2nd_tot>df_1stVS2nd_tot
if p_MLR_1stVS2nd_tot>=5/100
disp(char({' The value of the MLRC statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_1stVS2nd_tot", but the probability ',...
' "p_MLR_1stVS2nd_tot" for such an event isn''t lower than 5%, ',...
' so the 2nd null hypothesis is accepted by the MLR Criterion ',...
' at the 5% significance level. ',''}))    
elseif p_MLR_1stVS2nd_tot<5/100
disp(char({' The value of the MLRC statistics computed onto the data for the tested sample ',...
' is greater than number of degree of freedom "df_1stVS2nd_tot" and the probability ',...
' "p_MLR_1stVS2nd_tot" for such an event is lower than 5%, ',...
' so the 2nd null hypothesis is rejected by the MLR Criterion ',...
' at the 5% significance level. ',''}))
end
end
disp('    Press [Enter] to continue...  ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to go on...)')
response=input(' ','s');
end
% Computation of the mean permanence times, along with their theoretical standard deviations,
% for the variable in the datasheet estimated for a 1st ord., stationary Markov chain;
% of our two suggestions to estimate the mean permanence times for a 1st ord. not-stationary chain
% with their related theoretical standard deviations. 
% Finally, computation of the empirical mean permanence times mechanically counted in the datasheet,
% together with their corrected sample standard deviations.
for k=1:classes_num
mean_permtime_M1stat_k(k)=1./(1-p_stat_k_h(k,k));
sigma_mean_permtime_M1stat_k(k)=sqrt(p_stat_k_h(k,k))./(1-p_stat_k_h(k,k));
mean_permtime_M1nonstat_k(k)=1+sum(aggr_p_nonstat_k_h_tult(k,k,:),3);
for t_ult=1:dim_incumb(2)-1
coeffc(t_ult)=2*t_ult+1;
aggr_p_nonstat_k_k(k,t_ult)=aggr_p_nonstat_k_h_tult(k,k,t_ult);
end
end
sigma_mean_permtime_M1nonstat_k=((1+aggr_p_nonstat_k_k*(coeffc)')'-(mean_permtime_M1nonstat_k.^2)).^(1/2);
disp(char({'    "mean_permtime_M1stat_k" is the vector of average times ',...
'  spent by an individual in everyone of the k-th classes ',...
'  (defined as a geometric series whose common ration is the one-period ',...
'  transition probability from ''k'' to ''k'', "p_stat_k_h(k,k)"), ',...
'  and is computed along with the vector of the associated standard deviations, ',...
'    "sigma_mean_permtime_M1stat_k", for a stationary first-order Markov chain ',...
'  (see, for instance, the explanation in Appendix B of the 1955 Prais''s paper). ',''}))
fprintf(' If the state index ''k'' varies from 1 to %d: \n',classes_num)
disp(char({'','    mean_permtime_M1stat_k(k) = '}))
fprintf('      (all these values are expressed in ''%s'') \n',time_scan)
disp(mean_permtime_M1stat_k)
disp('    sigma_mean_permtime_M1stat_k(k) = ')
fprintf('                      (in ''%s'') \n',time_scan)
disp(sigma_mean_permtime_M1stat_k)
disp('    Press [Enter] to continue...  ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to go on...)')
response=input(' ','s');
end
disp(char({'    "mean_permtime_M1nonstat_k" is instead our proposal for the average time ',...
'  spent inside the k-th class, according to a 1st order not stationary Markov chain. ',...
'  We have defined it by analogy with the sum of the first ''dim_incumb(2)'' ',...
'  terms of the geometric series for the stationary case. ','',...
'    "mean_permtime_M1nonstat_k" is the sum of 1 and p_nonstat_k_h_tstart(k,k,1) ',...
'  with the products of the first two, three, four (and so on....)  ',...
'  not-stationary transition matrices from ''k'' to ''k'' for a 1st order Markov chain. ',...
'  It is computed together with the vector of the corresponding standard deviations, ',...
'    "sigma_mean_permtime_M1nonstat_k", for a not stationary first-order Markov chain ',''}))
fprintf(' For ''k'' varying from 1 to %d: \n',classes_num)
disp(char({'','    mean_permtime_M1nonstat_k(k) = '}))
fprintf('                      (in ''%s'') \n',time_scan)
disp(mean_permtime_M1nonstat_k)
disp('    sigma_mean_permtime_M1nonstat_k(k) = ')
fprintf('                      (in ''%s'') \n',time_scan)
disp(sigma_mean_permtime_M1nonstat_k)
disp('    Press [Enter] to continue...  ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to go on...)')
response=input(' ','s');
end
disp(char({'    "mean_permtime_observed_k" is the vector of the mean numbers ',...
'  of observation periods spent by any statistical unit ',...
'  to remain in the k-th classes: they are automatically counted ',...
'  scanning each row of the incumbents inside the spreadsheet. ',...
'  The vector of the corresponding corrected sample standard deviations ',...
'  is also computed: it is "sigma_mean_permtime_observed_k". ',...
'  Such third approach is not based on a predetermined underlying model: ',...
'  the mean permanence times and the associated sample standard deviations ',...
'  are counted row by row inside the incumbents'' matrix, ',...
'  from any date to the successive one. ',''}))
%   for k=1:classes_num
%   tot_permtime_obs_k(k)=sum(num_k_h_tstart(k,k,:),3);
%   end
%   mean_permtime_obs_k=tot_permtime_obs_k./dim_incumb(1);
%    %  "mean_permtime_obs_k" and "mean_permtime_observed_k" must coincide
%    %  if they are the same quantity computed in two different manners.
%
Agent_k_h_tstart=zeros(dim_incumb(1),classes_num,classes_num,dim_incumb(2)-1);
for Ag_id=1:dim_incumb(1)
for t_start=1:dim_incumb(2)-1
for k=1:classes_num
for h=1:classes_num
Agent_k_h_tstart(Ag_id,k,h,t_start)=(incumbData(Ag_id,t_start)>=partition(k) & incumbData(Ag_id,t_start)<partition(k+1))&...
    (incumbData(Ag_id,t_start+1)>=partition(h) & incumbData(Ag_id,t_start+1)<partition(h+1));
    %
    % "Agent_k_h_tstart(Ag_id,k,h,t_start)" is an array whose elements are logical scalar outputs:
    % if the individual identified by the integer 'Ag_id' is in state 'k' at 't_start'
    % and it will be in state 'h' at 't_start+1', such array "Agent_k_h_tstart" is worth '1',
    % otherwise it is worth zero.
end
end
end
end
for Ag_id=1:dim_incumb(1)
for k=1:classes_num
    Agent_permtime_observed_k(Ag_id,k)=sum(Agent_k_h_tstart(Ag_id,k,k,:));
    % "Agent_permtime_observed_k(Ag_id,k)" is the total number of observation periods  
    % measured for the individual, identified by the integer 'Ag_id', during which it remains 
    % in state 'k' from the time 't_start' to 't_start+1', whatever 't_start'. 
    mean_permtime_observed_k(k)=sum(Agent_permtime_observed_k(:,k))./dim_incumb(1);
    sigma_mean_permtime_observed_k(k)=((sum((Agent_permtime_observed_k(:,k)-mean_permtime_observed_k(k)).^2))...
        /(dim_incumb(1)-1)).^(1/2);
end
end
fprintf(' For ''k'' varying from 1 to %d: \n',classes_num)
disp(char({'','    mean_permtime_observed_k(k) = '}))
fprintf('                      (in ''%s'') \n',time_scan)
disp(mean_permtime_observed_k)
disp('    sigma_mean_permtime_observed_k(k) = ')
fprintf('                      (in ''%s'') \n',time_scan)
disp(sigma_mean_permtime_observed_k)
disp('    Press [Enter] to continue...  ')
response=input(' ','s');
while isempty(response)==0
disp('(press [Enter] to go on...)')
response=input(' ','s');
end
%
% Now the setting of the weights {omega_i} of the starting states
% and the function v(|i-j|) of the jumps begins, both included in the definitions
% of the family of mobility indices by Alcalde-Unzu and of the Directional Index.
disp(char({' Now it begins the setting of the weights {omega_i} of the starting states ',...
' and the function v(|i-j|) of the jumps, both included in the definitions ',... 
' of the family of mobility indices by Alcalde-Unzu and of the Directional Index. ',''}))
default=char({' By default the weights {omega_i}, are the rates of individuals ',...
' in each state "i" at the first period in the datasheet, t_start=1. ',...
' If you want to change, it is possible to assign a number, ',...
' different from the starting probability, ',...
' weighting every single initial state "i" inside such Indices; ',...
'',' Do you want to modify the default weights? ',...
' ( "Y" for ''Yes'' or "N" for ''Not'', then [Enter]) '});
Input_5ter=char({'','','        Write, in the command line, the handle  " @(x) " ',...
'        followed by the mathematical shape of your v(|i-j|), ',...
' for example   @(x)exp(x^(1/3))   or   @(x)log(2*x^3)  . ',...
'        You can express it respect to the auxiliary variable "x", ',...
'        so that x = |i-j| : anyway, to avoid conflicts, ',...
'        you must refer the handle "@" and the function to the same variable. ',...
' In any case, the Alcalde-Unzu Index and the Directional one will be computed ',...
' simultaneously for the same function v(|i-j|). ',''});
disp(default)
resp_w=input('  ','s');
while strcmp(resp_w,{'Y','y','N','n'})==[0,0,0,0]
disp('   (You must answer exclusively "Y" or "N" !) ')
resp_w=input('  ','s');
end
if sum(strcmp(resp_w,{'N','n'}))==1
for t_start=1:dim_incumb(2)-1
immob_num(t_start)=trace(num_k_h_tstart(:,:,t_start));   
for k=1:classes_num
mov_k_tstart(k,t_start)=(n_k_tstart(k,t_start)-num_k_h_tstart(k,k,t_start))./...
                      (dim_incumb(1)-immob_num(t_start));
end
end
w_k=(sum(mov_k_tstart,2)./(dim_incumb(2)-1))';
elseif sum(strcmp(resp_w,{'Y','y'}))==1
change=char({'   If you choose a scalar, as a weight for each state, ',...
'   BEWARE that THE SUM OVER ALL THE STATES MUST BE EQUAL TO 1. ',...
'','  If you dislike this option about {omega_i} ',...
'  you have to modify directly the part of this script ',...
'  concerning the definitions of  "w_k(k)". ',...
'','BE also CAREFUL NOT TO PINCH the ENTER button if you don''t want to quit !!'});
disp(change)
for k=1:classes_num
fprintf('Insert the numerical value of the %-1.0d° weight w(%-1.0d):\n',k,k)
w_k(k)=input(' ');
while w_k(k)<0 || w_k(k)>1 || isnan(w_k(k))==1 || isreal(w_k(k))==0 || isscalar(w_k(k))==0
if w_k(k)<0
    disp(' This input is a probability: you have to insert a positive number !! ')
    w_k(k)=input(' ');
elseif w_k(k)>1
    disp(char({' This input is a probability of a single event: ',...
        ' it cannot be a number greater than 1 !! '}))
    w_k(k)=input(' ');
elseif (isnan(w_k(k))==1 || isscalar(w_k(k))==0)
    disp(' This input is a probability: you have to insert a scalar !! ')
    w_k(k)=input(' ');
elseif isreal(w_k(k))==0
    disp(' This input is a probability: you have to insert a real scalar !! ')
    w_k(k)=input(' ');
end
end
end
sum_w_k=single(sum(w_k));
if (sum_w_k>1 || sum_w_k<1)
while (sum_w_k>1 || sum_w_k<1)
warn=char({' Pay attention: ',...
    ' THE SUM OF ALL THE PREVIOUS INPUTS IS NOT EQUAL TO 1 !! ',...
    ' YOU NEED TO INSERT ALL THE WEIGHTS AGAIN. ','',...
    ' Be careful not to pinch the enter button !'});
disp(warn)
for k=1:classes_num
    fprintf('Insert the numerical value of the %-1.0d° weight w(%-1.0d):\n',k,k)
    w_k(k)=input(' ');
    while w_k(k)<0 || w_k(k)>1 || isnan(w_k(k))==1 || isreal(w_k(k))==0 || isscalar(w_k(k))==0
        if w_k(k)<0
            disp(' This input is a probability: you have to insert a positive number !! ')
            w_k(k)=input(' ');
        elseif w_k(k)>1
            disp(char({' This input is a probability of a single event: ',...
                ' it cannot be a number greater than 1 !! '}))
            w_k(k)=input(' ');
        elseif (isnan(w_k(k))==1 || isscalar(w_k(k))==0)
            disp(' This input is a probability: you have to insert a scalar !! ')
            w_k(k)=input(' ');
        elseif isreal(w_k(k))==0
            disp(' This input is a probability: you have to insert a real scalar !! ')
            w_k(k)=input(' ');
        end
    end
end
sum_w_k=single(sum(w_k));
end
else
end
end
disp(Input_5ter)
clear resp_w change
%
%  We are going to introduce
%
%  V_k_h(h,k)  the suitable function of the difference (h-k), weighting
%              the significance of a jump from the k-th class to the h-th one
%                       V_k_h(h,k) = v(|h-k|)
%
%  Below an auxiliary variable, named “x”, will be used.
%
v=input('  v(|i-j|) = v(x)=  ');
%
%  The function handle “ v ”, inserted by you, corresponds to the weight function  v(|i-j|)
%  namely the handle in x=0 is the weight  v(|i-j|)  computed in |i-j|=0,
%  the handle in x=1 is the weight  v(|i-j|)  computed in |i-j|=1, and so on.
%
while isa(v,'function_handle')==0
disp(' Your input is not a function handle: ')
disp(' please follow the instructon above and retype your function. ')
v=input('  v(|i-j|) = v(x)=  ');
end
for h=1:classes_num
for k=1:classes_num
V_k_h(h,k)=feval(v,abs(h-k));
Vs_k_h(h,k)=sign(h-k)*feval(v,abs(h-k));
%
%  for the matrix V_k_h, 'h' is a row index and 'k' is a column index
%  and they appear, as parameters, also in the scalar ‘abs(h-k)’
%  which is a non negative integer.
%
end
end
V_nan=isnan(V_k_h);
V_infinite=isinf(V_k_h);
V_real=isreal(V_k_h);
while (isequal(V_nan,zeros(classes_num,classes_num))==0 | isequal(V_infinite,zeros(classes_num,classes_num))==0 |...
V_real==0)
disp(char({'','  Note that the shape  " v(|i-j|) "  you have chosen ',...
'  is not defined over the whole domain of the couple (i,j): ',...
'  it could invalidate the computation of the Index. ',...
'  You can check it by observing the position of all the possible ',...
'  NaN values, infinite values or non-real values ',...
'  in the following matrix  " V_k_h(h,k) " : ',''}))
disp(V_k_h)
disp(char({'','  So insert a more convenient functional " v(|i-j|) ", ',...
'  defined for each possible value of the couple (i,j). ',...
'  The name ''x'' for the auxiliary variable is only a suggestion: ',''}))
v=input('  v(|i-j|) = v(x)=  ');
while isa(v,'function_handle')==0
disp(' Your input is not a function handle: ')
disp(' please follow the instructon above and retype your function. ')
v=input('  v(|i-j|) = v(x)=  ');
end
for h=1:classes_num
for k=1:classes_num
V_k_h(h,k)=feval(v,abs(h-k));
end
end
V_nan=isnan(V_k_h);
V_infinite=isinf(V_k_h);
V_real=isreal(V_k_h);
end
InputAlfa=char({'','',' Lastly, it is necessary to fix the value of the power "alfa" ',...
' appearing, in the family of Mobility Indices by Alcalde-Unzu, ',...
' over the difference (P-Identity): ',...
' it is a real number, not smaller than 1, ',...
' allowing for different distance functions. '});
disp(InputAlfa)
alfa=input('    Power order "alfa" =  ');
while (isscalar(alfa)==0 || isnan(alfa)==1 || isinf(alfa)==1 || alfa<=1 || isreal(alfa)==0)
if isscalar(alfa)==0
disp(' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ')
disp('   BUT  you must insert exclusively a positive number ! ) ')
alfa=input('   ');
elseif (isscalar(alfa)==1 && (isnan(alfa)==1 || isinf(alfa)==1 || isreal(alfa)==0))
disp(' Your input is a NaN value or a not finite one or it is not a real number: WHY? ')
alfa=input('   ');
elseif alfa<=1
disp(' ( You must insert exclusively a positive number, not smaller than 1 ! ) ')
alfa=input('   ');
end
end
% Here it begins the computation of the Alcalde-Unzu Index based
% on a first-order, not-stationary Markov process
Id=eye(classes_num);
v_zero=feval(v,0);
for t_start=1:dim_incumb(2)-1
for k=1:classes_num
Max(k)=max([V_k_h(1,k),V_k_h(classes_num,k)]);
aux_Z(k)=(Max(k)+v_zero).^(1/alfa);
end
Z=aux_Z*(w_k)';
end
for t_start=1:dim_incumb(2)-1
for k=1:classes_num
Y1_nonstat_V_alfa_k_tstart(k,t_start)=(((abs(p_nonstat_k_h_tstart(k,:,t_start)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
Y1up_nonstat_V_alfa_k_tstart(k,t_start)=(((abs(p_up_nonstat_k_h_tstart(k,:,t_start)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
Y1dn_nonstat_V_alfa_k_tstart(k,t_start)=(((abs(p_dn_nonstat_k_h_tstart(k,:,t_start)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
end
Imov_w_V_alfa_tstart(t_start)=(w_k(:))'*Y1_nonstat_V_alfa_k_tstart(:,t_start);
Imovup_w_V_alfa_tstart(t_start)=(w_k(:))'*Y1up_nonstat_V_alfa_k_tstart(:,t_start);
Imovdn_w_V_alfa_tstart(t_start)=(w_k(:))'*Y1dn_nonstat_V_alfa_k_tstart(:,t_start);
end
for k=1:classes_num
Y1_stat_V_alfa_k(k)=(((abs(p_stat_k_h(k,:)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
Y1up_stat_V_alfa_k(k)=(((abs(p_up_stat_k_h(k,:)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
Y1dn_stat_V_alfa_k(k)=(((abs(p_dn_stat_k_h(k,:)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
end
Imov_stat_w_V_a=(w_k(:))'*Y1_stat_V_alfa_k(:);
Imovup_stat_w_V_a=(w_k(:))'*Y1up_stat_V_alfa_k(:);
Imovdn_stat_w_V_a=(w_k(:))'*Y1dn_stat_V_alfa_k(:);
%
%         Imov_w_V_alfa_tstart(t_start)  and  Imov_stat_w_V_a
%    are the not-normalized versions of the Mobility-as-a-Movement Index,
%    respectively for a not-stationary 1st-order Markov chain and for the stationary one.
%    The values of the Normalized Indices 
%    Imovnorm_w_V_alfa_tstart(t_start)  and  Imovnorm_stat_w_V_a
%    are:
%
Imovnorm_w_V_alfa_tstart=Imov_w_V_alfa_tstart./Z;
Imovnorm_up_w_V_alfa_tstart=Imovup_w_V_alfa_tstart./Z;
Imovnorm_dn_w_V_alfa_tstart=Imovdn_w_V_alfa_tstart./Z;
Imovnorm_stat_w_V_a=Imov_stat_w_V_a./Z;
Imovnorm_up_stat_w_V_a=Imovup_stat_w_V_a./Z;
Imovnorm_dn_stat_w_V_a=Imovdn_stat_w_V_a./Z;
StayShares_k_tstart=zeros(classes_num,dim_incumb(2)-1);
Ratemov_k_tstart=zeros(classes_num,dim_incumb(2)-1);
for k=1:classes_num
for t_start=1:dim_incumb(2)-1
StayShares_k_tstart(k,t_start)=num_k_h_tstart(k,k,t_start)./dim_incumb(1);
if (w_k(k)*Y1_nonstat_V_alfa_k_tstart(k,t_start)==0 && Imov_w_V_alfa_tstart(t_start)==0)
Ratemov_k_tstart(k,t_start)=0;
elseif (w_k(k)*Y1_nonstat_V_alfa_k_tstart(k,t_start)~=0 && Imov_w_V_alfa_tstart(t_start)==0)
Ratemov_k_tstart(k,t_start)=Inf;
else
Ratemov_k_tstart(k,t_start)=(w_k(k).*Y1_nonstat_V_alfa_k_tstart(k,t_start))./Imov_w_V_alfa_tstart(t_start);
end
end
end
Ratemov_stat_k=zeros(1,classes_num);
for k=1:classes_num
if (w_k(k)*Y1_stat_V_alfa_k(k)==0 && Imov_stat_w_V_a==0)
Ratemov_stat_k(k)=0;
elseif (w_k(k)*Y1_stat_V_alfa_k(k)~=0 && Imov_stat_w_V_a==0)
Ratemov_stat_k(k)=Inf;
else
Ratemov_stat_k(k)=(w_k(k).*Y1_stat_V_alfa_k(k))./Imov_stat_w_V_a;
end
end
% Here the computation of the Alcalde-Unzu Index ends;
% the calculus of the Directional Mobility Index, based
% on a first-order, not-stationary Markov process, begins.
clear default Input0 Input_5ter InputAlfa
m1=Vs_k_h(1,:)*w_k(:);
m2=Vs_k_h(classes_num,:)*w_k(:);
%
%  m1 and m2  are, respectively, the global minimum and the global maximum
%  theoretically possible for the not-normalized index.
%  They are necessary to carry out the normalization of the Directional Mobility Index
%  to the target [-1;+1] in the following manner:
%        if the not normalized index is negative then
%                       normalized index =  -1*(not normalized index) / m1;
%        else if the not normalized index is non negative then
%                       normalized index =  (not normalized index) / m2;
%
for t_start=1:dim_incumb(2)-1
for k=1:classes_num
Y2_nonstat_V_k_tstart(k,t_start)=p_nonstat_k_h_tstart(k,:,t_start)*Vs_k_h(:,k);
Y2up_nonstat_V_k_tstart(k,t_start)=p_up_nonstat_k_h_tstart(k,:,t_start)*Vs_k_h(:,k);
Y2dn_nonstat_V_k_tstart(k,t_start)=p_dn_nonstat_k_h_tstart(k,:,t_start)*Vs_k_h(:,k);
end
Idir_w_V_tstart(t_start)=w_k(:)'*Y2_nonstat_V_k_tstart(:,t_start);
Idirup_w_V_tstart(t_start)=w_k(:)'*Y2up_nonstat_V_k_tstart(:,t_start);
Idirdn_w_V_tstart(t_start)=w_k(:)'*Y2dn_nonstat_V_k_tstart(:,t_start);
end
%
%         Idir_w_V_tstart
%    is the not-normalized version of the Directional Mobility Index.
%         Idirup_w_V_tstart
%    and  Idirdn_w_V_tstart
%    are respectively its upward and downward components.
%    The new value of the Normalized Directional Mobility Index
%         Idirnorm_w_V_tstart
%    along with its respective upward and downward components
%    are:
%
for t_start=1:dim_incumb(2)-1
if Idir_w_V_tstart(t_start)>=0
Idirnorm_w_V_tstart(t_start)=Idir_w_V_tstart(t_start)./m2;
Idirnorm_up_w_V_tstart(t_start)=Idirup_w_V_tstart(t_start)./m2;
Idirnorm_dn_w_V_tstart(t_start)=Idirdn_w_V_tstart(t_start)./m2;
elseif Idir_w_V_tstart(t_start)<0
Idirnorm_w_V_tstart(t_start)=-Idir_w_V_tstart(t_start)./m1;
Idirnorm_up_w_V_tstart(t_start)=-Idirup_w_V_tstart(t_start)./m1;
Idirnorm_dn_w_V_tstart(t_start)=-Idirdn_w_V_tstart(t_start)./m1;
end
end
Ratio_dirup_w_V_tstart=ones(1,dim_incumb(2)-1);
for t_start=1:dim_incumb(2)-1
if (Idirnorm_up_w_V_tstart(t_start)==0 && Idirnorm_w_V_tstart(t_start)==0)
Ratio_dirup_w_V_tstart(t_start)=0;
elseif (Idirnorm_up_w_V_tstart(t_start)~=0 && Idirnorm_w_V_tstart(t_start)==0)
Ratio_dirup_w_V_tstart(t_start)=Inf;
else Ratio_dirup_w_V_tstart(t_start)=Idirnorm_up_w_V_tstart(t_start)./Idirnorm_w_V_tstart(t_start);
end
end
Ratio_dirdn_w_V_tstart=ones(1,dim_incumb(2)-1);
for t_start=1:dim_incumb(2)-1
if (Idirnorm_dn_w_V_tstart(t_start)==0 && Idirnorm_w_V_tstart(t_start)==0)
Ratio_dirdn_w_V_tstart(t_start)=0;
elseif (Idirnorm_dn_w_V_tstart(t_start)~=0 && Idirnorm_w_V_tstart(t_start)==0)
Ratio_dirdn_w_V_tstart(t_start)=Inf;
else
Ratio_dirdn_w_V_tstart(t_start)=Idirnorm_dn_w_V_tstart(t_start)/Idirnorm_w_V_tstart(t_start);
end
end
Ratiodir_k_tstart=zeros(classes_num,dim_incumb(2)-1);
for k=1:classes_num
for t_start=1:dim_incumb(2)-1
if (w_k(k)*Y2_nonstat_V_k_tstart(k,t_start)==0 && Idir_w_V_tstart(t_start)==0)
Ratiodir_k_tstart(k,t_start)=0;
elseif (w_k(k)*Y2_nonstat_V_k_tstart(k,t_start)~=0 && Idir_w_V_tstart(t_start)==0)
Ratiodir_k_tstart(k,t_start)=Inf;
else
Ratiodir_k_tstart(k,t_start)=(w_k(k).*Y2_nonstat_V_k_tstart(k,t_start))/Idir_w_V_tstart(t_start);
end
end
end
%     Computation of the Directional Mobility Index based
%     on a stationary 1st-order Markov chain:
% The not normalized version of the Index is:
for k=1:classes_num
Y2_stat_V_k(k)=p_stat_k_h(k,:)*Vs_k_h(:,k);
Y2up_stat_V_k(k)=p_up_stat_k_h(k,:)*Vs_k_h(:,k);
Y2dn_stat_V_k(k)=p_dn_stat_k_h(k,:)*Vs_k_h(:,k);
end
Idir_stat_w_V=w_k(:)'*Y2_stat_V_k(:);
Idirup_stat_w_V=w_k(:)'*Y2up_stat_V_k(:);
Idirdn_stat_w_V=w_k(:)'*Y2dn_stat_V_k(:);
% The normalized version of the Directional Index is:
if Idir_stat_w_V>=0
Idirnorm_stat_w_V=Idir_stat_w_V./m2;
Idirnorm_up_stat_w_V=Idirup_stat_w_V./m2;
Idirnorm_dn_stat_w_V=Idirdn_stat_w_V./m2;
elseif Idir_stat_w_V<0
Idirnorm_stat_w_V=-Idir_stat_w_V./m1;
Idirnorm_up_stat_w_V=-Idirup_stat_w_V./m1;
Idirnorm_dn_stat_w_V=-Idirdn_stat_w_V./m1;
end
%
% Here the calculus of the Directional Mobility Index ends.
% It begins the computation of some other mobility measures
% based on a first-order, not-stationary Markov process.
InputBet=char({'','   Also one of the indices, proposed by Shorrocks in 1978 ',...
'   and based on the determinant of the transition matrix  ',...
'   as a period invariant measure, depends on a parameter "b", ',...
'   the power over the absolute value of the determinant: ',...
'   "b" is a positive real number ',...
'   which must be fixed by the user, like the previous power order "alfa". '});
disp(InputBet)
bet=input('    Power order "b" =  ');
while (isscalar(bet)==0 || isnan(bet)==1 || isinf(bet)==1 || bet<=0 || isreal(bet)==0)
if isscalar(bet)==0
disp(' ( You pressed the [ENTER] button or inserted a vector or an array ! ) ')
disp('   BUT  you must insert exclusively a positive number ! ) ')
bet=input('   ');
elseif (isscalar(bet)==1 && (isnan(bet)==1 || isinf(bet)==1 || isreal(bet)==0))
disp(' Your input is a NaN value or a not finite one or it is not a real number: WHY? ')
bet=input('   ');
elseif bet<=0
disp(' ( You must insert exclusively a positive number ! ) ')
bet=input('   ');
end
end
for k=1:classes_num
for h=1:classes_num
Intens_k_h(k,h)=abs(h-k);
end
end
%       Here it begins the calculation of remaining seven indices 
%       based on a not-stationary, 1st-order, Markov Chain:
for t_start=1:dim_incumb(2)-1
Itrace_tstart(t_start)=(classes_num-trace(p_nonstat_k_h_tstart(:,:,t_start)))/(classes_num-1);
% Trace index of mobility
Idet_b_tstart(t_start)=1-(abs(det(p_nonstat_k_h_tstart(:,:,t_start))))^bet;
% Determinant index of mobility
% The weakness of this index when it is applied to actual raw data:
% it is sufficient that the estimated transition matrix has a null column
% or a null row to make the index be worth 1 !
Ipredict_tstart(t_start)=classes_num*(sum(sum((p_nonstat_k_h_tstart(:,:,t_start)).^2,2))-1)/(classes_num-1);
% Index of predictability
I2ndBarthol_tstart(t_start)=(classes_num/(classes_num-1))*(sum(sum(p_nonstat_k_h_tstart(:,:,t_start).*Intens_k_h)));
% Second Index of Bartholomew
eigvals_tstart(:,t_start)=eig(p_nonstat_k_h_tstart(:,:,t_start));
abs_eigvals_tstart=abs(eigvals_tstart);
sorted_abs_eigvals_tstart=sort(abs_eigvals_tstart,'descend');
I2ndeigval_tstart(t_start)=1-sorted_abs_eigvals_tstart(2,t_start); % Mobility index of the 2nd biggest eigenvalue
Ialleigvals_tstart(t_start)=(classes_num-sum(abs_eigvals_tstart(:,t_start)))/(classes_num-1);
% All eigenvalues mobility index
h_tstart(t_start)=-log(2)/log(sorted_abs_eigvals_tstart(2,t_start)); % Asymptotic half life
Ihalflife_tstart(t_start)=exp(-h_tstart(t_start)); % Asymptotic-half-life mobility Index
end
%       The alternative versions of the remaining seven indices, 
%       based on a stationary 1st-order Markov Chain, are:
Itrace_stat=(classes_num-trace(p_stat_k_h(:,:)))/(classes_num-1); % Trace index
Idet_b_stat=1-(abs(det(p_stat_k_h(:,:))))^bet; % Determinant index
Ipredict_stat=classes_num*(sum(sum((p_stat_k_h(:,:)).^2,2))-1)/(classes_num-1);
% Index of predictability
I2ndBarthol_stat=(classes_num/(classes_num-1))*(sum(sum(p_stat_k_h(:,:).*Intens_k_h)));
% Second Index of Bartholomew
eigvals_stat(:)=eig(p_stat_k_h(:,:));
abs_eigvals_stat=abs(eigvals_stat);
sorted_abs_eigvals_stat=sort(abs_eigvals_stat,'descend');
I2ndeigval_stat=1-sorted_abs_eigvals_stat(2); % Index of the 2nd biggest eigenvalue
Ialleigvals_stat=(classes_num-sum(abs_eigvals_stat(:)))/(classes_num-1);
% All eigenvalues index
h_stat=-log(2)/log(sorted_abs_eigvals_stat(2)); % Asymptotic half life
Ihalflife_stat=exp(-h_stat); % Asymptotic-half-life Index
clear InputBet
% The following section makes the calculus of the aggregate mobility indices,
% that is all the mobility indices of interest defined on a multiperiod time interval,
% both in the case of a not-stationary first-order Markov chain and for a stationary one.
%
% The Alcalde-Unzu family of aggregated mobility measures:
% when the compound (multistep) transition matrix is not-stationary
for t_ult=1:dim_incumb(2)-1
for k=1:classes_num
Y1nonstat_V_alfa_k_tult(k,t_ult)=(((abs(aggr_p_nonstat_k_h_tult(k,:,t_ult)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
end
Aggr_nonstat_Imov_w_V_alfa_tult(t_ult)=(w_k(:))'*Y1nonstat_V_alfa_k_tult(:,t_ult);
end
Aggr_nonstat_Imovnorm_w_V_alfa_tult=Aggr_nonstat_Imov_w_V_alfa_tult./Z(1);
% when the compound transition matrix is stationary   
for t_ult=1:dim_incumb(2)-1
for k=1:classes_num
Y1_stat_V_alfa_k_tult(k,t_ult)=(((abs(power_p_stat_k_h_tult(k,:,t_ult)-Id(k,:))).^alfa)*V_k_h(:,k)).^(1/alfa);
end
Aggr_stat_Imov_w_V_alfa_tult(t_ult)=(w_k(:))'*Y1_stat_V_alfa_k_tult(:,t_ult);
end
Aggr_stat_Imovnorm_w_V_alfa_tult=Aggr_stat_Imov_w_V_alfa_tult./Z(1);
%
% The aggregated Directional Mobility Index:
% when the compound (multistep) transition matrix is not-stationary
%
aggr_pup_nonstat_k_h_tult=zeros(classes_num,classes_num,dim_incumb(2)-1);
aggr_pdn_nonstat_k_h_tult=zeros(classes_num,classes_num,dim_incumb(2)-1);
for t_ult=1:dim_incumb(2)-1
aggr_pup_nonstat_k_h_tult(:,:,t_ult)=triu(aggr_p_nonstat_k_h_tult(:,:,t_ult),1);
aggr_pdn_nonstat_k_h_tult(:,:,t_ult)=tril(aggr_p_nonstat_k_h_tult(:,:,t_ult),-1);
end
for t_ult=1:dim_incumb(2)-1
for k=1:classes_num
Y2_nonstat_V_k_tult(k,t_ult)=aggr_p_nonstat_k_h_tult(k,:,t_ult)*Vs_k_h(:,k);
Y2up_nonstat_V_k_tult(k,t_ult)=aggr_pup_nonstat_k_h_tult(k,:,t_ult)*Vs_k_h(:,k);
Y2dn_nonstat_V_k_tult(k,t_ult)=aggr_pdn_nonstat_k_h_tult(k,:,t_ult)*Vs_k_h(:,k);
end
Aggr_nonstat_Idir_w_V_tult(t_ult)=w_k(:)'*Y2_nonstat_V_k_tult(:,t_ult);
Aggr_nonstat_Idirup_w_V_tult(t_ult)=w_k(:)'*Y2up_nonstat_V_k_tult(:,t_ult);
Aggr_nonstat_Idirdn_w_V_tult(t_ult)=w_k(:)'*Y2dn_nonstat_V_k_tult(:,t_ult);
end
for t_ult=1:dim_incumb(2)-1
if Aggr_nonstat_Idir_w_V_tult(t_ult)>=0
Aggr_nonstat_Idirnorm_w_V_tult(t_ult)=Aggr_nonstat_Idir_w_V_tult(t_ult)./m2;
Aggr_nonstat_Idirnorm_up_w_V_tult(t_ult)=Aggr_nonstat_Idirup_w_V_tult(t_ult)./m2;
Aggr_nonstat_Idirnorm_dn_w_V_tult(t_ult)=Aggr_nonstat_Idirdn_w_V_tult(t_ult)./m2;
elseif Aggr_nonstat_Idir_w_V_tult(t_ult)<0
Aggr_nonstat_Idirnorm_w_V_tult(t_ult)=-Aggr_nonstat_Idir_w_V_tult(t_ult)./m1;
Aggr_nonstat_Idirnorm_up_w_V_tult(t_ult)=-Aggr_nonstat_Idirup_w_V_tult(t_ult)./m1;
Aggr_nonstat_Idirnorm_dn_w_V_tult(t_ult)=-Aggr_nonstat_Idirdn_w_V_tult(t_ult)./m1;
end
end
%
% when the compound (multistep) transition matrix is stationary
aggr_pup_stat_k_h_tult=zeros(classes_num,classes_num,dim_incumb(2)-1);
aggr_pdn_stat_k_h_tult=zeros(classes_num,classes_num,dim_incumb(2)-1);
for t_ult=1:dim_incumb(2)-1
aggr_pup_stat_k_h_tult(:,:,t_ult)=triu(power_p_stat_k_h_tult(:,:,t_ult),1);
aggr_pdn_stat_k_h_tult(:,:,t_ult)=tril(power_p_stat_k_h_tult(:,:,t_ult),-1);
end
for t_ult=1:dim_incumb(2)-1
for k=1:classes_num
Y2_stat_V_k_tult(k,t_ult)=power_p_stat_k_h_tult(k,:,t_ult)*Vs_k_h(:,k);
Y2up_stat_V_k_tult(k,t_ult)=aggr_pup_stat_k_h_tult(k,:,t_ult)*Vs_k_h(:,k);
Y2dn_stat_V_k_tult(k,t_ult)=aggr_pdn_stat_k_h_tult(k,:,t_ult)*Vs_k_h(:,k);
end
Aggr_stat_Idir_w_V_tult(t_ult)=w_k(:)'*Y2_stat_V_k_tult(:,t_ult);
Aggr_stat_Idirup_w_V_tult(t_ult)=w_k(:)'*Y2up_stat_V_k_tult(:,t_ult);
Aggr_stat_Idirdn_w_V_tult(t_ult)=w_k(:)'*Y2dn_stat_V_k_tult(:,t_ult);
end
for t_ult=1:dim_incumb(2)-1
if Aggr_stat_Idir_w_V_tult(t_ult)>=0
Aggr_stat_Idirnorm_w_V_tult(t_ult)=Aggr_stat_Idir_w_V_tult(t_ult)./m2;
Aggr_stat_Idirnorm_up_w_V_tult(t_ult)=Aggr_stat_Idirup_w_V_tult(t_ult)./m2;
Aggr_stat_Idirnorm_dn_w_V_tult(t_ult)=Aggr_stat_Idirdn_w_V_tult(t_ult)./m2;
elseif Aggr_stat_Idir_w_V_tult(t_ult)<0
Aggr_stat_Idirnorm_w_V_tult(t_ult)=-Aggr_stat_Idir_w_V_tult(t_ult)./m1;
Aggr_stat_Idirnorm_up_w_V_tult(t_ult)=-Aggr_stat_Idirup_w_V_tult(t_ult)./m1;
Aggr_stat_Idirnorm_dn_w_V_tult(t_ult)=-Aggr_stat_Idirdn_w_V_tult(t_ult)./m1;
end
end
% The other mobility measures, based on the transition matrix
% of a first order Markov chain:
% when the multistep transition matrix is not-stationary
for t_ult=1:dim_incumb(2)-1
Aggr_nonstat_Itrace_tult(t_ult)=(classes_num-trace(aggr_p_nonstat_k_h_tult(:,:,t_ult)))/(classes_num-1);
% Aggregated Trace index
Aggr_nonstat_Idet_b_tult(t_ult)=1-(abs(det(aggr_p_nonstat_k_h_tult(:,:,t_ult))))^(bet/t_ult);
% Aggregated Determinant Index
Aggr_nonstat_Ipredict_tult(t_ult)=classes_num*(sum(sum((aggr_p_nonstat_k_h_tult(:,:,t_ult)).^2,2))-1)/(classes_num-1);
% Aggregated Index of predictability
Aggr_nonstat_I2ndBarthol_tult(t_ult)=(classes_num/(classes_num-1))*(sum(sum(aggr_p_nonstat_k_h_tult(:,:,t_ult).*Intens_k_h)));
% Aggregated Second Index of Bartholomew
eigvals_aggrnonstat_tult(:,t_ult)=eig(aggr_p_nonstat_k_h_tult(:,:,t_ult));
abs_eigvals_aggrnonstat_tult=abs(eigvals_aggrnonstat_tult);
sorted_abs_eigvals_aggrnonstat_tult=sort(abs_eigvals_aggrnonstat_tult,'descend');
Aggr_nonstat_I2ndeigval_tult(t_ult)=1-sorted_abs_eigvals_aggrnonstat_tult(2,t_ult); % Aggregated index of the 2nd biggest eigenvalue
Aggr_nonstat_Ialleigvals_tult(t_ult)=(classes_num-sum(abs_eigvals_aggrnonstat_tult(:,t_ult)))/(classes_num-1);
% Aggregated All eigenvalues mobility index
h_aggrnonstat_tult(t_ult)=-log(2)/log(sorted_abs_eigvals_aggrnonstat_tult(2,t_ult)); % Asymptotic half life
Aggr_nonstat_Ihalflife_tult(t_ult)=exp(-h_aggrnonstat_tult(t_ult)*t_ult); % Aggregated Asymptotic-half-life mobility Index
end
% when the multistep transition matrix is stationary
for t_ult=1:dim_incumb(2)-1
Aggr_stat_Itrace_tult(t_ult)=(classes_num-trace(power_p_stat_k_h_tult(:,:,t_ult)))/(classes_num-1);
% Aggregated Trace index
Aggr_stat_Idet_b_tult(t_ult)=1-(abs(det(power_p_stat_k_h_tult(:,:,t_ult))))^(bet/t_ult);
% Aggregated Determinant Index
Aggr_stat_Ipredict_tult(t_ult)=classes_num*(sum(sum((power_p_stat_k_h_tult(:,:,t_ult)).^2,2))-1)/(classes_num-1);
% Aggregated Index of predictability
Aggr_stat_I2ndBarthol_tult(t_ult)=(classes_num/(classes_num-1))*(sum(sum(power_p_stat_k_h_tult(:,:,t_ult).*Intens_k_h)));
% Aggregated Second Index of Bartholomew
eigvals_aggrstat_tult(:,t_ult)=eig(power_p_stat_k_h_tult(:,:,t_ult));
abs_eigvals_aggrstat_tult=abs(eigvals_aggrstat_tult);
sorted_abs_eigvals_aggrstat_tult=sort(abs_eigvals_aggrstat_tult,'descend');
Aggr_stat_I2ndeigval_tult(t_ult)=1-sorted_abs_eigvals_aggrstat_tult(2,t_ult); % Aggregated index of the 2nd biggest eigenvalue
Aggr_stat_Ialleigvals_tult(t_ult)=(classes_num-sum(abs_eigvals_aggrstat_tult(:,t_ult)))/(classes_num-1);
% Aggregated All eigenvalues mobility index
h_aggrstat_tult(t_ult)=-log(2)/log(sorted_abs_eigvals_aggrstat_tult(2,t_ult)); % Asymptotic half life
Aggr_stat_Ihalflife_tult(t_ult)=exp(-h_aggrstat_tult(t_ult)*t_ult); % Aggregated Asymptotic-half-life mobility Index
end
%
% In the following lines are defined the bidimensional arrays
% where the values of all the mobilities indices are gathered
% (one particular index for each row) for every observation period
% (one column for each date).
nonstat_Mobilities=[Idirnorm_w_V_tstart;Imovnorm_w_V_alfa_tstart;Itrace_tstart;Idet_b_tstart;...
    Ipredict_tstart;I2ndBarthol_tstart;I2ndeigval_tstart;Ialleigvals_tstart;Ihalflife_tstart];
stat_Mobilities=[Aggr_stat_Idirnorm_w_V_tult(1);Aggr_stat_Imovnorm_w_V_alfa_tult(1);...
    Aggr_stat_Itrace_tult(1);Aggr_stat_Idet_b_tult(1);Aggr_stat_Ipredict_tult(1);...
    Aggr_stat_I2ndBarthol_tult(1);Aggr_stat_I2ndeigval_tult(1);Aggr_stat_Ialleigvals_tult(1);...
    Aggr_stat_Ihalflife_tult(1)];
Aggr_nonstat_Mobilities=[Aggr_nonstat_Idirnorm_w_V_tult;Aggr_nonstat_Imovnorm_w_V_alfa_tult;...
    Aggr_nonstat_Itrace_tult;Aggr_nonstat_Idet_b_tult;...
    Aggr_nonstat_Ipredict_tult;Aggr_nonstat_I2ndBarthol_tult;...
    Aggr_nonstat_I2ndeigval_tult;Aggr_nonstat_Ialleigvals_tult;Aggr_nonstat_Ihalflife_tult];
Aggr_stat_Mobilities=[Aggr_stat_Idirnorm_w_V_tult;Aggr_stat_Imovnorm_w_V_alfa_tult;...
    Aggr_stat_Itrace_tult;Aggr_stat_Idet_b_tult;Aggr_stat_Ipredict_tult;...
    Aggr_stat_I2ndBarthol_tult;Aggr_stat_I2ndeigval_tult;Aggr_stat_Ialleigvals_tult;...
    Aggr_stat_Ihalflife_tult];
clear Id coeffc aggr_p_nonstat_k_k
%
% GRAPHICS
for t_start=1:dim_incumb(2)-1
t(t_start)=t_start+1;
end
% One-period mobility indices vs the final time
% of any possible one-step, not-stationary, 1st ord. Markovian transition 
% (from t_start to t_start+1) from the datasheet.
nonstat_dir_mov_trace=figure('Name',...
'1st ord., not-stationary Markov chain: Directional Mobility, Mobility as a movement, Trace Index',...
'NumberTitle','off');
scatter(t,Idirnorm_w_V_tstart,40,'d','filled')
hold on
scatter(t,Imovnorm_w_V_alfa_tstart,40,'o','filled')
scatter(t,Itrace_tstart,50,'g','s','filled')
line(t,Idirnorm_w_V_tstart)
line(t,Imovnorm_w_V_alfa_tstart)
line(t,Itrace_tstart)
hold off
legend('Directional Mobility','Mobility as a movement','Trace Mobility Index')
nonstat_halflife_predic_Barthol=figure('Name',...
'1st ord., not-stationary Markov chain: Asymptotic-half-life Mobility, Index of predictability, 2nd Bartholomew Ind.',...
'NumberTitle','off');
scatter(t,Ihalflife_tstart,50,'g','s','filled')
hold on
scatter(t,Ipredict_tstart,40,'o','filled')
scatter(t,I2ndBarthol_tstart,40,'b','d','filled')
line(t,Ihalflife_tstart)
line(t,Ipredict_tstart)
line(t,I2ndBarthol_tstart)
hold off
legend('Asymptotic-half-life Ind.','Ind. of predictability','2nd Bartholomew Ind.')
nonstat_2ndEig_allEig_det=figure('Name',...
'1st ord., not-stationary Markov chain: 2nd Eigenvalue Mobility Index, All-eigenvalue Index, Determinant Index',...
'NumberTitle','off');
scatter(t,I2ndeigval_tstart,40,'d','filled')
hold on
scatter(t,Ialleigvals_tstart,40,'o','filled')
scatter(t,Idet_b_tstart,50,'g','s','filled')
line(t,I2ndeigval_tstart)
line(t,Ialleigvals_tstart)
line(t,Idet_b_tstart)
hold off
legend('2nd-Eigenvalue Mobil. Ind.','All-eigenvalue Ind.','Determinant Mobility Ind.')
%
% Aggregated-in-time mobility indices vs the final time t_ult
% (from the first period to each of the successive ones) 
% of any possible multistep, not-stationary transition 
% based on a 1st ord. Markov chain from the datasheet.
aggnonstat_dir_mov_trace=figure('Name',...
'1st ord., not-stationary Markov ch.: aggregation in time-Directional Mobility, Mob. as a movement, Trace Ind.',...
'NumberTitle','off');
scatter(t,Aggr_nonstat_Idirnorm_w_V_tult,40,'d','filled')
hold on
scatter(t,Aggr_nonstat_Imovnorm_w_V_alfa_tult,40,'o','filled')
scatter(t,Aggr_nonstat_Itrace_tult,50,'g','s','filled')
line(t,Aggr_nonstat_Idirnorm_w_V_tult)
line(t,Aggr_nonstat_Imovnorm_w_V_alfa_tult)
line(t,Aggr_nonstat_Itrace_tult)
hold off
legend('Aggregated Directional Ind.','Aggregat. Mobility as a movement','Aggregated Trace Ind.')
aggnonstat_halflife_predic_Barthol=figure('Name',...
'1st ord., not-stationary Markov ch.: aggregation in time-Asymptotic-half-life Ind., Ind. of predictability, 2nd Bartholomew Ind.',...
'NumberTitle','off');
scatter(t,Aggr_nonstat_Ihalflife_tult,50,'g','s','filled')
hold on
scatter(t,Aggr_nonstat_Ipredict_tult,40,'o','filled')
scatter(t,Aggr_nonstat_I2ndBarthol_tult,40,'d','filled')
line(t,Aggr_nonstat_Ihalflife_tult)
line(t,Aggr_nonstat_Ipredict_tult)
line(t,Aggr_nonstat_I2ndBarthol_tult)
hold off
legend('Aggregated Asymptotic-half-life Ind.','Aggreg. Ind. of predictability','Aggreg. 2nd Bartholomew Ind.')
aggnonstat_2ndEig_allEig_det=figure('Name',...
'1st ord., not-stationary Markov ch.: aggregation in time-2nd Eigenvalue Mobil. Ind., All-eigenvalue Ind., Determinant Ind.',...
'NumberTitle','off');
scatter(t,Aggr_nonstat_I2ndeigval_tult,40,'b','d','filled')
hold on
scatter(t,Aggr_nonstat_Ialleigvals_tult,40,'o','filled')
scatter(t,Aggr_nonstat_Idet_b_tult,40,'g','s','filled')
line(t,Aggr_nonstat_I2ndeigval_tult)
line(t,Aggr_nonstat_Ialleigvals_tult)
line(t,Aggr_nonstat_Idet_b_tult)
hold off
legend('Aggregated 2nd-Eigenvalue Ind.','Aggreg. All-eigenvalue Ind.','Aggreg. Determinant Ind.')
%
% Aggregated-in-time mobility indices vs the final time t_ult
% (from the first period to each of the successive ones) 
% of any possible multistep, stationary transition 
% based on a 1st ord. Markov chain from the datasheet.
aggstat_dir_mov_trace=figure('Name',...
'1st ord., stationary Markov ch.: aggregation in time-Directional Mobility, Mob. as a movement, Trace Ind.',...
'NumberTitle','off');
scatter(t,Aggr_stat_Idirnorm_w_V_tult,40,'d','filled')
hold on
scatter(t,Aggr_stat_Imovnorm_w_V_alfa_tult,40,'o','filled')
scatter(t,Aggr_stat_Itrace_tult,50,'g','s','filled')
line(t,Aggr_stat_Idirnorm_w_V_tult)
line(t,Aggr_stat_Imovnorm_w_V_alfa_tult)
line(t,Aggr_stat_Itrace_tult)
hold off
legend('Aggregated Directional Ind.','Aggregat. Mobility as a movement','Aggregated Trace Ind.')
aggstat_halflife_predic_Barthol=figure('Name',...
'1st ord., stationary Markov ch.: aggregation in time-Asymptotic-half-life Ind., Ind. of predictability, 2nd Bartholomew Ind.',...
'NumberTitle','off');
scatter(t,Aggr_stat_Ihalflife_tult,50,'g','s','filled')
hold on
scatter(t,Aggr_stat_Ipredict_tult,40,'o','filled')
scatter(t,Aggr_stat_I2ndBarthol_tult,40,'d','filled')
line(t,Aggr_stat_Ihalflife_tult)
line(t,Aggr_stat_Ipredict_tult)
line(t,Aggr_stat_I2ndBarthol_tult)
hold off
legend('Aggregated Asymptotic-half-life Ind.','Aggreg. Ind. of predictability','Aggreg. 2nd Bartholomew Ind.')
aggstat_2ndEig_allEig_det=figure('Name',...
'1st ord., stationary Markov ch.: aggregation in time-2nd Eigenvalue Mobil. Ind., All-eigenvalue Ind., Determinant Ind.',...
'NumberTitle','off');
scatter(t,Aggr_stat_I2ndeigval_tult,40,'b','d','filled')
hold on
scatter(t,Aggr_stat_Ialleigvals_tult,40,'o','filled')
scatter(t,Aggr_stat_Idet_b_tult,40,'g','s','filled')
line(t,Aggr_stat_I2ndeigval_tult)
line(t,Aggr_stat_Ialleigvals_tult)
line(t,Aggr_stat_Idet_b_tult)
hold off
legend('Aggregated 2nd-Eigenvalue Ind.','Aggreg. All-eigenvalue Ind.','Aggreg. Determinant Ind.')
% End of graphical section
clear response k h l t_start t_ult Ag_id
