% Simulation Analysis Utility (Created 09/23/18; Last Modified 03/22/19)

%% Introduction

% NOTE 061619: Implemented a more efficient approach to the
% microsimulation process, which should, in general, reduce processing
% times by 90% or more.
% 
% NOTE 060604: The script now initializes the random number generator using
% the same seed, which should promote reproducibility by ensuring the same
% numbers are drawn each time - and, consequently, that we should always
% see the same results (provided, of course, nothing else changes from one
% simulation to the next).
% 
% NOTE 032219: Incorporated individual rewards. The intent of this is to
% calculate the present value of expected lifetime well-being for the
% cohort and for population sub-groups.
% 
% NOTE 121618: Redefined premature mortality in terms of death prior to
% expected retirement among individuals already diagnosed with T2DM -
% rather than ALL individuals, as had previously been the case.
% 
% This script implements a microsimulation analysis that compares economic
% losses attributable to absenteeism, presenteeism and premature mortality
% between two scenarios, namely 
% a) the optimal policy resulting from the application of dynamic 
% programming, and; 
% b) an alternate user-specified policy in which the decision-maker 
%       i) does not adhere to ANY health-promoting behaviour;
%       ii) uniformly adheres to ALL health-promoting behaviours;
%       iii) does not adhere to any health-promoting behaviours while
%       healthy, but adopts the optimal policy thereafter.
% 
% * It should be noted that the alternate policy is specified in the same
% script as the program that solves the DP problem; in the future, it may
% be advisable to migrate that segment of the script here, to facilitate
% ease of use.

%% Initialization

% 090918: Stop execution of code by typing CTRL + C.
% 100418: For information regarding (1) conversion of data from one type to
% another or (2) data formatting, see: 
% 1) https://www.mathworks.com/help/matlab/data-type-conversion.html
% 2) https://www.mathworks.com/help/matlab/ref/sprintf.html
clearvars; clc;

% BASE refers to the directory in which the model is stored. The following
% line returns the full path to the current script, but then excises
% everything after "Diabetes_Model_090618". The intent is to simplify model
% implementation for the user by avoiding the need to adjust file pathways
% manually.
% tic
BASE = regexprep(matlab.desktop.editor.getActiveFilename, ...
    'Diabetes_Model_042019.*$', 'Diabetes_Model_042019');
% NOTE 061619: The next line acts as a switch determining whether MATLAB's
% random number generator should use the nonnegative integer seed as the
% author does in his thesis. If equal to 1, MATLAB will initialize the
% microsimulation using the integer seed '1234'; otherwise, it will run
% using whatever the current seed may be.
Thesis_Results = 1;

Ref_Path = strcat(BASE, "\R_Objects\");
MAT_Path = strcat(BASE, "\MATLAB\");

% R_alpha, R_beta and R_gamma indicate the integer ranges over which all
% loops in this script (of which there are several) should run. Specifying
% the ranges here prevents us from having to modify these manually if, for
% example, we wish to focus on a specific subset of the cohort.
% 
% 091518: "Alpha" denotes gender (i.e., male and female); "beta" denotes 
% ethnicity (i.e., Afro-Caribbean, Asian-Indian, and Caucasian); and 
% "gamma" denotes education level (i.e., no high school degree, high school
% degree or equivalent, certificate or diploma below bachelor level, 
% certificate, diploma or degree at bachelor level or above).
R_alpha = 1:2; R_beta = 1:3; R_gamma = 1:4;

%% Designing and Implementing the Simulation

% Now we carry out a series of simluations, using the MDPSolve package: 
% Fackler, P. L. (2014). MDPSOLVE: MATLAB tools for solving Markov 
% Decision Problems (Version 2014a). Raleigh, NC. Retrieved from 
% https://sites.google.com/site/mdpsolve/

% 100918: Load starting cohort used in the simulations
Starting_Cohort = readtable(strcat(MAT_Path, ...
    "Simulation_Results\", "Starting_Cohort", ".dat"), ...
    'TextType', 'String');

% Create reference tables to link integer values (such as we will be using
% in our for-loops) to individual characteristics.
Gend_Ref = array2table([transpose([2 1]) unique(Starting_Cohort.Gender)]);
Ethn_Ref = array2table([transpose([2 3 1]) ...
    unique(Starting_Cohort.Ethnicity)]);
Educ_Ref = array2table([transpose(1:4) unique(Starting_Cohort.Education)]);

% Var1 in each table is the integer value, but it is currently stored as a 
% string. We need to convert it to a number before we proceed. Note that
% the function STR2NUM seems not to work as intended when applied below
% ("Input must be a character vector or string scalar").
[Gend_Ref.Var1, Ethn_Ref.Var1, Educ_Ref.Var1] = ...
    deal(str2double(Gend_Ref.Var1), str2double(Ethn_Ref.Var1), ...
    str2double(Educ_Ref.Var1));

% Generate labels for each of the tables produced above.
Gend_Labels = cellstr({"ID_Gend", "Gender"}); 
Ethn_Labels = cellstr({"ID_Ethn", "Ethnicity"});
Educ_Labels = cellstr({"ID_Educ", "Education"});

% Now apply the labels by calling (and modifying) the properties of each
% table.
Gend_Ref.Properties.VariableNames = Gend_Labels;
Ethn_Ref.Properties.VariableNames = Ethn_Labels;
Educ_Ref.Properties.VariableNames = Educ_Labels;

Ref_Tables = {Gend_Ref Ethn_Ref Educ_Ref};
% Finally, carry out a series of outer joins to merge the starting cohort
% with the reference tables.
for i = 1:3
    Starting_Cohort = outerjoin(Starting_Cohort, Ref_Tables{1,i}, ...
        'Type', 'Left', 'MergeKeys', true);
end
% Sort the resultant table by gender, then by ethnicity, and finally by
% level of educational attainment.
Starting_Cohort = sortrows(Starting_Cohort, ...
    {'ID_Gend', 'ID_Ethn', 'ID_Educ'});

% 101218: Create a series of cell arrays to hold the simulation results.
% We will effectively run a series of 24 simulations (one for each possible
% combination of individual characteristics), and then integrate the
% results.
Sim_Results_Size = cell(size(Gend_Ref, 1), size(Ethn_Ref, 1), ...
    size(Educ_Ref, 1));
% Create several new cell arrays with the same dimensions as
% Sim_Results_Size, through the application of the DEAL function.
% http://stackoverflow.com/questions/5158032/define-multiple-variables-at-the-same-time-in-matlab/5158206#5158206
[SR_Opt, SR_Opt_Cohort, RWD_Opt, SR_Alt_A, SR_Alt_A_Cohort, RWD_Alt_A, ...
    SR_Alt_DA, SR_Alt_DA_Cohort, RWD_Alt_DA, SR_Alt_HB, ...
    SR_Alt_HB_Cohort, RWD_Alt_HB, SR_ID, SR_PM_Opt, SR_PM_Alt_A, ...
    SR_PM_Alt_DA, SR_PM_Alt_HB, EL_Alive, EL_Dead_00, EL_Dead_01, ...
    LM_Part] = deal(Sim_Results_Size);

% - T denotes the number of time periods we wish to simulate. The final
% simulations run for a total of 70 periods.
T = 70;

% 0606019: In order to promote reproducibility of simulation results, the
% script now initializes the random number generator using the same seed,
% meaning we should obtain the same set of random numbers - and, therefore,
% the same results - every time anyone performs the simulation. Note that
% the choice of seeds depends on the value of the variable Thesis_Results. 
% If set equal to 1, the program uses the seed chosen by the author; if set
% to anything else, the microsimulation will run using whatever the current
% seed may be. See:
% https://www.mathworks.com/help/matlab/ref/rng.html
% https://www.mathworks.com/help/matlab/math/generate-random-numbers-that-are-repeatable.html
if Thesis_Results == 1
    rng(1234, 'twister');
else
end

tic
for alpha = R_alpha
   for beta = R_beta
       for gamma = R_gamma

Person_File = [alpha, beta, gamma];
% Alpha, beta, and gamma are all integers; before we can concatenate them,
% we must therefore convert each one from an integer to a string.
Case = strcat(num2str(Person_File(1)), num2str(Person_File(2)), ...
    num2str(Person_File(3)));
           
% Load the optimal probability transition matrices - one for the optimal
% policy, and a second for the alternative policy.
pstar_opt = load(strcat(MAT_Path, "Simulation_Results\", ...
    "pstar_opt_cum_", Case, ".mat"));
pstar_alt_A = load(strcat(MAT_Path, "Simulation_Results\", ...
    "pstar_alt_cum_", Case, "_1", ".mat"));
pstar_alt_DA = load(strcat(MAT_Path, "Simulation_Results\", ...
    "pstar_alt_cum_", Case, "_8", ".mat"));
pstar_alt_HB = load(strcat(MAT_Path, "Simulation_Results\", ...
    "pstar_alt_cum_", Case, "_hybrid", ".mat"));

% For whatever reason, the two matrices are embedded in a structure array;
% we extract the data we want (and, incidentally, the only data in the
% array) through the use of dot notation.
pstar_opt = pstar_opt.pstar_opt_cum;
pstar_alt_A = pstar_alt_A.pstar_alt_cum_temp;
pstar_alt_DA = pstar_alt_DA.pstar_alt_cum_temp;
pstar_alt_HB = pstar_alt_HB.pstar_alt_cum_temp;

% - s0 denotes the initial state of our simulation patients. The size of s0
% (i.e., the number of patients) is tied to the composition of the starting
% cohort, which is in turn based upon the number of 25 year-olds in the
% Province of Manitoba in July 2017, and the distribution of individual
% characteristics contained in the 25-29 age bracket in the 2011 NHS PUMF
% for all of Canada.
% 110418: I observed that Start_Cohort.ID_Gend was set equal to Person_File
% (2) while Starting_Cohort.ID_Ethn was set equal to Person_File(1). I
% believe this to have been an error - Person_File = [alpha, beta, gamma],
% where alpha, beta and gamma denote gender, ethnicity, and education
% level, respectively (i.e., gender and ethnicity appear to have been
% accidentially reversed).
s0 = ones( ...
    Starting_Cohort.Number(Starting_Cohort.ID_Gend == Person_File(1) & ...
    Starting_Cohort.ID_Ethn == Person_File(2)& ...
    Starting_Cohort.ID_Educ == Person_File(3)), 1);

% 101218: Here is where the simulation is actually performed. Originally,
% the intent was to accomplish this using the MDPSIM function included in
% Paul Fackler's MDP Solve function; however, this was ultimately found to
% be unworkable. This is because the simulation process works by producing
% an array, each row of which consists of probabilities accumulated over
% previous columns, which implies that every cell must be filled. However,
% the array in question is enormous (i.e., 52.1 billion cells) and would
% not fit in memory. 
% 
% Therefore, we were required to seek out an alternative approach that 
% modified the original function. This was accomplished by calculating
% cumulative transition probabilities in such a way that cell values are
% non-empty only when there is, in fact, something to add. This sidesteps
% the need to fill every cell inbetween, freeing up a substantial amount of
% memory.

% k is simply the number of units to be simulated.
k = length(s0);
% SI holds the simulation data. It has as many rows as there are patients,
% and as many columns as there are periods in the simulation - plus one to
% reflect starting conditions.
SI = zeros(k,T+1);
% ns = length(s0);
% u  = ones(ns,1);
s = s0(:);
% The first column of the simulation results 
SI_opt(:,1) = s; SI_alt_A(:,1) = s; SI_alt_DA(:,1) = s; SI_alt_HB(:,1) = s;

% 061519: This is the core of the simulation process. The current approach
% is the third I've employed, and the most efficient to date (see:
% https://www.mathworks.com/matlabcentral/answers/365710-how-to-find-first-1-in-every-row).
% It involves the following steps:
% 1) Generate a matrix that pulls the requisite number of random draws from
% the uniform distribution.
% 2) As already noted, the cumulative transition probability matrix is an
% nxn array (where n is the number of states in the decision model), where
% each row and column denote starting and ending states, respectively. Each
% entry, therefore, denotes the likelihood of transitioning from state n_x
% to n_y or any "earlier" state. The SI_opt matrix, meanwhile indicates the
% starting state for each person in the cohort. The expression 
% pstar_opt(SI_opt(:, t-1), :), therefore, pulls out a single row of the
% cumulative transition probability matrix for every individual in the
% cohort, where the row in question is the state they occupy at a point in
% time. The rows extracted in this way will consist of only a few entries,
% each of which denotes the likelihood of transitioning to that state or
% any "earlier" state (states are said to be "earlier" solely in the sense
% that they are associated with a lower index number).
% 3) Each vector of possible transitions is then compared against a single
% random draw from the matrix "r" for every patient in the cohort. If the
% cumulative transition probably is larger than the random draw, it is
% assigned a value of 1; otherwise, it is assigned a value of 0. Each row
% of the resultant matrix will therefore consist of a series of 0's and
% 1's, where the first 1 indicates the first cumulative transition
% probability larger than the random draw, and, therefore, the transition
% the individual will actually make.
% 4) Apply the MAX function, the syntax for which tells us that the result
% should be a column vector indicating the maximum value in each row. As
% just noted, each row consists of 0's and 1's, so the maximum value will
% of course be 1. The beauty of the MAX function is that it appear to
% select only the first such value, which, again, is the first cumulative
% transition probability larger than the random draw. Moreover, the LHS of
% the expression tells us we are not interested in the maximum value itself
% (we already know what that is), but rather the index number associated
% with the column in which it appears.
% 5) Finally, the expression loops over every time period included in the
% simulation. Unfortunately, it is not possible to do away with the loops
% altogether, because which transition probabilities should be compared
% against the random draws at time t requires us to know which state they
% occupied at time t-1. As such, we are required to proceed iteratively.
r = rand(k,T);
tic
for t = 2:T+1
    [~, SI_opt(:, t)] = ...
        max(pstar_opt(SI_opt(:, t-1), :) > r(:, t-1), [], 2);
    [~, SI_alt_A(:, t)] = ...
        max(pstar_alt_A(SI_alt_A(:, t-1), :) > r(:, t-1), [], 2);
    [~, SI_alt_DA(:, t)] = ...
        max(pstar_alt_DA(SI_alt_DA(:, t-1), :) > r(:, t-1), [], 2);
    [~, SI_alt_HB(:, t)] = ...
        max(pstar_alt_HB(SI_alt_HB(:, t-1), :) > r(:, t-1), [], 2);
end
toc

SI_opt = array2table(SI_opt); SI_alt_A = array2table(SI_alt_A);
SI_alt_DA = array2table(SI_alt_DA); SI_alt_HB = array2table(SI_alt_HB);
% Create a list of identifiers that indicate the patients' gender,
% ethnicity, and education. We do this so that we can continue to
% distinguish between patients after we concatenate the data associated
% with each sub-simulation. Note the use of the REPMAT function, which
% enables us to create a column containing as many of each value as the
% number of patients incorporated into the simulation.
SI_ID = array2table([repmat(alpha, k, 1), repmat(beta, k, 1), ...
    repmat(gamma, k, 1)]);
