%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Analisis de la Economia %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; clc; close;

filename = 'Data.xlsx';
sheet = 'SAM_Desagregada';
SAM = xlsread(filename, sheet);

SAM(1,:) = [];  % Elimina la primera fila
SAM(:,1) = [];  % Elimina la primera columna

SAM(isnan(SAM)) = 0;

%% Estadistica Descriptiva

% Pagos Salariales
% Suma de todos los pagos salariales que se distribuyen entre los quintiles
total_wage_payments_to_quintiles = sum(SAM(28:32,25));

% Cálculo de las proporciones para cada quintil
proportions_from_wages = SAM(28:32,25) / total_wage_payments_to_quintiles;

% Suma total de las proporciones
sum_proportions = sum(proportions_from_wages);

% Mostrar el resultado
disp(['Suma total de las proporciones: ', num2str(sum_proportions)]);
disp('Proporciones distribuidas entre los quintiles:');
disp(proportions_from_wages);


% Pagos por Capital
% Suma de todos los por Capital que se distribuyen entre los quintiles
total_capital_payments_to_quintiles = sum(SAM(28:32,26));

% Cálculo de las proporciones para cada quintil
proportions_from_capital = SAM(28:32,26) / total_capital_payments_to_quintiles;

% Suma total de las proporciones
sum_proportions = sum(proportions_from_capital);

% Mostrar el resultado
disp(['Suma total de las proporciones: ', num2str(sum_proportions)]);
disp('Proporciones distribuidas entre los quintiles:');
disp(proportions_from_capital);

% Transferencias de Gobierno

% Extraer las transferencias del gobierno a cada quintil
government_transfers = SAM(28:32, 33);

% Calcular la proporción de transferencias a cada quintil
transfer_proportions = government_transfers / sum(government_transfers);

