function generate_brief_graphics 

%% Setting recurring values
base_dir = fileparts(mfilename('fullpath'));
x_aspect = 1.5; %#ok<NASGU>

% Setting years
year_values = 1989:3:2019;
year_values_short = 1989:6:2019;
n_years = length(year_values);


% Preallocating
data(11).weight   = [];



%% Importing Data
fprintf('Importing data.\n')
for iy = 1:n_years
    
    fprintf('Importing year %i.\n', year_values(iy))
    
    data_file = ['SCFP',num2str(year_values(iy)),'.csv'];


    opts = detectImportOptions(data_file);

    opts.SelectedVariableNames = {'WGT'};
    data(iy).weight            = readmatrix(data_file,opts);
    
    % Net Worth
    
    opts.SelectedVariableNames = {'NETWORTH'};
    data(iy).networth          = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'ASSET'};
    data(iy).assets            = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'DEBT'};
    data(iy).debt            = readmatrix(data_file,opts);
    
    % Assets
    
    opts.SelectedVariableNames = {'FIN'};
    data(iy).fin               = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'NFIN'};
    data(iy).nfin              = readmatrix(data_file,opts);
    
    % Financial
    
    opts.SelectedVariableNames  = {'LIQ'};
    data(iy).liquid             = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'CDS'};
    data(iy).cds                = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'NMMF'};
    data(iy).pooled_investments = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'SAVBND'};
    data(iy).savings_bonds      = readmatrix(data_file,opts);

    opts.SelectedVariableNames  = {'STOCKS'};
    data(iy).stocks             = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'BOND'};
    data(iy).bonds              = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'CASHLI'};
    data(iy).life_insurance     = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'OTHMA'};
    data(iy).other_managed      = readmatrix(data_file,opts);

    opts.SelectedVariableNames  = {'RETQLIQ'};
    data(iy).qliquid_retirement = readmatrix(data_file,opts);

    opts.SelectedVariableNames  = {'OTHFIN'};
    data(iy).other_fin          = readmatrix(data_file,opts);
    
    
    
    % Additional Financial Items
    
    opts.SelectedVariableNames  = {'STMUTF'};
    data(iy).stock_mf           = readmatrix(data_file,opts);

    opts.SelectedVariableNames  = {'COMUTF'};
    data(iy).combination_mf     = readmatrix(data_file,opts);


    

    % Non-financial
    
    opts.SelectedVariableNames = {'VEHIC'};
    data(iy).vehicle           = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'HOUSES'};
    data(iy).houses            = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'ORESRE'};
    data(iy).other_houses      = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'NNRESRE'};
    data(iy).nonresidential    = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames = {'BUS'};
    data(iy).business          = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'OTHNFIN'};
    data(iy).other_nfin          = readmatrix(data_file,opts);
    
    
    % Debt
    
    opts.SelectedVariableNames  = {'MRTHEL'};
    data(iy).residential_debt   = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'RESDBT'};
    data(iy).other_res_debt     = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'OTHLOC'};
    data(iy).other_credit       = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'CCBAL'};
    data(iy).credit_card        = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'INSTALL'};
    data(iy).installment        = readmatrix(data_file,opts);
    
    opts.SelectedVariableNames  = {'ODEBT'};
    data(iy).other_debt         = readmatrix(data_file,opts);
    
    
    % Income
    
    opts.SelectedVariableNames = {'INCOME'};
    data(iy).income            = readmatrix(data_file,opts);

end



%% Calculating Ginis and shares

% Preallocating structures
gini(11).wealth = [];
gini(11).income = [];
fprintf('Calculating ginis.\n')
for iy = 1:n_years
    
    gini(iy).wealth = calculate_gini( data(iy).networth, data(iy).weight );
    gini(iy).income = calculate_gini( data(iy).income  , data(iy).weight );
    
end



% Calculating percentiles and corresponding conditional median values

xtiles = 0:0.1:1;
wealth_ginis = zeros(1,n_years);
income_ginis = zeros(1,n_years);

wealth_wealth_shares = [];
income_income_shares = [];

wealth_wealth_medians = [];
income_income_medians = [];