% I'm not altogether sure why the CELLSTR function is necessary.
SI_Names = cellstr({"Gender", "Ethnicity", "Education"});
% Rename the three identifier columns appropriately.
SI_ID.Properties.VariableNames = SI_Names;
% Concatenate the identifier columns with the two sets of simulation
% results.
SI_opt = [SI_ID SI_opt]; SI_alt_A = [SI_ID SI_alt_A];
SI_alt_DA = [SI_ID SI_alt_DA]; SI_alt_HB = [SI_ID SI_alt_HB];
% Assign the simulation results to the appropriate slot within the cell
% arrays.
SR_Opt{alpha, beta, gamma} = SI_opt;
SR_Alt_A{alpha, beta, gamma} = SI_alt_A;
SR_Alt_DA{alpha, beta, gamma} = SI_alt_DA;
SR_Alt_HB{alpha, beta, gamma} = SI_alt_HB;

% 101518: When running a simulation involving people with two distinct sets
% of personal characteristics, I encountered the following error: 
% "To assign to or create a variable in a table, the number of rows must 
% match the height of the table." The issue appears to be that data
% constructs continue to accumulate information across loops. Some of this
% is desirable, and indeed necessary. However, issues can emerge if we do
% not "clear the slate". For example (and to the point), it appears that 
% tables are not in fact overwritten from one loop to the next.
% Consequently, MATLAB seems to attempt to paste new content from each such
% loop into a table that already contains information from the previous
% loop. In this case, I believe MATLAB was trying to integrate data from
% two sub-populations and was unable to do so.
% 
% 101918: The notation below (i.e., "-except", followed by a list of
% variables) enables us to specify which variables to keep.

clearvars -except alpha beta gamma R_alpha R_beta R_gamma Ref_Path ...
    MAT_Path T Starting_Cohort SR_Opt SR_Alt_A SR_Alt_DA SR_Alt_HB ...
    RWD_Opt RWD_Alt_A RWD_Alt_DA RWD_Alt_HB SR_ID SR_PM_Opt SR_PM_Alt_A ...
    SR_PM_Alt_DA SR_PM_Alt_HB EL_Alive EL_Dead_00 EL_Dead_01 LM_Part rng

       end
   end
end
toc

% This integrates data from all of the sub-simulations together into a
% single dataset.
SR_Opt_Cohort = vertcat(SR_Opt{:}); SR_Alt_A_Cohort = vertcat(SR_Alt_A{:});
SR_Alt_DA_Cohort = vertcat(SR_Alt_DA{:});
SR_Alt_HB_Cohort = vertcat(SR_Alt_HB{:});
% NOTE: 121618 - Save the workspace to disk so it is not necessary to
% constantly re-run the simulations (this was not down later in the code
% because the operations ultimately consume a significant amount of disk
% space).
save(strcat(MAT_Path, "Simulation_Results\", "Sim_Results.mat"), '-v7.3');

%% Interpreting Simulation Results

% Call up the RefStates table, which links the integer state ID to a
% qualitative description that faciliates interpretation. Note that
% 1) The TreatAsEmpty argument converts NUMERIC values of specific types
% (here, "." or "NA") to NaN. Note that given the formatting applied to the
% table (i.e., only the first column is really numeric), the argument has 
% no practical implications.
% 2) The Format argument tells us how each column in the table should be
% formatted (where "%" distinguishes between columns). For example, "u"
% instructs MATLAB to format as a Base 10 integer, while "s" instructs it
% to format as a character vector or string array ("f" is for floating-
% point numbers).
% https://www.mathworks.com/help/matlab/matlab_prog/clean-messy-and-missing-data-in-tables.html
% https://www.mathworks.com/help/matlab/ref/sprintf.html
% https://www.mathworks.com/help/matlab/ref/readtable.html
RefStates = readtable(strcat(Ref_Path, "Ref_States.dat"), ... 
    'TreatAsEmpty',{'.','NA'}, 'Format', '%u%s%s%s%s%s%s%s%s%s%s');
% EffLosses_Alive was tabulated directly from 
% RWD&TP_Matrices_Construction_090618, and (as the name suggests) is a
% reference table that ties health status to efficiency losses for
% individuals who are alive.

for alpha = R_alpha
   for beta = R_beta
       for gamma = R_gamma

Case = strcat(string(alpha), string(beta), string(gamma));

% Import estimates of efficiency losses attributable to absenteeism and
% presenteeism from the generation of the reward and transition probability
% matrices.
EL_Alive{alpha, beta, gamma} = ...
    readtable(strcat(MAT_Path, "Simulation_Results\",...
    "Efficiency_Loss_Calculations\", "Eff_Losses_Alive_", Case, ...
    ".dat"), 'TreatAsEmpty',{'.','NA'}, 'Format', ...
    '%u%s%s%s%s%s%s%s%s%s%s%f%f%f');

% Before importing choices relating to labour market participation from the
% optimal policies, we first set import options. Note that we could have
% done this with the previous dataset, but simply chose not to.
% 
% The detectImportOptions function enables us to control how tabular data
% is imported into MATLAB. First we need to know what the options are.
% https://www.mathworks.com/help/matlab/ref/detectimportoptions.html
opts = detectImportOptions(strcat(MAT_Path, "Baseline_Results\", ...
    "MATLAB_Output\", "Optimal_Policies\", "Case_", Case, ...
    "\OptPolicy_", Case, ".dat"));
% Next, we change the import data type for the variables
% ID_States_Opt_Pol_01 and LMP
% https://www.mathworks.com/help/matlab/ref/matlab.io.text.delimitedtextimportoptions.setvartype.html
opts = setvartype(opts,{'ID_States_Opt_Pol_01'},'int32');
opts = setvartype(opts,{'LMP'},'char');
% We can always verify that variables will be in the format we want upon
% import by applying the GETVAROPTS function.
getvaropts(opts,{'ID_States_Opt_Pol_01','LMP'});
% We don't need or want the entire table, only two variables of interest.
opts.SelectedVariableNames = {'ID_States_Opt_Pol_01', 'LMP'};
% Having set the import parameters the way we wanted them, we now read the
% table into MATLAB.
LM_Part{alpha, beta, gamma} = readtable(strcat(MAT_Path, ...
    "Baseline_Results\", "MATLAB_Output\", "Optimal_Policies\", ...
    "Case_", Case, "\OptPolicy_", Case, ".dat"), opts);
% https://www.mathworks.com/matlabcentral/answers/179129-is-it-possible-to-rename-columns-in-a-table
LM_Part{alpha, beta, gamma}.Properties.VariableNames(1) = {'ID_States'};
clear opts

% Now we do essentially the same thing to calculate efficiency losses for 
% individuals who are deceased. This is necessary because losses depend not
% simply on the individual's health state, but also on their individual 
% characteristics - because this determines what they WOULD have earned had
% they still been alive.
opts = detectImportOptions(strcat(MAT_Path, "Simulation_Results\",...
    "Efficiency_Loss_Calculations\", "Eff_Losses_Dead_", Case, ".dat"));
opts = setvartype(opts,{'Age_1'},'char');
opts = setvartype(opts,{'Net_Emp_Inc_FT_00'},'char');
getvaropts(opts,{'Age_1','Net_Emp_Inc_FT_00'});
opts.SelectedVariableNames = {'Age_1', 'Net_Emp_Inc_FT_00'};
EL_Dead_00{alpha, beta, gamma} = readtable(strcat(MAT_Path, ...
    "Simulation_Results\", "Efficiency_Loss_Calculations\", ...
    "Eff_Losses_Dead_", Case, ".dat"), opts);
clear opts

% NOTE 032219: 
rwd_opt_temp = load(strcat(MAT_Path, "Simulation_Results\", ...
    "rwdstar_opt_", Case, ".mat"));
RWD_Opt{alpha, beta, gamma} = [(1:size(RefStates, 1))' ...
    rwd_opt_temp.rwdstar_opt];

rwd_alt_A_temp = load(strcat(MAT_Path, "Simulation_Results\", ...
    "rwdstar_alt_", Case, "_1", ".mat"));
RWD_Alt_A{alpha, beta, gamma} = [(1:size(RefStates, 1))' ...
    rwd_alt_A_temp.rwdstar_alt_temp];

rwd_alt_DA_temp = load(strcat(MAT_Path, "Simulation_Results\", ...
    "rwdstar_alt_", Case, "_8", ".mat"));
RWD_Alt_DA{alpha, beta, gamma} = [(1:size(RefStates, 1))' ...
    rwd_alt_DA_temp.rwdstar_alt_temp];

rwd_alt_HB_temp = load(strcat(MAT_Path, "Simulation_Results\", ...
    "rwdstar_alt_", Case, "_hybrid", ".mat"));
RWD_Alt_HB{alpha, beta, gamma} = [(1:size(RefStates, 1))' ...
    rwd_alt_HB_temp.rwdstar_alt_temp];

       end
   end
end

% 100418: By default, SI_opt and SI_alt exist as tables. Some of the
% computations we are about to perform cannot, however, be done with a 
% table, so we convert these to arrays first.
SR_Opt_Cohort = table2array(SR_Opt_Cohort); 
SR_Alt_A_Cohort = table2array(SR_Alt_A_Cohort);
SR_Alt_DA_Cohort = table2array(SR_Alt_DA_Cohort);
SR_Alt_HB_Cohort = table2array(SR_Alt_HB_Cohort);

for alpha = R_alpha
   for beta = R_beta
       for gamma = R_gamma
           
% Assign the table characterizing efficiency losses for an individual with 
% a certain set of characteristics to a temporary "host".
           EL_temp = EL_Alive{alpha, beta, gamma};
% Load data related to labour market participation from the optimal policy.
           LM_Part_Temp = LM_Part{alpha, beta, gamma};
% Link information on LMP with information about efficiency losses for
% living members of the cohort, using the common key ID_States.
% https://www.mathworks.com/help/matlab/ref/table.join.html
           EL_temp = join(EL_temp, LM_Part_Temp, 'Keys', 'ID_States');
% Identify the states in which the individual who chosen to retire.
           Rows = EL_temp.LMP == "RT";
% If the individual chooses to retire, we assume there are no efficiency 
% losses during that year or thereafter.
           EL_temp{EL_temp.LMP == "RT", 12:14} = 0;
% Store the table headers in a cell array so we can reference them later.
           EL_temp_names = EL_temp.Properties.VariableNames;
           
% For every column of EL_temp ...
           for delta = 1:size(EL_temp, 2)
% ... if that column is classified as numeric data type ...
               if isnumeric(EL_temp.(EL_temp_names{delta}))
% ... convert data type to cell.
                   EL_temp.(EL_temp_names{delta}) = ...
                       num2cell(EL_temp.(EL_temp_names{delta}));
               end
           end
           
% We only require a subset of the variables from Eff_Losses_Alive, namely 
% the identifier, and the three variables that record efficiency losses
% corresponding to each state (i.e., losses attributable to absenteeism and
% presenteeism, and total losses).
           EL_Alive{alpha, beta, gamma} = ...
               table2array(EL_temp(:, [1, 12:14]));
% Convert the contents of Sim_Results_Opt and Sim_Results_Alt (which 
% consist of tables) to cell arrays. This way, we'll be able to merge data 
% presently stored in different constructs.
           SR_Opt{alpha, beta, gamma} = ...
               table2array(SR_Opt{alpha, beta, gamma});
           SR_Alt_A{alpha, beta, gamma} = ...
               table2array(SR_Alt_A{alpha, beta, gamma});
           SR_Alt_DA{alpha, beta, gamma} = ...
               table2array(SR_Alt_DA{alpha, beta, gamma});
           SR_Alt_HB{alpha, beta, gamma} = ...
               table2array(SR_Alt_HB{alpha, beta, gamma});
           
       end
   end
end

% NOTE 032219: Added rewards
% Duplicate the contents of Sim_Results_Opt and Sim_Results_Alt. This is
% done because different sets of computations can be performed using the
% same information.
[SR_Opt_EL_Abs, SR_Opt_EL_Pre, SR_Opt_EL_Tot, SR_Opt_HS, SR_Opt_Age, ...
    SR_Opt_LMP, SR_Opt_RWD] = deal(SR_Opt);
[SR_Alt_EL_Abs_A, SR_Alt_EL_Pre_A, SR_Alt_EL_Tot_A, SR_Alt_HS_A, ...
    SR_Alt_Age_A, SR_Alt_LMP_A, SR_Alt_RWD_A] = deal(SR_Alt_A);
[SR_Alt_EL_Abs_DA, SR_Alt_EL_Pre_DA, SR_Alt_EL_Tot_DA, SR_Alt_HS_DA, ...
    SR_Alt_Age_DA, SR_Alt_LMP_DA, SR_Alt_RWD_DA] = deal(SR_Alt_DA);
[SR_Alt_EL_Abs_HB, SR_Alt_EL_Pre_HB, SR_Alt_EL_Tot_HB, SR_Alt_HS_HB, ...
    SR_Alt_Age_HB, SR_Alt_LMP_HB, SR_Alt_RWD_HB] = deal(SR_Alt_HB);

% NOTE 032219: 
SR_Opt_RWD_disc = cell(2, 3, 4); SR_Alt_RWD_A_disc = cell(2, 3, 4);
SR_Alt_RWD_DA_disc = cell(2, 3, 4); SR_Alt_RWD_HB_disc = cell(2, 3, 4);

% NOTE 032219:
Discounts = table((0:T)', 'VariableNames', {'Elapsed_Time'}); 
DF = 0.95; DR = (1/DF) - 1;
Discounts.DiscFact = 1./power((1+DR), Discounts.Elapsed_Time);
Discounts_Matrix = diag(Discounts.DiscFact);

% 100418: This is probably the most unusual line in this script. The ideal
% would be to create a new table that inserts the information we're looking
% for in the same place as it appeared in the original table (i.e.,
% SimResults). The issue with the next line of code is that results are 
% presented in a single columns which we must then reshape. It is not a 
% huge inconvenience, but it is puzzling that there is no convenient
% workaround. I attempted a couple of things which did not work well:
% 1) The same operation can be accomplished with a for-loop, but it takes
% substantially longer. By contrast, the approach here appears to benefit
% from vectorization.
% 2) I also explored the possibility of declaring the desired size of
% SimRes_HS_00. When I attempt this, I get the error message, "cannot 
% create a table variable with a discontiguous index."
% 101518: Calculating efficiency losses is a little more involved, because
% the magnitude of the losses depends on the individual's characteristics.
% Hence, we'll need to reference the contents of Eff_Losses_Alive and
% Eff_Losses_Dead to proceed.
for alpha = R_alpha
   for beta = R_beta
       for gamma = R_gamma

% NOTE 112418: Caught error in code, in which Opt_Temp was referred to
% later in each loop but not defined here. This did not result in
% noticeable issues because we were dealing with only a small set of the
% cohort - but it would have become problematic if we had been dealing with
% two or more subsets of the population. The error has now been corrected.
% Load data related to labour market participation from the optimal policy.
           LM_Part_Temp = LM_Part{alpha, beta, gamma};
% Extract the information we're looking for from each of these three cell 
% arrays, and assign it to a temporary construct that we can recycle on 
% each loop.
           EL_A_Temp = EL_Alive{alpha, beta, gamma};
           EL_D_Temp = table2array(EL_Dead_00{alpha, beta, gamma});
           SR_Opt_Temp = SR_Opt{alpha, beta, gamma};
           SR_Alt_Temp_A = SR_Alt_A{alpha, beta, gamma};
           SR_Alt_Temp_DA = SR_Alt_DA{alpha, beta, gamma};
           SR_Alt_Temp_HB = SR_Alt_HB{alpha, beta, gamma};
% NOTE 032219: 
           RWD_Opt_Temp = RWD_Opt{alpha, beta, gamma};
           RWD_Alt_Temp_A = RWD_Alt_A{alpha, beta, gamma};
           RWD_Alt_Temp_DA = RWD_Alt_DA{alpha, beta, gamma};
           RWD_Alt_Temp_HB = RWD_Alt_HB{alpha, beta, gamma};
           NumRow = size(SR_Opt_Temp(:, 4:size(SR_Opt_Temp, 2)), 1);
           NumCol = size(SR_Opt_Temp(:, 4:size(SR_Opt_Temp, 2)), 2);

% Create three columns containing the individual's characteristics, and 
% store it in the appropriate slot in the Sim_Results_ID cell array.
           SR_ID{alpha, beta, gamma} = ...
               horzcat(...
               repmat(alpha, size(SR_Opt_Temp, 1), 1), ...
               repmat(beta, size(SR_Opt_Temp, 1), 1), ...
               repmat(gamma, size(SR_Opt_Temp, 1), 1));
           
           EL_Dead_01{alpha, beta, gamma} = ...
               horzcat(EL_D_Temp, ...
               num2cell(repmat(alpha, size(EL_D_Temp, 1), 1)), ...
               num2cell(repmat(beta, size(EL_D_Temp, 1), 1)), ...
               num2cell(repmat(gamma, size(EL_D_Temp, 1), 1)));
           
% 101918: Note that this refers to TOTAL efficiency losses. We'll need to 
% replicate this exercise to decompose total efficiency losses into losses 
% atributable to absenteeism, presenteeism, and premature mortality.
           SR_Opt_EL_Abs{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 2), NumRow, NumCol);
           SR_Alt_EL_Abs_A{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 2), NumRow, NumCol);
           SR_Alt_EL_Abs_DA{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 2), NumRow, NumCol);     
           SR_Alt_EL_Abs_HB{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 2), NumRow, NumCol);
                      
           SR_Opt_EL_Pre{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 3), NumRow, NumCol);
           SR_Alt_EL_Pre_A{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 3), NumRow, NumCol);
           SR_Alt_EL_Pre_DA{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 3), NumRow, NumCol);
           SR_Alt_EL_Pre_HB{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 3), NumRow, NumCol);
           
           SR_Opt_EL_Tot{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 4), NumRow, NumCol);
           SR_Alt_EL_Tot_A{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 4), NumRow, NumCol);
           SR_Alt_EL_Tot_DA{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 4), NumRow, NumCol);
           SR_Alt_EL_Tot_HB{alpha, beta, gamma} = ...
               reshape(EL_A_Temp(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 4), NumRow, NumCol);