% Crear una tabla para visualizar las proporciones y cantidades absolutas
quintile_names = {'Quintile 1', 'Quintile 2', 'Quintile 3', 'Quintile 4', 'Quintile 5'};
T = table(quintile_names', government_transfers, transfer_proportions, 'VariableNames', {'Quintile', 'Amount', 'Proportion'});
disp(T);


% Consumo Hogares

% Extract the consumption of each quintile for all the productive sectors
consumption_by_quintile = SAM(1:24, 28:32);

% Pre-allocate a matrix to store the top 4 sectors for each quintile
top_sectors = zeros(4, 5); 

% For each quintile, find the top 4 sectors of consumption
for q = 1:5
    [~, sorted_indices] = sort(consumption_by_quintile(:, q), 'descend');
    top_sectors(:, q) = sorted_indices(1:4);
end

% Define the sector names for ease of reading
sector_names = {
    'Agricultural forestry and fishing';
    'Mining';
    'Manufacturing industry';
    'Electricity, gas, water, and waste management';
    'Building';
    'Commerce, hotels, and restaurants';
    'Transport, communications, and information services';
    'Financial intermediation';
    'Real estate and housing services';
    'Business services';
    'Personal services';
    'Public administration';
    'Activity—Agricultural forestry and fishing';
    'Activity—Mining';
    'Activity—Manufacturing Industry';
    'Activity—Electricity, gas, water, and waste management';
    'Activity—Building';
    'Activity—Commerce, hotels, and restaurants';
    'Activity—Transport, communications, and information services';
    'Activity—Financial intermediation';
    'Activity—Real estate and housing services';
    'Activity—Business services';
    'Activity—Personal services';
    'Activity—Public administration'
};

% Display the top 4 sectors for each quintile
for q = 1:5
    disp(['Top 4 sectors for Quintile ', num2str(q), ':']);
    for s = 1:4
        disp(['  ', sector_names{top_sectors(s, q)}]);
    end
    disp('-----------------------------');
end


%% Ordenar la matriz: Sacar fila y columna de gobierno
totCtas = 41;
cExo = [33, 41];                          % ctas exo: 33 gobierno, 41 err
cEndo = setdiff(1:totCtas, cExo);
nCtasEx = length(cExo);
nCtasEn = length(cEndo);
Dnew = zeros(totCtas);

% Resultado: Ctas exo en pos 40 y 41
Dnew(1:nCtasEn, 1:nCtasEn) = SAM(cEndo, cEndo);
for i = 1:nCtasEx
    Dnew(nCtasEn + i, 1:nCtasEn) = SAM(cExo(i), cEndo);
    Dnew(nCtasEn + i, nCtasEn + 1:totCtas) = SAM(cExo(i), cExo);
    Dnew(1:nCtasEn, nCtasEn + i) = SAM(cEndo, cExo(i));
    Dnew(nCtasEn + 1:totCtas, nCtasEn + i) = SAM(cExo, cExo(i));
end

xf = sum(Dnew, 2); % Suma por columnas

AnularColumna = find(xf == 0);

% Matriz A de coeficientes técnicos
A = Dnew * diag(1./xf);
A(:, AnularColumna) = 0;

Ann = A(1:nCtasEn, 1:nCtasEn);
Ank = A(1:nCtasEn, nCtasEn+1:totCtas);
Akn = A(nCtasEn+1:totCtas, 1:nCtasEn);
Akk = A(nCtasEn+1:totCtas, nCtasEn+1:totCtas);
yn = xf(1:nCtasEn);
yk = xf(nCtasEn+1:totCtas);

% Efectos totales
if abs(sum(yn - (Ann * yn + Ank * yk))) < 1e-6
    disp('Test de consistencia 1 OK');
else 
    disp('No')
end

Ma = inv(eye(nCtasEn) - Ann);
x = Ank * yk;

if abs(sum(Ma * x - (Ann * yn + Ank * yk))) < 1e-6
    disp('Test de consistencia 2 OK');
end

Ec2 = Ma * x;

% Matriz redistributiva
e = ones(nCtasEn, 1);
zn = yn ./ (e' * yn);
R1 = 1 / (e' * yn);
R2 = eye(nCtasEn) - zn * e';
R = R1 * R2 * Ma;

% Impacto: Tabla 5, col 3; Garrido y Morales (2023)
dx1 = [1; 0];
dx = Ank * dx1;
impacto = R * dx; 

% Efectos de difusión del aumento en una unidad de la cuenta exógena
dYn = Ma * dx; % Muestra el efecto difusión de los hogares, igual a tesis
               % Tabla 5, col 4; Garrido y Morales (2023)

% sectores con mayor encadenamiento:
% Tomar solamente los valores de las cuentas del 13 al 24
subsetValues = dYn(13:24);

% Ordenar esos valores
[sortedSubsetValues, sortedSubsetIndices] = sort(subsetValues, 'descend');

% Tomar los primeros 5 valores y sus índices
top5SubsetValues = sortedSubsetValues(1:5);
top5SubsetIndices = sortedSubsetIndices(1:5) + 12; % Ajustar los índices para reflejar la posición en dYn_prod

% Mostrar los resultados
disp('Los 5 mayores valores de las cuentas 13 a 24 y sus índices son:')
disp([top5SubsetValues, top5SubsetIndices]);

%% Network: Sectores Productivos

% Submatriz de coef tec para solo los sectores productivos (1 a 24)
A_prod = A(1:24, 1:24);

% Efectos multiplicadores para solo los sectores productivos (1 a 24)
dYn_prod = dYn(1:24);

% Convertir A_prod en una matriz simétrica
adjMatrix_prod = A_prod | A_prod';

% Crear el gráfico
G_prod = graph(adjMatrix_prod);

% Normalizar y escalar el tamaño de los nodos basado en el efecto multiplicador de los sectores productivos
nodeSize_prod = 20 * (dYn_prod - min(dYn_prod)) / (max(dYn_prod) - min(dYn_prod)) + 10;

% Definir los colores de los nodos
nodeColors = [repmat([0.5 0.5 0.5], 12, 1); repmat([0 0 0], 12, 1)]; % Gris para cuentas 1 a 12 y negro para cuentas 13 a 24

% Dibujar el gráfico
figure('Position', [10 10 1200 800]);
h_prod = plot(G_prod, 'MarkerSize', nodeSize_prod, 'NodeColor', nodeColors, 'NodeLabel', {});
layout(h_prod, 'force', 'UseGravity', true, 'Iterations', 150);
h_prod.EdgeColor = [0 0 0]; % Aristas en color negro

% Añadir etiquetas con colores apropiados
for i = 1:12 % Para nodos grises, etiquetas negras
    text(h_prod.XData(i), h_prod.YData(i), num2str(i), 'Color', 'black', 'FontSize', 10, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
end
for i = 13:24 % Para nodos negros, etiquetas blancas
    text(h_prod.XData(i), h_prod.YData(i), num2str(i), 'Color', 'white', 'FontSize', 10, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
end

title('Red de Sectores Productivos (Sectores 1 a 24)');
xlabel('Sectores Productivos');
ylabel('Conexiones');