fprintf('Calculating shares.\n')
for iy = 1:n_years
    
    wealth_ginis(iy) = gini(iy).wealth ;
    income_ginis(iy) = gini(iy).income ;
    
    [ww_shares, ww_medians] = calculate_xtile(xtiles, data(iy).networth, data(iy).networth, data(iy).weight, 'method2');
    wealth_wealth_shares  = [wealth_wealth_shares ; ww_shares ]; %#ok<*AGROW>
    wealth_wealth_medians = [wealth_wealth_medians; ww_medians];    
    
    [ii_shares, ii_medians] = calculate_xtile(xtiles, data(iy).income  , data(iy).income  , data(iy).weight, 'method2');
    income_income_shares  = [income_income_shares ; ii_shares ];
    income_income_medians = [income_income_medians; ii_medians];
    
end
fprintf('Generating graphics.\n')


% Ginis with shares


% Wealth 
figure
subplot(2,2,1)
plot(year_values, wealth_ginis, 'LineWidth', 2)
grid on
xlim([year_values(1),year_values(end)])
ylim([0.78,0.865])
xticks(year_values_short)
ylabel('Wealth Gini Coefficient')
pbaspect([4.0,2.5,1])

subplot(2,2,3)
plot(year_values,100.*[sum(wealth_wealth_shares(:,1:5),2),...
                       sum(wealth_wealth_shares(:,6:9),2),...
                       sum(wealth_wealth_shares(:, 10),2)]', 'LineWidth', 2)
grid on
xlim([year_values(1),year_values(end)])
ylim(100.*[0,0.9])
xticks(year_values_short)
ylabel('% of Total Wealth')
pbaspect([4.0,2.5,1])
legend('Bottom 50%','Top 50-90%', 'Top 10%','location','southoutside')


% Income

subplot(2,2,2)
plot(year_values, income_ginis, 'LineWidth', 2)
grid on
xlim([year_values(1),year_values(end)])
ylim([0.49,0.61])
xticks(year_values_short)
ylabel('Income Gini Coefficient')
pbaspect([4.0,2.5,1])

subplot(2,2,4)
plot(year_values,100.*[sum(income_income_shares(:,1:5),2),...
                       sum(income_income_shares(:,6:9),2),...
                       sum(income_income_shares(:, 10),2)]', 'LineWidth', 2)
grid on
xlim([year_values(1),year_values(end)])
ylim(100.*[0.1,0.52])
xticks(year_values_short)
ylabel('% of Total Income')
pbaspect([4.0,2.5,1])
legend('Bottom 50%','Top 50-90%', 'Top 10%','location','southoutside')

print('-djpeg',fullfile(base_dir,'gini_shares_combined.jpeg'), '-r300');
print('-djpeg',fullfile(base_dir,             'figure1.jpeg'), '-r300');





%% Wealth by wealth

labels= categorical({'0-10','10-20','20-30','30-40','40-50','50-60','60-70','70-80','80-90','90-100'});

% Wealth shares by wealth deciles and median wealth by wealth decile
[wealth_wealth_shares_2016, wealth_wealth_medians_2016] = calculate_xtile(xtiles, data(10).networth, data(10).networth, data(10).weight);
[wealth_wealth_shares_2019, wealth_wealth_medians_2019] = calculate_xtile(xtiles, data(11).networth, data(11).networth, data(11).weight);


% Changes in the wealth distribution as a function of wealth
figure
subplot(2,2,1)
bar(labels,100.*[wealth_wealth_shares_2016;wealth_wealth_shares_2019])
ylim([0,82])
legend('2016','2019','location','northwest')
ylabel('% of Total Wealth')
xlabel('Wealth Percentiles')
pbaspect([4.0,2.5,1])

subplot(2,2,3)
bar(labels,100.*(wealth_wealth_shares_2019-wealth_wealth_shares_2016))
ylim([-0.65,0.3])
ylabel('Change, 2016-2019')
xlabel('Wealth Percentiles')
pbaspect([4.0,2.5,1])

subplot(2,2,2)
bar(labels,[wealth_wealth_medians_2016;wealth_wealth_medians_2019])
legend('2016','2019','location','northwest')
ylabel('Median Wealth')
xlabel('Wealth Percentiles')
pbaspect([4.0,2.5,1])

subplot(2,2,4)
bar(labels,100.*(-1+wealth_wealth_medians_2019./wealth_wealth_medians_2016))
ylim([0,63])
ylabel('% Change, 2016-2019')
xlabel('Wealth Percentiles')
pbaspect([4.0,2.5,1])

print('-djpeg',fullfile(base_dir,'wealth_wealth_shares.jpeg'), '-r300');
print('-djpeg',fullfile(base_dir,             'figure2.jpeg'), '-r300');



%% Distributions by type of net worth

% Wealth shares by wealth deciles and median wealth by wealth decile
[fin_wealth_shares_2016] = calculate_xtile(xtiles, data(10).networth, data(10).fin, data(10).weight);
[fin_wealth_shares_2019] = calculate_xtile(xtiles, data(11).networth, data(11).fin, data(11).weight);

[nfin_wealth_shares_2016] = calculate_xtile(xtiles, data(10).networth, data(10).nfin, data(10).weight);
[nfin_wealth_shares_2019] = calculate_xtile(xtiles, data(11).networth, data(11).nfin, data(11).weight);


figure
subplot(2,2,1)
bar(labels,100.*[fin_wealth_shares_2016;fin_wealth_shares_2019])
ylim(100.*[0,0.8])
xlabel('Wealth Percentiles')
ylabel('% of Total Fin. Wealth')
legend('2016','2019','location','northwest')
pbaspect([4.0,2.5,1])

subplot(2,2,3)
bar(labels,100.*(fin_wealth_shares_2019-fin_wealth_shares_2016) )
ylim([-2,1])
xlabel('Wealth Percentiles')
ylabel('Change, 2016-2019')
pbaspect([4.0,2.5,1])

subplot(2,2,2)
bar(labels,100.*([nfin_wealth_shares_2016;nfin_wealth_shares_2019]))
ylim(100.*[0,0.7])
legend('2016','2019','location','northwest')
xlabel('Wealth Percentiles')
ylabel('% of Total Nonfin. Wealth')
pbaspect([4.0,2.5,1])

subplot(2,2,4)
bar(labels,100.*(nfin_wealth_shares_2019-nfin_wealth_shares_2016))
xlabel('Wealth Percentiles')
ylabel('Change, 2016-2019')
pbaspect([4.0,2.5,1])

print('-djpeg',fullfile(base_dir,'fin_nfin_wealth.jpeg'), '-r300');
print('-djpeg',fullfile(base_dir,        'figure3.jpeg'), '-r300');



%% Income by income

% Wealth shares by income deciles and median wealth by income decile
[income_income_shares_2016, income_income_medians_2016] = calculate_xtile(xtiles, data(10).income, data(10).income, data(10).weight);
[income_income_shares_2019, income_income_medians_2019] = calculate_xtile(xtiles, data(11).income, data(11).income, data(11).weight);

% Changes in the income distribution, as a function of income
figure
subplot(2,2,1)
bar(labels,100.*[income_income_shares_2016;income_income_shares_2019])
ylim([0,60])
legend('2016','2019','location','northwest')
ylabel('% of Total Income')
xlabel('Income Percentiles')
pbaspect([4.0,2.5,1])

subplot(2,2,3)
bar(labels,100.*(income_income_shares_2019-income_income_shares_2016))
ylim([-4,1.5])
yticks(-4:1)
ylabel('Changes, 2016 to 2019')
xlabel('Income Percentiles')
pbaspect([4.0,2.5,1])

subplot(2,2,2)
bar(labels,[income_income_medians_2016;income_income_medians_2019])
legend('2016','2019','location','northwest')
ylabel('Median Income')
xlabel('Income Percentiles')
pbaspect([4.0,2.5,1])

subplot(2,2,4)
bar(labels,100.*(-1+income_income_medians_2019./income_income_medians_2016))
ylabel('% Change, 2016-2019')
xlabel('Income Percentiles')
pbaspect([4.0,2.5,1])

print('-djpeg',fullfile(base_dir,'income_income_shares.jpeg'), '-r300');
print('-djpeg',fullfile(base_dir,'figure4.jpeg'), '-r300');



%% Additional information

% Median wealth, entire sample
total_median_wealth_2016 = find_median( data(10).networth , data(10).weight );
total_median_wealth_2019 = find_median( data(11).networth , data(11).weight );

fprintf('Median wealth in 2016: $%0.4f.\n', total_median_wealth_2016)
fprintf('Median wealth in 2019: $%0.4f.\n', total_median_wealth_2019)

fprintf('Change in median wealth 2016-2019: %0.4f%%.\n', 100*(-1+total_median_wealth_2019/total_median_wealth_2016) )

% Median income, entire sample
total_median_income_2016 = find_median( data(10).income , data(10).weight );
total_median_income_2019 = find_median( data(11).income , data(11).weight );

fprintf('Median income in 2016: $%0.4f.\n', total_median_income_2016)
fprintf('Median income in 2019: $%0.4f.\n', total_median_income_2019)

fprintf('Change in median income 2016-2019: %0.4f%%.\n', 100*(-1+total_median_income_2019/total_median_income_2016) )



% Nonfinancial share of total wealth
total_nonfinancial_2019 = sum( data(11).weight .* data(11).nfin     );
total_networth_2019     = sum( data(11).weight .* data(11).networth );
fprintf('Nonfinancial share of total wealth in 2019:          %0.4f%%.\n', 100*total_nonfinancial_2019/total_networth_2019)

% Housing share of total wealth
total_houses_2019   = sum( data(11).weight .* data(11).houses   );
total_networth_2019 = sum( data(11).weight .* data(11).networth );
fprintf('Housing share of total wealth in 2019:               %0.4f%%.\n', 100*total_houses_2019/total_networth_2019)

% Housing share of total nonfinancial wealth
fprintf('Housing share of total nonfinancial wealth in 2019:  %0.4f%%.\n', 100*total_houses_2019/total_nonfinancial_2019)


end



function [gini,domain,lorenz_curve] = calculate_gini( data, weight  )


if ~exist('weight','var') || isempty(weight)
    % Equal weights if none supplied
    weight = ones(1,length(data))';
end

[data, indexes] = sort   ( data    )  ;
weight          = weight ( indexes )  ;
weight          = weight./sum(weight)';
pop_share       = cumsum( weight )    ;
data_share      = cumsum( weight.* data) ./sum( weight.*data ) ;


pop_share  = [0, pop_share'];
data_share = [0,data_share'];




f = @(x) interp1( pop_share, data_share, x, 'linear');

B = integral( f, 0, 1, 'AbsTol', 1e-12 );

A = 0.5 - B;

gini = A/(A+B);

domain       = pop_share  ;
lorenz_curve = data_share ;




end



function median_value = find_median( data, weight )


if ~exist('weight','var') || isempty(weight)
    weight = ones(size(data));
end

[data, indexes] = sort ( data );
weight = weight ( indexes );

cumulative_densities = cumsum( weight ./ sum( weight ) );

median_upper_index = find( cumulative_densities>0.5, 1, 'first' );

lower_prob = cumulative_densities( median_upper_index - 1 );
upper_prob = cumulative_densities( median_upper_index     );

median_value = data( median_upper_index - 1 ) * ( (upper_prob - 0.5)/(upper_prob - lower_prob) ) + ...
               data( median_upper_index     ) * ( (0.5 - lower_prob)/(upper_prob - lower_prob) );



end


function [shares, medians] = calculate_xtile(xtiles, data1, data2, weight, method)

% Gives shares and corresponding median values of data2 as a function of
% data1. 'xtiles' must give *endpoints* including zero for lower bound.

% This method has some small error corresponding to the imprecision of
% discrete index values. Can be reduced by local weighting.

if ~exist('weight','var')||isempty(weight)
    weight = zeros(1,length(data1));
end


if ~exist('method','var')
    method = 'method1';
end

[~, indexes] = sort ( data1 );
data2  = data2  ( indexes );
weight = weight ( indexes );
weight = weight./sum(weight)'; % Tranformed to probability measures

pop_share   = cumsum( weight  ); % recall: weight sorted according to data1
data2_share = cumsum(  weight.* data2) ./sum(  weight.*data2 );

% Finding upper indexes of xtile groups
n_xtiles = length(xtiles);    
    
shares  = zeros(1,n_xtiles-1);
medians = zeros(1,n_xtiles-1);

lower_index = 0;

% Adjusting for numerical precision
pop_share  (end) = 1;
data2_share(end) = 1;
for ix = 1:n_xtiles-1


    lower_index = lower_index+1;
    upper_index = find( pop_share <= xtiles(ix+1), 1, 'last');

    xtile_subset = lower_index:upper_index ;


    medians(ix) = find_median( data2(xtile_subset), weight(xtile_subset) ); 


    if strcmp(method,'method1')
        shares (ix) = data2_share(upper_index) - data2_share(lower_index);
    elseif strcmp(method,'method2')
        shares(ix) =  interp1([0,pop_share'  ],[0,data2_share'],xtiles(ix+1)) - ...
                      interp1([0,pop_share'  ],[0,data2_share'],xtiles(ix  ));
    end

    lower_index = upper_index;
end

    



end