% This is very similar to the code presented immediately above, except that
% we are using a different lookup table - Eff_Losses_Temp doesn't include 
% any other information about the individual's health state, but RefStates 
% does. Note, however, that we convert the information extracted from the 
% RefStates table to an array before we reshape it.
           SR_Opt_HS{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 4)), NumRow, NumCol);
           SR_Alt_HS_A{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 4)), NumRow, NumCol);
           SR_Alt_HS_DA{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 4)), NumRow, NumCol);
           SR_Alt_HS_HB{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 4)), NumRow, NumCol);

% The next piece of information we need is the individual's age. This is 
% important because efficiency losses attributable to premature mortality 
% depends upon the age an individual would have been had they still been 
% alive. Curiously, though, we won't actually use the variable to 
% accomplish this, but will rather derive this information from the state 
% ID. As such, the variable will serve primarily to faciliate validation of 
% the results.
           SR_Opt_Age{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 2)), NumRow, NumCol);
           SR_Alt_Age_A{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 2)), NumRow, NumCol);
           SR_Alt_Age_DA{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 2)), NumRow, NumCol);
           SR_Alt_Age_HB{alpha, beta, gamma} = ...
               reshape(table2array(RefStates(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 2)), NumRow, NumCol);
% The last thing we need to know is the individual's labour market choices.
% This will tell us when the individual retired, and whether retirement
% preceded death. This is important, in that it will enable us to ascertain
% whether the individual would likely have continued to participate in the
% labour market had they not died - which, in turn, facilitates calculation
% of economic losses attributable to premature mortality.
           SR_Opt_LMP{alpha, beta, gamma} = ...
               reshape(table2array(LM_Part_Temp(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 2)), NumRow, NumCol);
           SR_Alt_LMP_A{alpha, beta, gamma} = ...
               reshape(table2array(LM_Part_Temp(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 2)), NumRow, NumCol);
           SR_Alt_LMP_DA{alpha, beta, gamma} = ...
               reshape(table2array(LM_Part_Temp(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 2)), NumRow, NumCol);
           SR_Alt_LMP_HB{alpha, beta, gamma} = ...
               reshape(table2array(LM_Part_Temp(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 2)), NumRow, NumCol);

% NOTE 032219: Calculate rewards enjoyed by simulated agents.
           SR_Opt_RWD{alpha, beta, gamma} = ...
               reshape(RWD_Opt_Temp(SR_Opt_Temp ...
               (:, 4:size(SR_Opt_Temp, 2)), 2), NumRow, NumCol);
           SR_Alt_RWD_A{alpha, beta, gamma} = ...
               reshape(RWD_Alt_Temp_A(SR_Alt_Temp_A ...
               (:, 4:size(SR_Alt_Temp_A, 2)), 2), NumRow, NumCol);
           SR_Alt_RWD_DA{alpha, beta, gamma} = ...
               reshape(RWD_Alt_Temp_DA(SR_Alt_Temp_DA ...
               (:, 4:size(SR_Alt_Temp_DA, 2)), 2), NumRow, NumCol);
           SR_Alt_RWD_HB{alpha, beta, gamma} = ...
               reshape(RWD_Alt_Temp_HB(SR_Alt_Temp_HB ...
               (:, 4:size(SR_Alt_Temp_HB, 2)), 2), NumRow, NumCol);
           
           SR_Opt_RWD_disc{alpha, beta, gamma} = ...
               SR_Opt_RWD{alpha, beta, gamma} * Discounts_Matrix;
           SR_Alt_RWD_A_disc{alpha, beta, gamma} = ...
               SR_Alt_RWD_A{alpha, beta, gamma} * Discounts_Matrix;
           SR_Alt_RWD_DA_disc{alpha, beta, gamma} = ...
               SR_Alt_RWD_DA{alpha, beta, gamma} * Discounts_Matrix;
           SR_Alt_RWD_HB_disc{alpha, beta, gamma} = ...
               SR_Alt_RWD_HB{alpha, beta, gamma} * Discounts_Matrix;
           
       end
   end
end

%% Identifying Instances of Premature Mortality

% 112418: This content is new. It is intended to enable the simulation
% model to estimate economic losses resulting from premature mortality,
% which we define as mortality occurring before the average age of
% retirement. This is actually not a straightforward exercise, because it
% requires us to combine information on health status with information on
% labour market participation.
for alpha = R_alpha
   for beta = R_beta
       for gamma = R_gamma
           
% Draw pertinent information related to health status and labour market
% participation.
            SR_Ret_opt = SR_Opt_LMP{alpha, beta, gamma};
            SR_HS_opt = SR_Opt_HS{alpha, beta, gamma};
            SR_Ret_alt_A = SR_Alt_LMP_A{alpha, beta, gamma};
            SR_HS_alt_A = SR_Alt_HS_A{alpha, beta, gamma};
            SR_Ret_alt_DA = SR_Alt_LMP_DA{alpha, beta, gamma};
            SR_HS_alt_DA = SR_Alt_HS_DA{alpha, beta, gamma};
            SR_Ret_alt_HB = SR_Alt_LMP_HB{alpha, beta, gamma};
            SR_HS_alt_HB = SR_Alt_HS_HB{alpha, beta, gamma};

% We track changes in health status and labour market participation by
% generating two sets of sub-matrices - one from the starting state to the
% second-to-last period, and one from the second period to the final
% period.
            SR_Ret_opt_01 = SR_Ret_opt(:, 1:(end-1));
            SR_Ret_opt_02 = SR_Ret_opt(:, 2:end);
            SR_HS_opt_01 = SR_HS_opt(:, 1:(end-1));
            SR_HS_opt_02 = SR_HS_opt(:, 2:end);

% Retirement is said to occur when an individual is working full- or
% part-time in one period, and retires in the next (recalling that
% retirement in the model is a choice, and not a state).
            SR_Ret_opt_03 = 1.*(or(SR_Ret_opt_01 == "FT",...
                SR_Ret_opt_01 == "FT") & SR_Ret_opt_02 == "RT");
% Mortality occurs when an individual is healthy or diabetic in one period,
% and dies in the next.
            SR_HS_opt_03 = 1.*(and(or(SR_HS_opt_01 == "Healthy", ...
                SR_HS_opt_01 == "T2DM"), SR_HS_opt_02 == "Dead"));
% 121618: Mortality attributable to T2DM occurs when an individual is
% diabetic in one period and subsequently dies. This is important to track
% because economic losses resulting from premature mortality will 
% ultimately be defined in terms of death among individuals who have been
% diagnosed with T2DM (while not trivial, this study does not focus
% productivity losses associated with mortality that results from other
% acute or chronic health conditions or from natural aging).
            SR_HS_opt_03DM = 1.*(and(SR_HS_opt_01 == "T2DM", ...
                SR_HS_opt_02 == "Dead"));

% Do the same thing for the alternate policies.
            SR_Ret_alt_A_01 = SR_Ret_alt_A(:, 1:(end-1));
            SR_Ret_alt_A_02 = SR_Ret_alt_A(:, 2:end);
            SR_HS_alt_A_01 = SR_HS_alt_A(:, 1:(end-1));
            SR_HS_alt_A_02 = SR_HS_alt_A(:, 2:end);
            
            SR_Ret_alt_A_03 = 1.*(or(SR_Ret_opt_01 == "FT",...
                SR_Ret_opt_01 == "FT") & SR_Ret_opt_02 == "RT");
            SR_HS_alt_A_03 = 1.*(and(or(SR_HS_alt_A_01 == "Healthy", ...
                SR_HS_alt_A_01 == "T2DM"), SR_HS_alt_A_02 == "Dead"));
            SR_HS_alt_A_03DM = 1.*(and(SR_HS_alt_A_01 == "T2DM", ...
                SR_HS_alt_A_02 == "Dead"));
            
            SR_Ret_alt_DA_01 = SR_Ret_alt_DA(:, 1:(end-1));
            SR_Ret_alt_DA_02 = SR_Ret_alt_DA(:, 2:end);
            SR_HS_alt_DA_01 = SR_HS_alt_DA(:, 1:(end-1));
            SR_HS_alt_DA_02 = SR_HS_alt_DA(:, 2:end);
            
            SR_Ret_alt_DA_03 = 1.*(or(SR_Ret_opt_01 == "FT",...
                SR_Ret_opt_01 == "FT") & SR_Ret_opt_02 == "RT");
            SR_HS_alt_DA_03 = 1.*(and(or(SR_HS_alt_DA_01 == "Healthy", ...
                SR_HS_alt_DA_01 == "T2DM"), SR_HS_alt_DA_02 == "Dead"));
            SR_HS_alt_DA_03DM = 1.*(and(SR_HS_alt_DA_01 == "T2DM", ...
                SR_HS_alt_DA_02 == "Dead"));
            
            SR_Ret_alt_HB_01 = SR_Ret_alt_HB(:, 1:(end-1));
            SR_Ret_alt_HB_02 = SR_Ret_alt_HB(:, 2:end);
            SR_HS_alt_HB_01 = SR_HS_alt_HB(:, 1:(end-1));
            SR_HS_alt_HB_02 = SR_HS_alt_HB(:, 2:end);
            
            SR_Ret_alt_HB_03 = 1.*(or(SR_Ret_opt_01 == "FT",...
                SR_Ret_opt_01 == "FT") & SR_Ret_opt_02 == "RT");
            SR_HS_alt_HB_03 = 1.*(and(or(SR_HS_alt_HB_01 == "Healthy", ...
                SR_HS_alt_HB_01 == "T2DM"), SR_HS_alt_HB_02 == "Dead"));
            SR_HS_alt_HB_03DM = 1.*(and(SR_HS_alt_HB_01 == "T2DM", ...
                SR_HS_alt_HB_02 == "Dead"));
            
% This next step is central to our approach. People who die in the model
% are classified as having retired (since it wouldn't make sense to
% categorize them as working full-time or part-time). However, this creates
% a problem, because we now do not know whether a person who "retires" does
% so of their own volition or not (overlooking the fact that people who are
% unemployed upon reaching the standard age of retirement, or thereafter,
% cannot strictly be said to have CHOSEN to retire in this model). To 
% address this, we create a new array, the cells of which must meet the 
% following conditions:
% a) The individual opts to retire in the current period, after having
% worked in the previous period.
% b) AND, the individual DID NOT experience mortality between the preceding
% period and this one.
            SR_Ret_opt_04 = 1.*(and(SR_Ret_opt_03 == 1, ...
                SR_HS_opt_03 == 0));
            SR_Ret_alt_A_04 = 1.*(and(SR_Ret_alt_A_03 == 1, ...
                SR_HS_alt_A_03 == 0));
            SR_Ret_alt_DA_04 = 1.*(and(SR_Ret_alt_DA_03 == 1, ...
                SR_HS_alt_DA_03 == 0));
            SR_Ret_alt_HB_04 = 1.*(and(SR_Ret_alt_HB_03 == 1, ...
                SR_HS_alt_HB_03 == 0));
% Rows that do not contain a non-zero value denote individuals who died 
% after having retired. We now use the FIND function to pinpoint which
% individuals these were, noting that the first and second values on the 
% LHS refer to the row and column indices associated with each non-zero
% value in the array.
            [a_opt, b_opt] = find(SR_Ret_opt_04);
            [a_alt_A, b_alt_A] = find(SR_Ret_alt_A_04);
            [a_alt_DA, b_alt_DA] = find(SR_Ret_alt_DA_04);
            [a_alt_HB, b_alt_HB] = find(SR_Ret_alt_HB_04);
% Finally, we can calculate the average age of retirement, excluding those
% who 'retired' because they experienced mortality, by calculating the mean
% of the column "coordinate" of each value, which denotes the year in which
% retirement occurs.
            Ret_Exp_Opt = round(mean(b_opt));
            Ret_Exp_Alt_A = round(mean(b_alt_A));
            Ret_Exp_Alt_DA = round(mean(b_alt_DA));
            Ret_Exp_Alt_HB = round(mean(b_alt_HB));

% Now identify rows in which the transition to retirement actually
% signifies the individual's death. This is done by changing a single value
% in the statement above (i.e., SR_HS_opt_03/SR_HS_alt_03 changes from 0 to
% 1). This time, we are identifying the row and column indexes denoting
% cases where an individual dies before choosing to leave the labour market
% of their own volition - which, as noted above, is how we have chosen to
% define premature mortality.
% 
% 121618: These next few lines of code can be implemented in at least two
% ways. Originally, the statements below referred to the variables
% SR_HS_opt_03/SR_HS_alt_A_03/SR_HS_alt_DA_03; however, I subsequently
% determined this was not appropriate for the purposes of the analysis
% because it included premature mortality among people who were healthy
% when they died, which had the effect of increasing estimated productivity
% losses stemming from this source. This new formulation instead focuses
% specifically on death among individuals already diagnosed with T2DM.
            SR_Ret_opt_05 = 1.*(and(SR_Ret_opt_03 == 1, ...
                SR_HS_opt_03DM == 1));
            [c_opt, d_opt] = find(SR_Ret_opt_05);
            Mort_Opt = [c_opt d_opt]; Mort_Opt = sortrows(Mort_Opt);
            
            SR_Ret_alt_A_05 = 1.*(and(SR_Ret_alt_A_03 == 1, ...
                SR_HS_alt_A_03DM == 1));
            [c_alt_A, d_alt_A] = find(SR_Ret_alt_A_05);
            Mort_Alt_A = [c_alt_A d_alt_A]; 
            Mort_Alt_A = sortrows(Mort_Alt_A);
            
            SR_Ret_alt_DA_05 = 1.*(and(SR_Ret_alt_DA_03 == 1, ...
                SR_HS_alt_DA_03DM == 1));
            [c_alt_DA, d_alt_DA] = find(SR_Ret_alt_DA_05);
            Mort_Alt_DA = [c_alt_DA d_alt_DA]; 
            Mort_Alt_DA = sortrows(Mort_Alt_DA);
            
            SR_Ret_alt_HB_05 = 1.*(and(SR_Ret_alt_HB_03 == 1, ...
                SR_HS_alt_HB_03DM == 1));
            [c_alt_HB, d_alt_HB] = find(SR_Ret_alt_HB_05);
            Mort_Alt_HB = [c_alt_HB d_alt_HB]; 
            Mort_Alt_HB = sortrows(Mort_Alt_HB);
% Create matrices to denote periods in which we will be tracking economic
% losses attributable to premature mortality.
            SR_PM_Opt_Temp = zeros(size(SR_Ret_opt, 1), T);
            SR_PM_Alt_Temp_A = zeros(size(SR_Ret_alt_A, 1), T);
            SR_PM_Alt_Temp_DA = zeros(size(SR_Ret_alt_DA, 1), T);
            SR_PM_Alt_Temp_HB = zeros(size(SR_Ret_alt_HB, 1), T);
            
% Periods "marked" by premature mortality are those including every year
% between when the person is initially classified as deceased, and when
% they would have been expected to retire had they still been alive.
            for i = 1:size(Mort_Opt,1)
                SR_PM_Opt_Temp(Mort_Opt(i,1), ...
                    Mort_Opt(i,2):Ret_Exp_Opt) = 1;
            end
% Now we do the same thing for the alternate policy. Note that we cannot
% use the same for-loop to carry out this operation for the optimal and
% alternate policies, as we cannot be sure the average age of retirement
% will be equivalent.
            for i = 1:size(Mort_Alt_A,1)
                SR_PM_Alt_Temp_A(Mort_Alt_A(i,1), ...
                    Mort_Alt_A(i,2):Ret_Exp_Alt_A) = 1;
            end
            
            for i = 1:size(Mort_Alt_DA,1)
                SR_PM_Alt_Temp_DA(Mort_Alt_DA(i,1), ...
                    Mort_Alt_DA(i,2):Ret_Exp_Alt_DA) = 1;
            end
            
            for i = 1:size(Mort_Alt_HB,1)
                SR_PM_Alt_Temp_HB(Mort_Alt_HB(i,1), ...
                    Mort_Alt_HB(i,2):Ret_Exp_Alt_HB) = 1;
            end

% We've got one more hurdle to overcome. Strictly speaking, the simulation
% model consists of 71 periods, not 70 periods, since we must include "Year
% 0", which denotes the state of the cohort upon entering the model. To do
% this, we create a column of zero values with as many rows as there are
% people in each segment of the cohort, and append it to the "front" of the
% array (if we do not do this, the reshape command in the next section will
% return an error).
            SR_PM_Opt_Temp = ...
                [zeros(size(SR_Ret_opt, 1), 1) SR_PM_Opt_Temp];
            SR_PM_Alt_Temp_A = ...
                [zeros(size(SR_Ret_alt_A, 1), 1) SR_PM_Alt_Temp_A];
            SR_PM_Alt_Temp_DA = ...
                [zeros(size(SR_Ret_alt_DA, 1), 1) SR_PM_Alt_Temp_DA];
            SR_PM_Alt_Temp_HB = ...
                [zeros(size(SR_Ret_alt_HB, 1), 1) SR_PM_Alt_Temp_HB];

% Now, convert the temporary array from numeric data type to cell (to
% preserve consistency with the other data we integrate this with below),
% and store it in the proper place in the premature mortality array.
            SR_PM_Opt{alpha, beta, gamma} = num2cell(SR_PM_Opt_Temp);
            SR_PM_Alt_A{alpha, beta, gamma} = num2cell(SR_PM_Alt_Temp_A);
            SR_PM_Alt_DA{alpha, beta, gamma} = num2cell(SR_PM_Alt_Temp_DA);
            SR_PM_Alt_HB{alpha, beta, gamma} = num2cell(SR_PM_Alt_Temp_HB);
            
       end
   end
end

%% Calculating Efficiency Losses

Size_SimRes = size(SR_Opt_Cohort, 1)*...
    size(SR_Opt_Cohort(:, 4:size(SR_Opt_Cohort, 2)), 2);
% 102018: This code doesn't work - when run, I obtain the following error:
% To RESHAPE the number of elements must not change.
% See "~" in scrapyard - START
Sim_Res = {...
    SR_Opt_EL_Abs; SR_Alt_EL_Abs_A; SR_Alt_EL_Abs_DA; SR_Alt_EL_Abs_HB; ...
    SR_Opt_EL_Pre; SR_Alt_EL_Pre_A; SR_Alt_EL_Pre_DA; SR_Alt_EL_Pre_HB; ... 
    SR_Opt_EL_Tot; SR_Alt_EL_Tot_A; SR_Alt_EL_Tot_DA; SR_Alt_EL_Tot_HB; ...
    SR_Opt_HS;     SR_Alt_HS_A;     SR_Alt_HS_DA;     SR_Alt_HS_HB; ...
    SR_Opt_Age;    SR_Alt_Age_A;    SR_Alt_Age_DA;    SR_Alt_Age_HB; ... 
    SR_PM_Opt;     SR_PM_Alt_A;     SR_PM_Alt_DA;     SR_PM_Alt_HB};
func_sim = @(x) reshape(vertcat(x{:,:,:}), Size_SimRes, 1);
Sim_Res = cellfun(func_sim, Sim_Res, 'UniformOutput', false);

% https://stackoverflow.com/questions/2337126/how-do-i-do-multiple-assignment-in-matlab
   [SR_opt_EL_Abs_00, SR_alt_EL_Abs_A_00, SR_alt_EL_Abs_DA_00, SR_alt_EL_Abs_HB_00, ...
    SR_opt_EL_Pre_00, SR_alt_EL_Pre_A_00, SR_alt_EL_Pre_DA_00, SR_alt_EL_Pre_HB_00, ...
    SR_opt_EL_Tot_00, SR_alt_EL_Tot_A_00, SR_alt_EL_Tot_DA_00, SR_alt_EL_Tot_HB_00,...
    SR_opt_HS_00,     SR_alt_HS_A_00,     SR_alt_HS_DA_00,     SR_alt_HS_HB_00, ...
    SR_opt_Age_00,    SR_alt_Age_A_00,    SR_alt_Age_DA_00,    SR_alt_Age_HB_00, ...
    SR_opt_PM_00,     SR_alt_PM_A_00,     SR_alt_PM_DA_00,     SR_alt_PM_HB_00]...
    = deal(Sim_Res{:, 1});
% See "~" in scrapyard - END

EL_Dead_02 = vertcat(EL_Dead_01{:,:,:});

SR_ID = vertcat(SR_ID{:,:,:});
SR_ID_01 = num2cell(repmat(SR_ID, ...
    size(SR_opt_EL_Tot_00, 1)/size(SR_ID, 1), 1));
SR_opt_EL_01 = ...
    cell2table([SR_ID_01 SR_opt_Age_00 SR_opt_HS_00 ...
    SR_opt_EL_Abs_00 SR_opt_EL_Pre_00 SR_opt_EL_Tot_00]);
SR_alt_EL_A_01 = ...
    cell2table([SR_ID_01 SR_alt_Age_A_00 SR_alt_HS_A_00 ...
    SR_alt_EL_Abs_A_00 SR_alt_EL_Pre_A_00 SR_alt_EL_Tot_A_00]);
SR_alt_EL_DA_01 = ...
    cell2table([SR_ID_01 SR_alt_Age_DA_00 SR_alt_HS_DA_00 ...
    SR_alt_EL_Abs_DA_00 SR_alt_EL_Pre_DA_00 SR_alt_EL_Tot_DA_00]);
SR_alt_EL_HB_01 = ...
    cell2table([SR_ID_01 SR_alt_Age_HB_00 SR_alt_HS_HB_00 ...
    SR_alt_EL_Abs_HB_00 SR_alt_EL_Pre_HB_00 SR_alt_EL_Tot_HB_00]);
EL_Dead_02 = cell2table(EL_Dead_02);

% Need to add a new column; furthermore, we require a label to distinguish
% the new variable from the other variables included in each table.
Labels = {'Index'};
Index_opt = table(transpose(1:size(SR_opt_EL_01, 1)));
Index_opt.Properties.VariableNames = Labels;
Index_alt_A = table(transpose(1:size(SR_alt_EL_A_01, 1)));
Index_alt_A.Properties.VariableNames = Labels;
Index_alt_DA = table(transpose(1:size(SR_alt_EL_DA_01, 1)));
Index_alt_DA.Properties.VariableNames = Labels;
Index_alt_HB = table(transpose(1:size(SR_alt_EL_HB_01, 1)));
Index_alt_HB.Properties.VariableNames = Labels;

SR_opt_EL_01 = [SR_opt_EL_01 Index_opt];
SR_alt_EL_A_01 = [SR_alt_EL_A_01 Index_alt_A];
SR_alt_EL_DA_01 = [SR_alt_EL_DA_01 Index_alt_DA];
SR_alt_EL_HB_01 = [SR_alt_EL_HB_01 Index_alt_HB];

Labels_SR = {'Gender', 'Ethnicity', 'Education', 'Age_00', 'HS', ...
    'EL_Abs', 'EL_Pre', 'EL_Tot', 'Index'};
SR_opt_EL_01.Properties.VariableNames = Labels_SR;
SR_alt_EL_A_01.Properties.VariableNames = Labels_SR;
SR_alt_EL_DA_01.Properties.VariableNames = Labels_SR;
SR_alt_EL_HB_01.Properties.VariableNames = Labels_SR;
Labels_ELD = {'Age_01', 'EL_Dec', 'Gender', 'Ethnicity', ...
    'Education'};
EL_Dead_02.Properties.VariableNames = Labels_ELD;

EL_Dead_02.EL_Dec = str2double(string(EL_Dead_02.EL_Dec));
EL_Dead_02.Age_01 = str2double(string(EL_Dead_02.Age_01));
SR_opt_EL_01.Age_00 = str2double(string(SR_opt_EL_01.Age_00));
SR_alt_EL_A_01.Age_00 = str2double(string(SR_alt_EL_A_01.Age_00));
SR_alt_EL_DA_01.Age_00 = str2double(string(SR_alt_EL_DA_01.Age_00));
SR_alt_EL_HB_01.Age_00 = str2double(string(SR_alt_EL_HB_01.Age_00));

SR_opt_EL_01.Age_01 = ...
    25 + floor(SR_opt_EL_01.Index./(size(SR_Opt_Cohort, 1)));
SR_opt_EL_01.Mod = mod(SR_opt_EL_01.Index, size(SR_Opt_Cohort,1)) == 0;
SR_opt_EL_01.Age_01 = SR_opt_EL_01.Age_01 - ...
    SR_opt_EL_01.Mod;
SR_opt_EL_01.Mod = [];

SR_alt_EL_A_01.Age_01 = ...
    25 + floor(SR_alt_EL_A_01.Index./(size(SR_Alt_A_Cohort, 1)));
SR_alt_EL_A_01.Mod = mod(SR_alt_EL_A_01.Index, size(SR_Alt_A_Cohort,1)) == 0;
SR_alt_EL_A_01.Age_01 = SR_alt_EL_A_01.Age_01 - ...
    SR_alt_EL_A_01.Mod;
SR_alt_EL_A_01.Mod = []; 

SR_alt_EL_DA_01.Age_01 = ...
    25 + floor(SR_alt_EL_DA_01.Index./(size(SR_Alt_DA_Cohort, 1)));
SR_alt_EL_DA_01.Mod = mod(SR_alt_EL_DA_01.Index, ...
    size(SR_Alt_DA_Cohort,1)) == 0;
SR_alt_EL_DA_01.Age_01 = SR_alt_EL_DA_01.Age_01 - ...
    SR_alt_EL_DA_01.Mod;
SR_alt_EL_DA_01.Mod = [];

SR_alt_EL_HB_01.Age_01 = ...
    25 + floor(SR_alt_EL_HB_01.Index./(size(SR_Alt_HB_Cohort, 1)));
SR_alt_EL_HB_01.Mod = mod(SR_alt_EL_HB_01.Index, ...
    size(SR_Alt_HB_Cohort,1)) == 0;
SR_alt_EL_HB_01.Age_01 = SR_alt_EL_HB_01.Age_01 - ...
    SR_alt_EL_HB_01.Mod;
SR_alt_EL_HB_01.Mod = []; 

% https://www.mathworks.com/help/matlab/ref/string.html
SR_opt_EL_01.HS = string(SR_opt_EL_01.HS);
SR_alt_EL_A_01.HS = string(SR_alt_EL_A_01.HS);
SR_alt_EL_DA_01.HS = string(SR_alt_EL_DA_01.HS);
SR_alt_EL_HB_01.HS = string(SR_alt_EL_HB_01.HS);
% SimRes_opt_EL_01.Index = []; SimRes_alt_EL_01.Index = [];

SR_opt_EL_02 = outerjoin(SR_opt_EL_01, EL_Dead_02, ...
    'Type', 'Left', 'MergeKeys', true);
SR_alt_EL_A_02 = outerjoin(SR_alt_EL_A_01, EL_Dead_02, ...
    'Type', 'Left', 'MergeKeys', true);
SR_alt_EL_DA_02 = outerjoin(SR_alt_EL_DA_01, EL_Dead_02, ...
    'Type', 'Left', 'MergeKeys', true);
SR_alt_EL_HB_02 = outerjoin(SR_alt_EL_HB_01, EL_Dead_02, ...
    'Type', 'Left', 'MergeKeys', true);

Full_Loss_opt = table(SR_opt_EL_02.HS == "Dead");
Full_Loss_alt_A = table(SR_alt_EL_A_02.HS == "Dead");
Full_Loss_alt_DA = table(SR_alt_EL_DA_02.HS == "Dead");
Full_Loss_alt_HB = table(SR_alt_EL_HB_02.HS == "Dead");

SR_ID_02 = cell2table([SR_ID_01 SR_opt_Age_00]);
Labels_SR_PM = {'Gender', 'Ethnicity', 'Education', 'Age_00'};
SR_ID_02.Properties.VariableNames = Labels_SR_PM;
SR_ID_02 = [SR_ID_02 Index_opt];
SR_ID_02.Age_00 = str2double(string(SR_ID_02.Age_00));

SR_ID_02.Age_01 = 25 + floor(SR_ID_02.Index./(size(SR_Opt_Cohort, 1)));
SR_ID_02.Mod = mod(SR_ID_02.Index, size(SR_Opt_Cohort,1)) == 0;
SR_ID_02.Age_01 = SR_ID_02.Age_01 - SR_ID_02.Mod; SR_ID_02.Mod = [];

PM_opt_00 = table(SR_opt_PM_00); PM_alt_A_00 = table(SR_alt_PM_A_00);
PM_alt_DA_00 = table(SR_alt_PM_DA_00); PM_alt_HB_00 = table(SR_alt_PM_HB_00);
PM_opt_00.SR_opt_PM_00 = double(string(PM_opt_00.SR_opt_PM_00));
PM_alt_A_00.SR_alt_PM_A_00 = double(string(PM_alt_A_00.SR_alt_PM_A_00));
PM_alt_DA_00.SR_alt_PM_DA_00 = double(string(PM_alt_DA_00.SR_alt_PM_DA_00));
PM_alt_HB_00.SR_alt_PM_HB_00 = double(string(PM_alt_HB_00.SR_alt_PM_HB_00));
PM_opt_01 = [SR_ID_02 PM_opt_00]; PM_alt_A_01 = [SR_ID_02 PM_alt_A_00];
PM_alt_DA_01 = [SR_ID_02 PM_alt_DA_00]; PM_alt_HB_01 = [SR_ID_02 PM_alt_HB_00];

PM_opt_02 = sortrows(PM_opt_01, [1:3, 6], 'ascend');
PM_alt_A_02 = sortrows(PM_alt_A_01, [1:3, 6], 'ascend');
PM_alt_DA_02 = sortrows(PM_alt_DA_01, [1:3, 6], 'ascend');
PM_alt_HB_02 = sortrows(PM_alt_HB_01, [1:3, 6], 'ascend');

PM_opt = table(PM_opt_02.SR_opt_PM_00, 'VariableNames', ...
    {'SR_opt_PM_00'}); 
PM_alt_A = table(PM_alt_A_02.SR_alt_PM_A_00, 'VariableNames', ...
    {'SR_alt_PM_A_00'}); 
PM_alt_DA = table(PM_alt_DA_02.SR_alt_PM_DA_00, 'VariableNames', ...
    {'SR_alt_PM_DA_00'});
PM_alt_HB = table(PM_alt_HB_02.SR_alt_PM_HB_00, 'VariableNames', ...
    {'SR_alt_PM_HB_00'});

SR_opt_EL_02 = [SR_opt_EL_02 Full_Loss_opt PM_opt];
SR_alt_EL_A_02 = [SR_alt_EL_A_02 Full_Loss_alt_A PM_alt_A];
SR_alt_EL_DA_02 = [SR_alt_EL_DA_02 Full_Loss_alt_DA PM_alt_DA];
SR_alt_EL_HB_02 = [SR_alt_EL_HB_02 Full_Loss_alt_HB PM_alt_HB];
% https://www.mathworks.com/matlabcentral/answers/179129-is-it-possible-to-rename-columns-in-a-table
SR_opt_EL_02.Properties.VariableNames([12 13]) = {'Full_Loss' 'Pre_Mort'};
SR_alt_EL_A_02.Properties.VariableNames([12 13]) = ...
    {'Full_Loss' 'Pre_Mort'};
SR_alt_EL_DA_02.Properties.VariableNames([12 13]) = ...
    {'Full_Loss' 'Pre_Mort'};
SR_alt_EL_HB_02.Properties.VariableNames([12 13]) = ...
    {'Full_Loss' 'Pre_Mort'};

% https://www.mathworks.com/help/matlab/matlab_prog/create-a-table.html
SR_opt_EL_02.EL_Dec = SR_opt_EL_02.Pre_Mort.*SR_opt_EL_02.EL_Dec;
SR_alt_EL_A_02.EL_Dec = SR_alt_EL_A_02.Pre_Mort.*SR_alt_EL_A_02.EL_Dec;
SR_alt_EL_DA_02.EL_Dec = SR_alt_EL_DA_02.Pre_Mort.*SR_alt_EL_DA_02.EL_Dec;
SR_alt_EL_HB_02.EL_Dec = SR_alt_EL_HB_02.Pre_Mort.*SR_alt_EL_HB_02.EL_Dec;

SR_opt_EL_02.EL_Tot = SR_opt_EL_02.Pre_Mort.*SR_opt_EL_02.EL_Dec + ...
    (1-SR_opt_EL_02.Pre_Mort).*SR_opt_EL_02.EL_Tot;
SR_alt_EL_A_02.EL_Tot = SR_alt_EL_A_02.Pre_Mort.*...
    SR_alt_EL_A_02.EL_Dec + (1-SR_alt_EL_A_02.Pre_Mort).*...
    SR_alt_EL_A_02.EL_Tot;
SR_alt_EL_DA_02.EL_Tot = SR_alt_EL_DA_02.Pre_Mort.*...
    SR_alt_EL_DA_02.EL_Dec + (1-SR_alt_EL_DA_02.Pre_Mort).*...
    SR_alt_EL_DA_02.EL_Tot;
SR_alt_EL_HB_02.EL_Tot = SR_alt_EL_HB_02.Pre_Mort.*...
    SR_alt_EL_HB_02.EL_Dec + (1-SR_alt_EL_HB_02.Pre_Mort).*...
    SR_alt_EL_HB_02.EL_Tot;

SR_opt_EL_02 = sortrows(SR_opt_EL_02, 'Index', 'ascend');
SR_alt_EL_A_02 = sortrows(SR_alt_EL_A_02, 'Index', 'ascend');
SR_alt_EL_DA_02 = sortrows(SR_alt_EL_DA_02, 'Index', 'ascend');
SR_alt_EL_HB_02 = sortrows(SR_alt_EL_HB_02, 'Index', 'ascend');

SR_00 = ...
   {SR_opt_EL_02.EL_Abs; SR_alt_EL_A_02.EL_Abs; SR_alt_EL_DA_02.EL_Abs; ...
    SR_alt_EL_HB_02.EL_Abs; ...
    SR_opt_EL_02.EL_Pre; SR_alt_EL_A_02.EL_Pre; SR_alt_EL_DA_02.EL_Pre; ...
    SR_alt_EL_HB_02.EL_Pre; ...
    SR_opt_EL_02.EL_Dec; SR_alt_EL_A_02.EL_Dec; SR_alt_EL_DA_02.EL_Dec; ...
    SR_alt_EL_HB_02.EL_Dec; ...
    SR_opt_EL_02.EL_Tot; SR_alt_EL_A_02.EL_Tot; SR_alt_EL_DA_02.EL_Tot; ...
    SR_alt_EL_HB_02.EL_Tot; ...
    SR_opt_EL_02.HS;     SR_alt_EL_A_02.HS;     SR_alt_EL_DA_02.HS; ...
    SR_alt_EL_HB_02.HS; ...
    SR_opt_EL_02.Age_01; SR_alt_EL_A_02.Age_01; SR_alt_EL_DA_02.Age_01; ...
    SR_alt_EL_HB_02.Age_01};

SR_ID = array2table(SR_ID);
Labels_Chars = {'Gender', 'Ethnicity', 'Education'};
SR_ID.Properties.VariableNames = Labels_Chars;

for i = 1:size(SR_Opt_Cohort(:, 4:size(SR_Opt_Cohort, 2)),2)
    
%   NOTE 093018: You cannot concatenate a string and a number; therefore,
%   to perform the operation below you must first convert each integer into
%   a string.
    Labels_00(i) = strcat("Time_", num2str(i-1));
    
end

% https://www.mathworks.com/help/matlab/ref/cellstr.html
Labels_00 = cellstr([Labels_Chars Labels_00]);

% https://stackoverflow.com/questions/27642838/reshape-matlab-table/27643127#27643127
for i = 1:size(SR_00, 1)
    
    SR_00{i} = array2table(reshape(SR_00{i}, ...
            [size(SR_Opt_Cohort(:, 4:size(SR_Opt_Cohort, 2)), 1), ...
            size(SR_Opt_Cohort(:, 4:size(SR_Opt_Cohort, 2)), 2)]));

    SR_00{i} = [SR_ID SR_00{i}];
    SR_00{i}.Properties.VariableNames = Labels_00;
    
end

%% Estimating T2DM Prevalence

Health_States = {'Healthy', 'T2DM', 'Dead'};
% Cell arrays can hold anything, evidently - including functions. This
% allows us to apply a slightly different function to each cell of a second
% an array with the same dimensions as the first.
Functions = cell(3, 1);
% The function STRCMP compares two strings, returning 1 if they are the
% same and 0 otherwise.
% https://www.mathworks.com/help/matlab/ref/strcmp.html
for i = 1:3
    % This function essentially sums the number of times we count a
    % particular health state in our simulation results. This is done by
    % comparing each observed occurrence of said state against a state of
    % interest found at some position within the cell array 
    % "Health_States", with each match contributing one to the total count.
    Functions{i} = @(x) sum(strcmp(x, Health_States{i}));
end

SR_Counts_opt = SR_00{17}; SR_Counts_alt_A = SR_00{18};
SR_Counts_alt_DA = SR_00{19}; SR_Counts_alt_HB = SR_00{20};
SR_HS_Counts_00 = cell(12, 1);
% https://www.mathworks.com/help/matlab/matlab_prog/calculations-on-tables.html
% https://www.mathworks.com/help/matlab/ref/varfun.html

% Every fourth cell of Sim_HS_Counts_00, beginning with the first, contains 
% data on the number of healthy, diabetic, and deceased individuals 
% observed at each point in time, when the optimal policy is applied.
% The remaining cells hold the same information for the three alternative
% policies under consideration.
% 
% Remember that we have four functions contained in a cell array. How do
% we determine which function should apply to each cell in Sim_HS_Counts?
% The find function, as applied here, tells us where i appears in the
% sequence [1, 5, 9], as denoted by an integer in [1, 4]. This indicates
% which function we are applying, and, by extension, which health state we
% are counting. As such, we are telling the program to count the number of
% healthy individuals when the optimal policy is applied, then when the
% first alternate policy is applied, etc. Then we do the same thing for
% individuals in the diabetic and deceased states.

SimCounts = cell(3, 1);
for i = [1, 5, 9]
    
    SR_HS_Counts_00{i} = varfun(Functions{find([1, 5, 9] == i)}, ...
        SR_00{17}(:, 4:size(SR_00{17}, 2)));
    SR_HS_Counts_00{i + 1} = varfun(Functions{find([1, 5, 9] == i)}, ...
        SR_00{18}(:, 4:size(SR_00{18}, 2)));
    SR_HS_Counts_00{i + 2} = varfun(Functions{find([1, 5, 9] == i)}, ...
        SR_00{19}(:, 4:size(SR_00{19}, 2)));
    SR_HS_Counts_00{i + 3} = varfun(Functions{find([1, 5, 9] == i)}, ...
        SR_00{20}(:, 4:size(SR_00{20}, 2)));
    SimCounts{find([1, 5, 9] == i), 1} = [SR_HS_Counts_00{i}; ...
       SR_HS_Counts_00{i + 1}; SR_HS_Counts_00{i + 2}; ...
       SR_HS_Counts_00{i + 3}];

end

func_01 = @(x) transpose(100*(x{:,:}/size(SR_Opt_Cohort, 1)));
SimCounts = cellfun(func_01, SimCounts, 'UniformOutput', false);
SimCounts = array2table([SimCounts{1} SimCounts{2} SimCounts{3}]);
SimCounts = SimCounts(:, [1 5 9 2 6 10 3 7 11 4 8 12]);

SimCounts_Labels = {'Opt_Healthy', 'Opt_T2DM', 'Opt_Dead', ...
    'Alt_A_Healthy', 'Alt_A_T2DM', 'Alt_A_Dead', ...
    'Alt_DA_Healthy', 'Alt_DA_T2DM', 'Alt_DA_Dead', ...
    'Alt_HB_Healthy', 'Alt_HB_T2DM', 'Alt_HB_Dead'};
SimCounts.Properties.VariableNames = SimCounts_Labels;
SimCounts.Opt_NHealthy = SimCounts.Opt_T2DM + SimCounts.Opt_Dead;
SimCounts.Alt_A_NHealthy = SimCounts.Alt_A_T2DM + SimCounts.Alt_A_Dead;
SimCounts.Alt_DA_NHealthy = SimCounts.Alt_DA_T2DM + SimCounts.Alt_DA_Dead;
SimCounts.Alt_HB_NHealthy = SimCounts.Alt_HB_T2DM + SimCounts.Alt_HB_Dead;
SimCounts = SimCounts(:, [1:3 13 4:6 14 7:9 15 10:12 16]);

SR_Pre = cell(4, 1);
SR_Pre{1, 1} = 100*(SimCounts.Opt_T2DM(2:end)./...
    (SimCounts.Opt_Healthy(2:end) + SimCounts.Opt_T2DM(2:end)));
SR_Pre{2, 1} = 100*(SimCounts.Alt_A_T2DM(2:end)./...
    (SimCounts.Alt_A_Healthy(2:end) + SimCounts.Alt_A_T2DM(2:end)));
SR_Pre{3, 1} = 100*(SimCounts.Alt_DA_T2DM(2:end)./...
    (SimCounts.Alt_DA_Healthy(2:end) + SimCounts.Alt_DA_T2DM(2:end)));
SR_Pre{4, 1} = 100*(SimCounts.Alt_HB_T2DM(2:end)./...
    (SimCounts.Alt_HB_Healthy(2:end) + SimCounts.Alt_HB_T2DM(2:end)));

TableNames_Pre = {'Elapsed_Time', 'Prev_Opt', 'Prev_Alt_A', ...
    'Prev_Alt_DA', 'Prev_Alt_HB'};
SR_Pre = array2table([(1:T)' SR_Pre{1} SR_Pre{2} SR_Pre{3} SR_Pre{4}]);
SR_Pre.Properties.VariableNames = TableNames_Pre;
SR_Pre{:,:}(isnan(SR_Pre{:,:})) = 0;

%% Estimating T2DM Incidence

SR_Inc_opt_00 = SR_00{17}; SR_Inc_alt_A_00 = SR_00{18};
SR_Inc_alt_DA_00 = SR_00{19}; SR_Inc_alt_HB_00 = SR_00{20};

SR_Inc_opt_01 = SR_Inc_opt_00{:, 4:(end-1)};
SR_Inc_opt_02 = SR_Inc_opt_00{:, 5:end};
SR_Inc_alt_A_01 = SR_Inc_alt_A_00{:, 4:(end-1)};
SR_Inc_alt_A_02 = SR_Inc_alt_A_00{:, 5:end};
SR_Inc_alt_DA_01 = SR_Inc_alt_DA_00{:, 4:(end-1)};
SR_Inc_alt_DA_02 = SR_Inc_alt_DA_00{:, 5:end};
SR_Inc_alt_HB_01 = SR_Inc_alt_HB_00{:, 4:(end-1)};
SR_Inc_alt_HB_02 = SR_Inc_alt_HB_00{:, 5:end};

% The cases of interest are those in which a person is healthy in period t 
% and is diagnosed with T2DM in period t + 1.
SR_Inc_opt_03 = 1.*(SR_Inc_opt_01 == "Healthy" & SR_Inc_opt_02 == "T2DM");
SR_Inc_alt_A_03 = 1.*(SR_Inc_alt_A_01 == "Healthy" & ...
    SR_Inc_alt_A_02 == "T2DM");
SR_Inc_alt_DA_03 = 1.*(SR_Inc_alt_DA_01 == "Healthy" & ...
    SR_Inc_alt_DA_02 == "T2DM");
SR_Inc_alt_HB_03 = 1.*(SR_Inc_alt_HB_01 == "Healthy" & ...
    SR_Inc_alt_HB_02 == "T2DM");

SR_Inc_opt_04 = sum(SR_Inc_opt_03, 1);
SR_Inc_alt_A_04 = sum(SR_Inc_alt_A_03, 1);
SR_Inc_alt_DA_04 = sum(SR_Inc_alt_DA_03, 1);
SR_Inc_alt_HB_04 = sum(SR_Inc_alt_HB_03, 1);
% This creates a table indicating the population at risk (i.e., who are 
% healthy) in the preceding period.
SR_Inc_opt_05 = table2array(SR_HS_Counts_00{1}(:, 2:end)); 
SR_Inc_alt_A_05 = table2array(SR_HS_Counts_00{2}(:, 2:end));
SR_Inc_alt_DA_05 = table2array(SR_HS_Counts_00{3}(:, 2:end));
SR_Inc_alt_HB_05 = table2array(SR_HS_Counts_00{4}(:, 2:end));

% Incidence is expressed in terms of the number of individuals with T2DM at
% time t + 1 who were at risk for - but did not have - T2DM in the 
% preceding period (i.e., period t).
SR_Inc_opt_06 = (SR_Inc_opt_04./SR_Inc_opt_05)*1000;
SR_Inc_alt_A_06 = (SR_Inc_alt_A_04./SR_Inc_alt_A_05)*1000;
SR_Inc_alt_DA_06 = (SR_Inc_alt_DA_04./SR_Inc_alt_DA_05)*1000;
SR_Inc_alt_HB_06 = (SR_Inc_alt_HB_04./SR_Inc_alt_HB_05)*1000;

SimTime_Inc = transpose(1:T);
TableNames_Inc = {'Elapsed_Time', 'Inc_Opt', 'Inc_Alt_A', ...
    'Inc_Alt_DA', 'Inc_Alt_HB'};
SR_Inc = array2table([SimTime_Inc transpose(SR_Inc_opt_06) ...
    transpose(SR_Inc_alt_A_06) transpose(SR_Inc_alt_DA_06) ...
    transpose(SR_Inc_alt_HB_06)]);
SR_Inc.Properties.VariableNames = TableNames_Inc;

%% Plotting Prevalence and Incidence of T2DM in the Cohort - Alt #1

figure(1);
subplot(2, 2, 1);
Prev_opt_HS = bar((1:T)', SimCounts{2:end, 1:3}, 'stacked');
title('Simulated Health States - Optimal Policy');
xlabel('Time Elapsed (years)');
ylabel('Percent of Cohort in Each State (%)');

% https://blogs.mathworks.com/loren/2015/09/16/individual-axis-customization/
ax = gca; % current axes
ax.YAxis.Exponent = 0;
ax.YAxis.TickLabelFormat = '%.2f%%';
ax.YLim = [0, 100];
legend('Healthy', 'T2DM', 'Dead')

subplot(2, 2, 2);
Prev_alt_A_HS = bar((1:T)', SimCounts{2:end, 5:7}, 'stacked');
title('Simulated Health States - Alternate Policy #1');
xlabel('Time Elapsed (years)');
ylabel('Percent of Cohort in Each State (%)');

ax = gca; % current axes
ax.YAxis.Exponent = 0;
ax.YAxis.TickLabelFormat = '%.2f%%';
ax.YLim = [0, 100];
legend('Healthy', 'T2DM', 'Dead')

subplot(2, 2, 3);
% https://plot.ly/matlab/multiple-axes/
title('T2DM Prevalence & Incidence - Optimal Policy');
xlabel('Time Elapsed (years)');
ax = gca;

yyaxis left
% https://www.mathworks.com/matlabcentral/answers/39221-how-do-i-get-the-color-gray
% https://www.mathworks.com/matlabcentral/answers/43949-how-to-change-the-colors-of-bars-in-a-bar-graph
Pre_opt_HS = bar((1:T)', SR_Pre.Prev_Opt, ...
    'FaceColor', [0.85 0.85 0.85], 'EdgeColor', [0.85 0.85 0.85]);
ylabel('Disease Prevalence (%)');
ax.YAxis(1).TickLabelFormat = '%.2f%%';
ax.YAxis(1).Color = 'k';
ylim(ax, [0, 100])

yyaxis right
Inc_opt_HS = plot((1:T)', SR_Inc.Inc_Opt, ...
    'LineWidth', 2, 'Color', [0.35 0.35 0.35]);
ylabel('Incidence Rates per 1,000 Manitobans');
ax.YAxis(2).Exponent = 3;
ax.YAxis(2).Color = 'k';
ylim(ax, [0, max(not(or(isnan(SR_Inc.Inc_Alt_A), ...
    isinf(SR_Inc.Inc_Alt_A))))*1000])
grid on; ax.YMinorGrid = 'on';
legend('Prevalence', 'Incidence')

subplot(2, 2, 4);
title('T2DM Prevalence & Incidence - Alternate Policy #1');
xlabel('Time Elapsed (years)');
ax = gca;

yyaxis left
% https://www.mathworks.com/matlabcentral/answers/39221-how-do-i-get-the-color-gray
% https://www.mathworks.com/matlabcentral/answers/43949-how-to-change-the-colors-of-bars-in-a-bar-graph
Pre_alt_HS = bar((1:T)', SR_Pre.Prev_Alt_A, ...
    'FaceColor', [0.85 0.85 0.85], 'EdgeColor', [0.85 0.85 0.85]);
ylabel('Disease Prevalence (%)');
ax.YAxis(1).TickLabelFormat = '%.2f%%';
ax.YAxis(1).Color = 'k';
ylim(ax, [0, 100])

yyaxis right
Inc_alt_HS = plot((1:T)', SR_Inc.Inc_Alt_A, ...
    'LineWidth', 2, 'Color', [0.35 0.35 0.35]);
ylabel('Incidence Rates per 1,000 Manitobans');
ax.YAxis(2).Exponent = 3;
ax.YAxis(2).Color = 'k';
ylim(ax, [0, max(not(or(isnan(SR_Inc.Inc_Alt_A), ...
    isinf(SR_Inc.Inc_Alt_A))))*1000])
grid on; ax.YMinorGrid = 'on';
legend('Prevalence', 'Incidence')

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
% https://www.mathworks.com/help/matlab/ref/print.html
print('-f1', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig1"), '-djpeg', '-r300');

%% Plotting Prevalence and Incidence of T2DM in the Cohort - Alt #2

figure(2);
subplot(2, 2, 1);
Prev_opt_HS = bar((1:T)', SimCounts{2:end, 1:3}, 'stacked');
title('Simulated Health States - Optimal Policy');
xlabel('Time Elapsed (years)');
ylabel('Percent of Cohort in Each State (%)');

% https://blogs.mathworks.com/loren/2015/09/16/individual-axis-customization/
ax = gca; % current axes
ax.YAxis.Exponent = 0;
ax.YAxis.TickLabelFormat = '%.2f%%';
ax.YLim = [0, 100];
legend('Healthy', 'T2DM', 'Dead')

subplot(2, 2, 2);
Prev_alt_A_HS = bar((1:T)', SimCounts{2:end, 9:11}, 'stacked');
title('Simulated Health States - Alternate Policy #2');
xlabel('Time Elapsed (years)');
ylabel('Percent of Cohort in Each State (%)');

ax = gca; % current axes
ax.YAxis.Exponent = 0;
ax.YAxis.TickLabelFormat = '%.2f%%';
ax.YLim = [0, 100];
legend('Healthy', 'T2DM', 'Dead')

subplot(2, 2, 3);
% https://plot.ly/matlab/multiple-axes/
title('T2DM Prevalence & Incidence - Optimal Policy');
xlabel('Time Elapsed (years)');
ax = gca;

yyaxis left
% https://www.mathworks.com/matlabcentral/answers/39221-how-do-i-get-the-color-gray
% https://www.mathworks.com/matlabcentral/answers/43949-how-to-change-the-colors-of-bars-in-a-bar-graph
Pre_opt_HS = bar((1:T)', SR_Pre.Prev_Opt, ...
    'FaceColor', [0.85 0.85 0.85], 'EdgeColor', [0.85 0.85 0.85]);
ylabel('Disease Prevalence (%)');
ax.YAxis(1).TickLabelFormat = '%.2f%%';
ax.YAxis(1).Color = 'k';
ylim(ax, [0, 100])

yyaxis right
Inc_opt_HS = plot((1:T)', SR_Inc.Inc_Opt, ...
    'LineWidth', 2, 'Color', [0.35 0.35 0.35]);
ylabel('Incidence Rates per 1,000 Manitobans');
ax.YAxis(2).Exponent = 3;
ax.YAxis(2).Color = 'k';
ylim(ax, [0, max(not(or(isnan(SR_Inc.Inc_Alt_DA), ...
    isinf(SR_Inc.Inc_Alt_DA))))*1000])
grid on; ax.YMinorGrid = 'on';
legend('Prevalence', 'Incidence')

subplot(2, 2, 4);
title('T2DM Prevalence & Incidence - Alternate Policy #2');
xlabel('Time Elapsed (years)');
ax = gca;

yyaxis left
% https://www.mathworks.com/matlabcentral/answers/39221-how-do-i-get-the-color-gray
% https://www.mathworks.com/matlabcentral/answers/43949-how-to-change-the-colors-of-bars-in-a-bar-graph
Pre_alt_HS = bar((1:T)', SR_Pre.Prev_Alt_DA, ...
    'FaceColor', [0.85 0.85 0.85], 'EdgeColor', [0.85 0.85 0.85]);
ylabel('Disease Prevalence (%)');
ax.YAxis(1).TickLabelFormat = '%.2f%%';
ax.YAxis(1).Color = 'k';
ylim(ax, [0, 100])

yyaxis right
Inc_alt_HS = plot((1:T)', SR_Inc.Inc_Alt_DA, ...
    'LineWidth', 2, 'Color', [0.35 0.35 0.35]);
ylabel('Incidence Rates per 1,000 Manitobans');
ax.YAxis(2).Exponent = 3;
ax.YAxis(2).Color = 'k';
ylim(ax, [0, max(not(or(isnan(SR_Inc.Inc_Alt_DA), ...
    isinf(SR_Inc.Inc_Alt_DA))))*1000])
grid on; ax.YMinorGrid = 'on';
legend('Prevalence', 'Incidence')

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
print('-f2', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig2"), '-djpeg', '-r300');

%% Plotting Prevalence and Incidence of T2DM in the Cohort - Alt #3

% 011219: Here we plot the new scenario. We would expect prevalence and
% incidence to be more or less the same as in Alt #2, since individuals in
% both scenarios opt not to adhere with LTPA or healthy eating habits when
% they are healthy - the scenarios only begin to diverge following
% diagnosis.
figure(3);
subplot(2, 2, 1);
Prev_opt_HS = bar((1:T)', SimCounts{2:end, 1:3}, 'stacked');
title('Simulated Health States - Optimal Policy');
xlabel('Time Elapsed (years)');
ylabel('Percent of Cohort in Each State (%)');

ax = gca;
ax.YAxis.Exponent = 0;
ax.YAxis.TickLabelFormat = '%.2f%%';
ax.YLim = [0, 100];
legend('Healthy', 'T2DM', 'Dead')

subplot(2, 2, 2);
Prev_alt_A_HS = bar((1:T)', SimCounts{2:end, 13:15}, 'stacked');
title('Simulated Health States - Alternate Policy #3');
xlabel('Time Elapsed (years)');
ylabel('Percent of Cohort in Each State (%)');

ax = gca;
ax.YAxis.Exponent = 0;
ax.YAxis.TickLabelFormat = '%.2f%%';
ax.YLim = [0, 100];
legend('Healthy', 'T2DM', 'Dead')

subplot(2, 2, 3);
title('T2DM Prevalence & Incidence - Optimal Policy');
xlabel('Time Elapsed (years)');
ax = gca;

yyaxis left
Pre_opt_HS = bar((1:T)', SR_Pre.Prev_Opt, ...
    'FaceColor', [0.85 0.85 0.85], 'EdgeColor', [0.85 0.85 0.85]);
ylabel('Disease Prevalence (%)');
ax.YAxis(1).TickLabelFormat = '%.2f%%';
ax.YAxis(1).Color = 'k';
ylim(ax, [0, 100])

yyaxis right
Inc_opt_HS = plot((1:T)', SR_Inc.Inc_Opt, ...
    'LineWidth', 2, 'Color', [0.35 0.35 0.35]);
ylabel('Incidence Rates per 1,000 Manitobans');
ax.YAxis(2).Exponent = 3;
ax.YAxis(2).Color = 'k';
ylim(ax, [0, max(not(or(isnan(SR_Inc.Inc_Alt_HB), ...
    isinf(SR_Inc.Inc_Alt_HB))))*1000])
grid on; ax.YMinorGrid = 'on';
legend('Prevalence', 'Incidence')

subplot(2, 2, 4);
title('T2DM Prevalence & Incidence - Alternate Policy #3');
xlabel('Time Elapsed (years)');
ax = gca;

yyaxis left
Pre_alt_HS = bar((1:T)', SR_Pre.Prev_Alt_HB, ...
    'FaceColor', [0.85 0.85 0.85], 'EdgeColor', [0.85 0.85 0.85]);
ylabel('Disease Prevalence (%)');
ax.YAxis(1).TickLabelFormat = '%.2f%%';
ax.YAxis(1).Color = 'k';
ylim(ax, [0, 100])

yyaxis right
Inc_alt_HS = plot((1:T)', SR_Inc.Inc_Alt_HB, ...
    'LineWidth', 2, 'Color', [0.35 0.35 0.35]);
ylabel('Incidence Rates per 1,000 Manitobans');
ax.YAxis(2).Exponent = 3;
ax.YAxis(2).Color = 'k';
ylim(ax, [0, max(not(or(isnan(SR_Inc.Inc_Alt_HB), ...
    isinf(SR_Inc.Inc_Alt_HB))))*1000])
grid on; ax.YMinorGrid = 'on';
legend('Prevalence', 'Incidence')

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
print('-f3', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig3"), '-djpeg', '-r300');

%% Estimating Efficiency Losses - Alt #1

Size_SimEL = size(SR_00{1}, 2);
SR_EL_abs_opt = transpose(table2array( ...
    varfun(@sum, SR_00{1}(:, 5:Size_SimEL))));
SR_EL_pre_opt = transpose(table2array( ...
    varfun(@sum, SR_00{5}(:, 5:Size_SimEL))));
SR_EL_dec_opt = transpose(table2array( ...
    varfun(@sum, SR_00{9}(:, 5:Size_SimEL))));
SR_EL_tot_opt = transpose(table2array( ...
    varfun(@sum, SR_00{13}(:, 5:Size_SimEL))));
SR_EL_opt = array2table([(1:T)' SR_EL_abs_opt SR_EL_pre_opt ...
    SR_EL_dec_opt SR_EL_tot_opt]);
TableNames = {'Elapsed_Time', 'EL_Absenteeism', 'EL_Presenteeism', ...
    'EL_Mortality', 'EL_Total'};
SR_EL_opt.Properties.VariableNames = TableNames;

SR_EL_abs_alt_A = transpose(table2array( ...
    varfun(@sum, SR_00{2}(:, 5:Size_SimEL))));
SR_EL_pre_alt_A = transpose(table2array( ...
    varfun(@sum, SR_00{6}(:, 5:Size_SimEL))));
SR_EL_dec_alt_A = transpose(table2array( ...
    varfun(@sum, SR_00{10}(:, 5:Size_SimEL))));
SR_EL_tot_alt_A = transpose(table2array( ...
    varfun(@sum, SR_00{14}(:, 5:Size_SimEL))));
SR_EL_alt_A = array2table([(1:T)' SR_EL_abs_alt_A SR_EL_pre_alt_A ...
    SR_EL_dec_alt_A SR_EL_tot_alt_A]);
SR_EL_alt_A.Properties.VariableNames = TableNames;

% 102718: Calculate the difference in productivity losses between the two
% scenarios.
SEL_Opt = array2table(SR_EL_opt{:,{'Elapsed_Time', 'EL_Total'}}); 
SEL_Alt_A = array2table(SR_EL_alt_A{:,{'Elapsed_Time', 'EL_Total'}});
TableNames_Opt = {'Elapsed_Time', 'EL_Total_Opt'};
TableNames_Alt_A = {'Elapsed_Time', 'EL_Total_Alt_A'};
SEL_Opt.Properties.VariableNames = TableNames_Opt;
SEL_Alt_A.Properties.VariableNames = TableNames_Alt_A;
SEL_Comp_1 = join(SEL_Opt, SEL_Alt_A);

SEL_Comp_1.EL_Diff = SEL_Comp_1.EL_Total_Alt_A - SEL_Comp_1.EL_Total_Opt;
Disc_Rate = 0.03;
SEL_Comp_1.DiscFact = 1./power((1+Disc_Rate), SEL_Comp_1.Elapsed_Time);
SEL_Comp_1{:,:}(isnan(SEL_Comp_1{:,:})) = 0;
SEL_Comp_1.EL_Diff_Disc = SEL_Comp_1.EL_Diff.*SEL_Comp_1.DiscFact;
SEL_Diff_Total_1 = sum(SEL_Comp_1.EL_Diff(1:70));
SEL_Diff_PV_Total_1 = sum(SEL_Comp_1.EL_Diff_Disc(1:70));

%% Visualizing Efficiency Losses - Alt #1

% https://www.mathworks.com/matlabcentral/answers/103188-how-do-i-format-tick-labels-before-r2016b
% https://www.mathworks.com/help/matlab/ref/gca.html
figure(4); 
subplot(2, 2, 1);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Absenteeism, 's--r');

hold on;

plot(SR_EL_alt_A.Elapsed_Time, SR_EL_alt_A.EL_Absenteeism, 'o:b');
title('Simulated Efficiency Losses: Absenteeism - Alt #1');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

% https://blogs.mathworks.com/loren/2015/09/16/individual-axis-customization/#2e1dff59-8864-497d-8e0a-84549d8004c9
ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
% https://www.mathworks.com/help/matlab/ref/matlab.graphics.axis.axes-properties.html
ax.YLim = [0, max(SR_EL_alt_A.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

subplot(2, 2, 2);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Presenteeism, 's--r');

hold on;

plot(SR_EL_alt_A.Elapsed_Time, SR_EL_alt_A.EL_Presenteeism, 'o:b');
title('Simulated Efficiency Losses: Presenteeism - Alt #1');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_A.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

subplot(2, 2, 3);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Mortality, 's--r');

hold on;

plot(SR_EL_alt_A.Elapsed_Time, SR_EL_alt_A.EL_Mortality, 'o:b');
title('Simulated Efficiency Losses: Mortality - Alt #1');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_A.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

subplot(2, 2, 4);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Total, 's--r');

hold on;

plot(SR_EL_alt_A.Elapsed_Time, SR_EL_alt_A.EL_Total, 'o:b');
title('Simulated Efficiency Losses: Total - Alt #1');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_A.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
print('-f4', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig4"), '-djpeg', '-r300');

%% Estimating Efficiency Losses - Alt #2

SR_EL_abs_alt_DA = transpose(table2array( ...
    varfun(@sum, SR_00{3}(:, 5:Size_SimEL))));
SR_EL_pre_alt_DA = transpose(table2array( ...
    varfun(@sum, SR_00{7}(:, 5:Size_SimEL))));
SR_EL_dec_alt_DA = transpose(table2array( ...
    varfun(@sum, SR_00{11}(:, 5:Size_SimEL))));
SR_EL_tot_alt_DA = transpose(table2array( ...
    varfun(@sum, SR_00{15}(:, 5:Size_SimEL))));
SR_EL_alt_DA = array2table([(1:T)' SR_EL_abs_alt_DA SR_EL_pre_alt_DA ...
    SR_EL_dec_alt_DA SR_EL_tot_alt_DA]);
SR_EL_alt_DA.Properties.VariableNames = TableNames;

SEL_Alt_DA = array2table(SR_EL_alt_DA{:,{'Elapsed_Time', 'EL_Total'}});
TableNames_Alt_DA = {'Elapsed_Time', 'EL_Total_Alt_DA'};
SEL_Alt_DA.Properties.VariableNames = TableNames_Alt_DA;
SEL_Comp_2 = join(SEL_Opt, SEL_Alt_DA);

SEL_Comp_2.EL_Diff = SEL_Comp_2.EL_Total_Alt_DA - SEL_Comp_2.EL_Total_Opt;
SEL_Comp_2.DiscFact = 1./power((1+Disc_Rate), SEL_Comp_2.Elapsed_Time);
SEL_Comp_2{:,:}(isnan(SEL_Comp_2{:,:})) = 0;
SEL_Comp_2.EL_Diff_Disc = SEL_Comp_2.EL_Diff.*SEL_Comp_2.DiscFact;
SEL_Diff_Total_2 = sum(SEL_Comp_2.EL_Diff(1:70));
SEL_Diff_PV_Total_2 = sum(SEL_Comp_2.EL_Diff_Disc(1:70));

%% Visualizing Efficiency Losses - Alt #2

figure(5); 
subplot(2, 2, 1);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Absenteeism, 's--r');

hold on;

plot(SR_EL_alt_DA.Elapsed_Time, SR_EL_alt_DA.EL_Absenteeism, 'o:b');
title('Simulated Efficiency Losses: Absenteeism - Alt #2');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

subplot(2, 2, 2);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Presenteeism, 's--r');

hold on;

plot(SR_EL_alt_DA.Elapsed_Time, SR_EL_alt_DA.EL_Presenteeism, 'o:b');
title('Simulated Efficiency Losses: Presenteeism - Alt #2');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

subplot(2, 2, 3);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Mortality, 's--r');

hold on;

plot(SR_EL_alt_DA.Elapsed_Time, SR_EL_alt_DA.EL_Mortality, 'o:b');
title('Simulated Efficiency Losses: Mortality - Alt #2');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

subplot(2, 2, 4);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Total, 's--r');

hold on;

plot(SR_EL_alt_DA.Elapsed_Time, SR_EL_alt_DA.EL_Total, 'o:b');
title('Simulated Efficiency Losses: Total - Alt #2');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate')

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
print('-f5', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig5"), '-djpeg', '-r300');

%% Estimating Efficiency Losses - Alt #3

SR_EL_abs_alt_HB = transpose(table2array( ...
    varfun(@sum, SR_00{4}(:, 5:Size_SimEL))));
SR_EL_pre_alt_HB = transpose(table2array( ...
    varfun(@sum, SR_00{8}(:, 5:Size_SimEL))));
SR_EL_dec_alt_HB = transpose(table2array( ...
    varfun(@sum, SR_00{12}(:, 5:Size_SimEL))));
SR_EL_tot_alt_HB = transpose(table2array( ...
    varfun(@sum, SR_00{16}(:, 5:Size_SimEL))));
SR_EL_alt_HB = array2table([(1:T)' SR_EL_abs_alt_HB SR_EL_pre_alt_HB ...
    SR_EL_dec_alt_HB SR_EL_tot_alt_HB]);
SR_EL_alt_HB.Properties.VariableNames = TableNames;

SEL_Alt_HB = array2table(SR_EL_alt_HB{:,{'Elapsed_Time', 'EL_Total'}});
TableNames_Alt_HB = {'Elapsed_Time', 'EL_Total_Alt_HB'};
SEL_Alt_HB.Properties.VariableNames = TableNames_Alt_HB;
SEL_Comp_3_Opt = join(SEL_Opt, SEL_Alt_HB);
SEL_Comp_3_Alt_DA = join(SEL_Alt_DA, SEL_Alt_HB);

SEL_Comp_3_Opt.EL_Diff = SEL_Comp_3_Opt.EL_Total_Alt_HB - ...
    SEL_Comp_3_Opt.EL_Total_Opt;
SEL_Comp_3_Opt.DiscFact = 1./power((1+Disc_Rate), ...
    SEL_Comp_3_Opt.Elapsed_Time);
SEL_Comp_3_Opt{:,:}(isnan(SEL_Comp_3_Opt{:,:})) = 0;
SEL_Comp_3_Opt.EL_Diff_Disc = SEL_Comp_3_Opt.EL_Diff.*...
    SEL_Comp_3_Opt.DiscFact;
SEL_Diff_Total_3_Opt = sum(SEL_Comp_3_Opt.EL_Diff(1:70));
SEL_Diff_PV_Total_3_Opt = sum(SEL_Comp_3_Opt.EL_Diff_Disc(1:70));

SEL_Comp_3_Alt_DA.EL_Diff = SEL_Comp_3_Alt_DA.EL_Total_Alt_HB - ...
    SEL_Comp_3_Alt_DA.EL_Total_Alt_DA;
SEL_Comp_3_Alt_DA.DiscFact = 1./power((1+Disc_Rate), ...
    SEL_Comp_3_Alt_DA.Elapsed_Time);
SEL_Comp_3_Alt_DA{:,:}(isnan(SEL_Comp_3_Alt_DA{:,:})) = 0;
SEL_Comp_3_Alt_DA.EL_Diff_Disc = SEL_Comp_3_Alt_DA.EL_Diff.*...
    SEL_Comp_3_Alt_DA.DiscFact;
SEL_Diff_Total_3_Alt_DA = sum(SEL_Comp_3_Alt_DA.EL_Diff(1:70));
SEL_Diff_PV_Total_3_Alt_DA = sum(SEL_Comp_3_Alt_DA.EL_Diff_Disc(1:70));

%% Visualizing Efficiency Losses - Alt #3

figure(6); 
subplot(2, 2, 1);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Absenteeism, 's--r');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_HB.EL_Absenteeism, 'o:b');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_DA.EL_Absenteeism, 'p-g');

title('Simulated Efficiency Losses: Absenteeism - Alt #3');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate (Hybrid)', 'Alternate (Non-Adher.)')

subplot(2, 2, 2);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Presenteeism, 's--r');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_HB.EL_Presenteeism, 'o:b');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_DA.EL_Presenteeism, 'p-g');

title('Simulated Efficiency Losses: Presenteeism - Alt #3');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate (Hybrid)', 'Alternate (Non-Adher.)')

subplot(2, 2, 3);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Mortality, 's--r');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_HB.EL_Mortality, 'o:b');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_DA.EL_Mortality, 'p-g');

title('Simulated Efficiency Losses: Mortality - Alt #3');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca;
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate (Hybrid)', 'Alternate (Non-Adher.)')

subplot(2, 2, 4);

plot(SR_EL_opt.Elapsed_Time, SR_EL_opt.EL_Total, 's--r');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_HB.EL_Total, 'o:b');

hold on;

plot(SR_EL_alt_HB.Elapsed_Time, SR_EL_alt_DA.EL_Total, 'p-g');

title('Simulated Efficiency Losses: Total - Alt #3');
xlabel('Time Elapsed (years)');
ylabel('Value (millions of dollars)');

ax = gca; % current axes
ax.YAxis.Exponent = 6;
ax.YAxis.TickLabelFormat = '\x0024 %,d M';
ax.YLim = [0, max(SR_EL_alt_DA.EL_Total)];
grid on; ax.YMinorGrid = 'on';

hold off;
legend('Optimal', 'Alternate (Hybrid)', 'Alternate (Non-Adher.)')

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
print('-f6', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig6"), '-djpeg', '-r300');

%% Sub-Group Analysis

SR_EL_SG = cell(4, 1);
SG_Labels = {'Gender', 'Ethnicity', 'Education'};
SG_Labels2 = {'Opt', 'Alt1' 'Alt2', 'Alt3'}';

for i = 13:16
    for j = 1:3
        SR_EL_SG{i - 12, j} = SR_00{i, 1};
    end
end

for x = 1:4
    for y = 1:3
        
    SR_EL_SG1{x, y} = ...
        varfun(@sum, SR_EL_SG{x, y}, 'GroupingVariables', SG_Labels{1, y});
    SR_EL_SG2{x, y} = SR_EL_SG1{x, y}(:, 1);
    SR_EL_SG3{x, y} = ...
        array2table(sum(table2array(SR_EL_SG1{x, y}(:, 5:74)), 2));
    SR_EL_SG3{x, y}.Properties.VariableNames{'Var1'} = ...
        SG_Labels2{x, 1};
    SR_EL_SG4{x, y} = varfun(@(z) 100*(z./sum(z)), SR_EL_SG3{x, y});
    
    end
end

% https://stackoverflow.com/questions/23334819/flattening-a-cell-array
SR_EL_SG_Gen = [SR_EL_SG2{1, 1} horzcat(SR_EL_SG4{:, 1})];
SR_EL_SG_Eth = [SR_EL_SG2{1, 2} horzcat(SR_EL_SG4{:, 2})];
SR_EL_SG_Edu = [SR_EL_SG2{1, 3} horzcat(SR_EL_SG4{:, 3})];

SR_EL_SG_Gen2 = table2array(SR_EL_SG_Gen)'; Gen = SR_EL_SG_Gen2(2:5, :);
SR_EL_SG_Eth2 = table2array(SR_EL_SG_Eth)'; Eth = SR_EL_SG_Eth2(2:5, :);
SR_EL_SG_Edu2 = table2array(SR_EL_SG_Edu)'; Edu = SR_EL_SG_Edu2(2:5, :);

figure(7);
subplot(3, 1, 1); bar(Gen,'stacked')
title('Total Indirect Costs by Gender');
set(gca,'xticklabel',SG_Labels2)
ylabel('Proportion of Costs (%)');
legend('Male', 'Female', 'Location', 'northwest')
% https://www.mathworks.com/matlabcentral/answers/328419-how-to-label-stacked-bar
for i=1:size(Gen,1)
    for j=1:size(Gen,2)
        if Gen(i,j)>0
        labels_stacked=num2str(Gen(i,j),'%0.1f%%');
        hText = text(i, sum(Gen(i,1:j),2), labels_stacked);
        set(hText, 'VerticalAlignment','top', 'HorizontalAlignment', ...
            'center','FontSize',10,'FontWeight','Bold','Color','w');
        end
    end
end

subplot(3, 1, 2); bar(Eth,'stacked')
title('Total Indirect Costs by Ethnicity');
set(gca,'xticklabel',SG_Labels2)
ylabel('Proportion of Costs (%)');
legend('Caucasian', 'Afro-Caribbean', 'Asian-Indian', ...
    'Location', 'northwest')
for i=1:size(Eth,1)
    for j=[1, 3]
        if Eth(i,j)>0
        labels_stacked=num2str(Eth(i,j),'%0.1f%%');
        hText = text(i, sum(Eth(i,1:j),2), labels_stacked);
        set(hText, 'VerticalAlignment','top', 'HorizontalAlignment', ...
            'center','FontSize',10,'FontWeight','Bold','Color','w');
        end
    end
end
for i=1:size(Eth,1)
    for j=2
        if Eth(i,j)>0
        labels_stacked=num2str(Eth(i,j),'%0.1f%%');
        hText = text(i, sum(Eth(i,1:j),2), labels_stacked);
        set(hText, 'VerticalAlignment','Middle', 'HorizontalAlignment', ...
            'center','FontSize',10,'FontWeight','Bold','Color','w');
        end
    end
end

subplot(3, 1, 3); bar(Edu,'stacked')
title('Total Indirect Costs by Level of Educational Attainment');
set(gca,'xticklabel',SG_Labels2)
ylabel('Proportion of Costs (%)');
legend('< HS', '   HS', '< Bach.', '+ Bach.', 'Location', 'northwest')
for i=1:size(Edu,1)
    for j=2:4
        if Edu(i,j)>0
        labels_stacked=num2str(Edu(i,j),'%0.1f%%');
        hText = text(i, sum(Edu(i,1:j),2), labels_stacked);
        set(hText, 'VerticalAlignment','top', 'HorizontalAlignment', ...
            'center','FontSize',10,'FontWeight','Bold','Color','w');
        end
    end
end
for i=1:size(Edu,1)
    for j=1
        if Edu(i,j)>0
        labels_stacked=num2str(Edu(i,j),'%0.1f%%');
        hText = text(i, sum(Edu(i,1:j),2), labels_stacked);
        set(hText, 'VerticalAlignment','Middle', 'HorizontalAlignment', ...
            'center','FontSize',10,'FontWeight','Bold','Color','w');
        end
    end
end

set(gcf, 'Units','centimeters', 'Position',[0 0 40 20])
print('-f7', strcat(MAT_Path, "Simulation_Results\Figures\", ...
    "Fig7"), '-djpeg', '-r300');

%% Prepare Detailed Data Tables

Table_1 = array2table([(1:T)' SR_Pre.Prev_Opt SR_Pre.Prev_Alt_A ...
    SR_Inc.Inc_Opt SR_Inc.Inc_Alt_A ...
    SR_EL_opt.EL_Absenteeism SR_EL_alt_A.EL_Absenteeism ...
    SR_EL_opt.EL_Presenteeism SR_EL_alt_A.EL_Presenteeism ...
    SR_EL_opt.EL_Mortality SR_EL_alt_A.EL_Mortality ...
    SR_EL_opt.EL_Total SR_EL_alt_A.EL_Total]);
Table_1_Names = {'Elapsed_Time', 'Prev_Opt', 'Prev_Alt_A', 'Inc_Opt', ...
    'Inc_Alt_A', 'Abs_Opt', 'Abs_Alt_A', 'Pre_Opt', 'Pre_Alt_A', ...
    'Dec_Opt', 'Dec_Alt_A', 'Tot_Opt', 'Tot_Alt_A'};
Table_1.Properties.VariableNames = Table_1_Names;
Table_1{:,2:3} = Table_1{:,2:3}/100;

Table_2 = array2table([(1:T)' SR_Pre.Prev_Opt SR_Pre.Prev_Alt_DA ...
    SR_Inc.Inc_Opt SR_Inc.Inc_Alt_DA ...
    SR_EL_opt.EL_Absenteeism SR_EL_alt_DA.EL_Absenteeism ...
    SR_EL_opt.EL_Presenteeism SR_EL_alt_DA.EL_Presenteeism ...
    SR_EL_opt.EL_Mortality SR_EL_alt_DA.EL_Mortality ...
    SR_EL_opt.EL_Total SR_EL_alt_DA.EL_Total]);
Table_2_Names = {'Elapsed_Time', 'Prev_Opt', 'Prev_Alt_DA', 'Inc_Opt', ...
    'Inc_Alt_DA', 'Abs_Opt', 'Abs_Alt_DA', 'Pre_Opt', 'Pre_Alt_DA', ...
    'Dec_Opt', 'Dec_Alt_DA', 'Tot_Opt', 'Tot_Alt_DA'};
Table_2.Properties.VariableNames = Table_2_Names;
Table_2{:,2:3} = Table_2{:,2:3}/100;

Table_3 = array2table([(1:T)' ...
    SR_Pre.Prev_Opt SR_Pre.Prev_Alt_DA SR_Pre.Prev_Alt_HB ...
    SR_Inc.Inc_Opt  SR_Inc.Inc_Alt_DA  SR_Inc.Inc_Alt_HB ...
    SR_EL_opt.EL_Absenteeism SR_EL_alt_DA.EL_Absenteeism ...
    SR_EL_alt_HB.EL_Absenteeism ...
    SR_EL_opt.EL_Presenteeism SR_EL_alt_DA.EL_Presenteeism ...
    SR_EL_alt_HB.EL_Presenteeism ...
    SR_EL_opt.EL_Mortality SR_EL_alt_DA.EL_Mortality ...
    SR_EL_alt_HB.EL_Mortality ...
    SR_EL_opt.EL_Total SR_EL_alt_DA.EL_Total SR_EL_alt_HB.EL_Total]);
Table_3_Names = {'Elapsed_Time', ...
    'Prev_Opt', 'Prev_Alt_DA', 'Prev_Alt_HB', ...
    'Inc_Opt',  'Inc_Alt_DA',  'Inc_Alt_HB', ...
    'Abs_Opt',  'Abs_Alt_DA',  'Abs_Alt_HB', ...
    'Pre_Opt',  'Pre_Alt_DA',  'Pre_Alt_HB', ...
    'Dec_Opt',  'Dec_Alt_DA',  'Dec_Alt_HB', ...
    'Tot_Opt',  'Tot_Alt_DA',  'Tot_Alt_HB'};
Table_3.Properties.VariableNames = Table_3_Names;
Table_3{:,2:4} = Table_3{:,2:4}/100;

T1 = strcat(MAT_Path, "Simulation_Results\Tables\", "Sim1_Table", ".xlsx");
T2 = strcat(MAT_Path, "Simulation_Results\Tables\", "Sim2_Table", ".xlsx");
T3 = strcat(MAT_Path, "Simulation_Results\Tables\", "Sim3_Table", ".xlsx");

% https://www.mathworks.com/matlabcentral/answers/91569-will-it-be-possible-to-use-xlswrite-to-write-formulas-to-the-excel-sheet-which-would-be-calculated-a
writetable(Table_1, T1); writetable(Table_2, T2); writetable(Table_3, T3);

% 012319: Now let's pre-format the workbooks in Excel so the user can begin
% using the tables immediately.
% https://www.mathworks.com/matlabcentral/answers/178989-open-edit-and-save-excel-file-with-matlab
% https://www.mathworks.com/matlabcentral/answers/96794-how-can-i-open-an-existing-excel-file-and-communicate-with-it-using-activex
% https://www.mathworks.com/matlabcentral/answers/216917-can-matlab-pre-format-individual-cells-when-writing-data-to-an-excel-spreadsheet
Excel = actxserver('Excel.Application'); Excel.Visible = 1; 
Workbooks = Excel.Workbooks; Workbook = Workbooks.Open(T1);
Sheet = Excel.Activesheet; Sheet.Range('B2:C71').NumberFormat = '0.0%';
Sheet.Range('D2:E71').NumberFormat = '0.0';
Sheet.Range('F2:M71').NumberFormat = '$0.0,,';
Workbook.Save; Workbook.Close;

Workbook = Workbooks.Open(T2);
Sheet = Excel.Activesheet; Sheet.Range('B2:C71').NumberFormat = '0.0%';
Sheet.Range('D2:E71').NumberFormat = '0.0';
Sheet.Range('F2:M71').NumberFormat = '$0.0,,';
Workbook.Save; Workbook.Close;

Workbook = Workbooks.Open(T3);
Sheet = Excel.Activesheet; Sheet.Range('B2:D71').NumberFormat = '0.0%';
Sheet.Range('E2:G71').NumberFormat = '0.0';
Sheet.Range('H2:S71').NumberFormat = '$0.0,,';
Workbook.Save; Workbook.Close;

Excel.Quit

%% Alternate Simulation Approach (DO NOT RUN)

% % 061519: The code below constitutes an alternative approach to carrying
% % out the microsimulation, preceding the current approach. It is
% % significantly less efficient than the approach I now use (~1500 seconds
% % vs. ~40 seconds for a single simulation involving Caucasian males with 
% % a high school diploma or equivalent), which I attribute largely to the 
% % need to loop over every simulated patient. The current approach is
% % vectorized over the cohort, and thus eliminates one of the loops
% % entirely.
% tic
% For every simulation period from t = 1 on (t = 0 is simply the starting
% state) ...
% for t = 2:T+1
% % ... generate a vector of random numbers - one for each simulated patient.
%     r = rand(k,1);
% %     Now, for each patient ...
%     for i = 1:length(s)
% %         ... assign the corresponding random number to a scalar value.
%         A = r(i,1);
% %         Identify the non-zero values in the row of the cumulative optimal
% %         transition probability matrix corresponding to the state the
% %         individual occupied in the previous period. B and C refer to the
% %         row and column indices associated with these values, while D
% %         refers to the values themselves.
%         [B_opt, C_opt, D_opt] = find(pstar_opt(SI_opt(i,t-1),:));
%         [B_alt_A, C_alt_A, D_alt_A] = ...
%             find(pstar_alt_A(SI_alt_A(i,t-1),:));
%         [B_alt_DA, C_alt_DA, D_alt_DA] = ...
%             find(pstar_alt_DA(SI_alt_DA(i,t-1),:));
%         [B_alt_HB, C_alt_HB, D_alt_HB] = ...
%             find(pstar_alt_HB(SI_alt_HB(i,t-1),:));
% %         Concatenate these three vectors together, as columns in a single
% %         array.
%         E_opt = horzcat(B_opt', C_opt', D_opt');
%         E_alt_A = horzcat(B_alt_A', C_alt_A', D_alt_A');
%         E_alt_DA = horzcat(B_alt_DA', C_alt_DA', D_alt_DA');
%         E_alt_HB = horzcat(B_alt_HB', C_alt_HB', D_alt_HB');
% %         The new state the individual occupies is the one associated with
% %         the first value (https://www.mathworks.com/help/matlab/ref/
% %         find.html) that is larger than the random value we drew at the 
% %         outset. It should be noted that these values will always fall in
% %         (0, 1), so we do not need to be concerned that no such value will
% %         meet this criterion, as might occur on the off-chance we drew a
% %         value exactly equal to 1 (https://www.mathworks.com/help/matlab/
% %         ref/rand.html). Also note that the first and second columns of E
% %         refer to the individual's initial and final states, respectively.
%         F_opt = E_opt(find(E_opt(:,3) > A, 1),2);
%         F_alt_A = E_alt_A(find(E_alt_A(:,3) > A, 1),2);
%         F_alt_DA = E_alt_DA(find(E_alt_DA(:,3) > A, 1),2);
%         F_alt_HB = E_alt_HB(find(E_alt_HB(:,3) > A, 1),2);
% %         Finally, we assign the individual's new state to the appropriate
% %         cell in SI.
%         SI_opt(i,t) = F_opt; SI_alt_A(i,t) = F_alt_A; 
%         SI_alt_DA(i,t) = F_alt_DA; SI_alt_HB(i,t) = F_alt_HB;
%     end
% end
% toc

% % 061519: This represented the second generation simulation mechanism. It
% % preserves the nested loop structure, but eliminated a lot of the
% % intermediate steps, making the code easier to read. I seem to recall
% % some improvement in operation time relative to the original approach, 
% % but would need to verify.
% tic
% SI_opt2(:,1) = s;
% for t = 2:T+1
%     for i = 1:length(s)
%         SI_opt2(i, t) = ...
%             find(pstar_opt(SI_opt2(i, t-1), :) > r(i, t-1), 1);
%     end
% end
% toc

% tic
% for t = 2:T+1
%     [~, SI_opt(:, t)] = ...
%         max(pstar_opt(SI_opt(:, t-1), :) > r(:, t-1), [], 2);
%     [~, SI_alt_A(:, t)] = ...
%         max(pstar_alt_A(SI_alt_A(:, t-1), :) > r(:, t-1), [], 2);
%     [~, SI_alt_DA(:, t)] = ...
%         max(pstar_alt_DA(SI_alt_DA(:, t-1), :) > r(:, t-1), [], 2);
%     [~, SI_alt_HB(:, t)] = ...
%         max(pstar_alt_HB(SI_alt_HB(:, t-1), :) > r(:, t-1), [], 2);
% end
% toc