clear all;
clc;
load("workspace_2020_07_01.mat");

months_long = repmat([1 4 7 10]',[ceil(134/4),1]);
months_long = months_long(1:end-3);
years = floor(1987:0.25:2020)';
Time_long=datetime(years,months_long,1);
Tlong = length(Time_long);
rng(123);
totalfinancialaccounts_matrix = vertcat(financialaccounts_matrix,zeros(1,T));
Liabilities = zeros(N,N,T);

for t=1:T
    for i=0:(N-1)
        for j=1:N
            Liabilities(i+1,j,t) = totalfinancialaccounts_matrix(N*i+j,t);
        end
    end
end



NetWorth_enhancedaccounts = zeros(N,T);
NetWorth_enhancedaccounts(:,:) = squeeze(sum(Liabilities,1))-squeeze(sum(Liabilities,2));    %assets - liabilities

NetWorth_enhancedaccounts(7,:)=NetWorth_enhancedaccounts(7,:)+NetWorth_enhancedaccounts(8,:);

NetWorth_enhancedaccounts(9,:)=NetWorth_enhancedaccounts(9,:)+NetWorth_enhancedaccounts(10,:);
NetWorth_enhancedaccounts([5,8,10,12],:)=[];

Flowoffundsbasednetworth(:,6)=Flowoffundsbasednetworth(:,6)+Flowoffundsbasednetworth(:,7);
Flowoffundsbasednetworth(:,8)=Flowoffundsbasednetworth(:,8)+Flowoffundsbasednetworth(:,9);
Flowoffundsbasednetworth(:,[7,9])=[];
Flowoffundsbasednetworth=Flowoffundsbasednetworth';

Flowoffundsbasedassets(:,6)=Flowoffundsbasedassets(:,6)+Flowoffundsbasedassets(:,7);
Flowoffundsbasedassets(:,8)=Flowoffundsbasedassets(:,8)+Flowoffundsbasedassets(:,9);
Flowoffundsbasedassets(:,[7,9])=[];
Flowoffundsbasedassets=Flowoffundsbasedassets';
totalassets = Flowoffundsbasedassets(:,1:Tlong);


Flowoffundsbasedliab(:,6)=Flowoffundsbasedliab(:,6)+Flowoffundsbasedliab(:,7);
Flowoffundsbasedliab(:,8)=Flowoffundsbasedliab(:,8)+Flowoffundsbasedliab(:,9);
Flowoffundsbasedliab(:,[7,9])=[];
Flowoffundsbasedliab=Flowoffundsbasedliab';
totalliabilities = Flowoffundsbasedliab(:,1:Tlong);
leverageratio = totalliabilities./totalassets;
debttoratio = totalliabilities./(totalliabilities+totalassets);


Liabilities_big = Liabilities;
%% put central government and local government together
row_hulp = Liabilities(8,:,:) + Liabilities(7,:,:);
Liabilities_cut = zeros(N-1,N,T);
Liabilities_cut(1:6,:,:) = Liabilities(1:6,:,:);
Liabilities_cut(7,:,:) = row_hulp;
Liabilities_cut(8:N-1,:,:) = Liabilities(9:N,:,:);

Liabilities = zeros(N-1,N-1,T);
Liabilities(:,1:6,:) = Liabilities_cut(:,1:6,:);
Liabilities(:,7,:) = Liabilities_cut(:,7,:)+Liabilities_cut(:,8,:);
Liabilities(:,8:N-1,:) = Liabilities_cut(:,9:N,:);

N=N-1;
%% remove unknown
Liabilities= Liabilities(1:N-1,1:N-1,:);
Liabilities(Liabilities<0)=0;
N=N-1;
%% put HH & NPO together
row_hulp = Liabilities(8,:,:) + Liabilities(9,:,:);
Liabilities_cut = zeros(N-1,N,T);
Liabilities_cut(1:7,:,:) = Liabilities(1:7,:,:);
Liabilities_cut(8,:,:) = row_hulp;
Liabilities_cut(9:N-1,:,:) = Liabilities(10:N,:,:);

Liabilities = zeros(N-1,N-1,T);
Liabilities(:,1:7,:) = Liabilities_cut(:,1:7,:);
Liabilities(:,8,:) = Liabilities_cut(:,8,:)+Liabilities_cut(:,9,:);
Liabilities(:,9:N-1,:) = Liabilities_cut(:,10:N,:);

N=N-1;
%% remove NMMF because it is a subset of OFI 
Liabilities_cut = zeros(N-1,N,T);
Liabilities_cut(1:4,:,:) = Liabilities(1:4,:,:);
Liabilities_cut(5:N-1,:,:)=Liabilities(6:N,:,:);
Liabilities=zeros(N-1,N-1,T);
Liabilities(:,1:4,:) = Liabilities_cut(:,1:4,:);
Liabilities(:,5:N-1,:) = Liabilities_cut(:,6:N,:);
N=N-1;
%%
Liabilities(1,7,32)=0.1;   %this is done to match the new series with the data downloaded in 2016
%


totalassets_enh = squeeze(sum(Liabilities,1));
totalliabilities_enh = squeeze(sum(Liabilities,2));
leverageratio_enh = totalliabilities./totalassets;
debttoratio_enh = totalliabilities./(totalliabilities+totalassets);





Sectors=["Public Corporations", "Private Non-Financial Corporations", "Monetary Financial Institutions", "Other Financial Intermediaries", "Insurance Corporations and Pension Funds", "Government", "Households & Non-Profit Institutions Serving Households", "Rest of the World"];
Sectors_ab = ["PC", "PNFC", "MFI", "OFI", "ICPF", "G", "HH+NPISH", "RoW"];
Sectors_semi_ab = ["Public Corp.", "Private Non-Fin.", "Monetary Fin. Inst.", "Other Fin. Interm.", "Insurance & Pensions", "Government", "Households & NPO", "Rest of the World"];

thresholds = 0.05:0.01:0.8;
%thresholds=0.1:0.02:5;

% kShells_T = ones(N,T);
pijs_T = ones(N,N,T);
node_strength = zeros(N,T);

for t=1:T
    network = Liabilities(:,:,t);
    [pijs] = AnalyzeNetwork(network, thresholds);
    pijs_T(:,:,t) = pijs;
end

%%
NetWorth_no_unkown = zeros(N,T);
% %%
NetWorth_no_unkown(:,:) = squeeze(sum(Liabilities,1))-squeeze(sum(Liabilities,2));    %assets - liabilities
%
NetWorthGrowth = 100*(Flowoffundsbasednetworth(1:N, 2:Tlong) - Flowoffundsbasednetworth(1:N, 1:Tlong-1))./ Flowoffundsbasednetworth(1:N, 1:Tlong-1);


Downturns = zeros(N,Tlong-1);
for i=1:N
    for t=1:(T-1)
        if NetWorthGrowth(i,t)<0   %quantile(NetWorthGrowth(i,:),.33)
            Downturns(i,t) = 1;
        end
    end
end


Two_level = Downturns;

for i=1:N
    for t=2:(Tlong-1)
        
        if Downturns(i,t-1)==1&&Downturns(i,t)==1
            Two_level(i,t)=2;
        end
    end
end

T1 = table(Time_long(2:Tlong), Downturns.');
T3 = table(Time_long(2:Tlong), Two_level.');

writetable(T1, "One level UK.csv");
writetable(T3, "Two level UK.csv");

numberOfCrises = zeros(1,N);
numberOfCrisesPeriods= zeros(1,N);
numberOfCrisesRecov = zeros(1,N);

for i=1:N
    for t=1:(Tlong-1)
        if Two_level(i,t)==2
            if t>1 && Two_level(i,t-1)==1
                numberOfCrises(i) = numberOfCrises(i)+1;
            end
            numberOfCrisesPeriods(i) = numberOfCrisesPeriods(i) + 1;
        end
        if Two_level(i,t)==0
            if t>1 && Two_level(i,t-1)==2
                numberOfCrisesRecov(i) = numberOfCrisesRecov(i)+1;
            end
        end
    end
end

MeanCrisisLength = numberOfCrisesPeriods./numberOfCrises;

number_of_downturns = zeros(1,N);
recovered_of_downturn = zeros(1,N);

for i=1:N
    for t=1:(Tlong-1)
        if Two_level(i,t)==1
            number_of_downturns(i) = number_of_downturns(i)+1;
            if t<(T-1)&&Two_level(i,t+1)==0
                recovered_of_downturn(i)=recovered_of_downturn(i)+1;
            end
        end
    end
end

recov_probability_downturn = recovered_of_downturn./number_of_downturns;

init_values = 0.1*ones(1,2*N+N*(N-1)+1);
init_values((N*(N-1)+1):(N*N)) = recov_probability_downturn;
init_values((N*N+1):N*(N+1)) = 1./MeanCrisisLength;
Init = table(init_values.');
writetable(Init, "init values UK.csv");

init_values_downturn = 0.1*ones(1,N+N*(N-1)+1);
init_values_downturn((N*(N-1)+1):(N*N)) = 1./MeanCrisisLength;
Init_downturn = table(init_values_downturn.');
writetable(Init_downturn, "init values downturn UK.csv");


% Compute Prior values -> gebruik pijs_T.
prior_pijs_alpha = zeros(N,N);
prior_pijs_beta = zeros(N,N);

prior_pijs_mean = zeros(N,N);
prior_pijs_se = zeros(N,N);
scale_factor_pijs = 1.0;

pijs_T = pijs_T*scale_factor_pijs;

for i=1:N
    for j=1:N
        if i~=j
            phat = betafit(reshape(pijs_T(i,j,:),1,T));
            prior_pijs_alpha(i,j) = phat(1);
            prior_pijs_beta(i,j) = phat(2);
            prior_pijs_mean(i,j) = prior_pijs_alpha(i,j)/(prior_pijs_alpha(i,j)+prior_pijs_beta(i,j));
            prior_pijs_se(i,j) = sqrt((prior_pijs_alpha(i,j)*prior_pijs_beta(i,j))/((prior_pijs_alpha(i,j)+prior_pijs_beta(i,j))^2*(prior_pijs_alpha(i,j)+prior_pijs_beta(i,j)+1))); 
        end
    end
end

Number_of_successful_infections = zeros(N,N);
Number_of_failed_infections = zeros(N,N);

for i=1:N
    for j=1:N
        if i~=j
            for t=2:Tlong-1
                if and(and(Two_level(i,t)==1, Two_level(i,t-1)==0), Two_level(j,t-1)>=1)
                    Number_of_successful_infections(j,i)=Number_of_successful_infections(j,i)+1;
                elseif and(and(Two_level(i,t)==0, Two_level(i,t-1)==0), Two_level(j,t-1)>=1)
                    Number_of_failed_infections(j,i)=Number_of_failed_infections(j,i)+1;
                end
            end
        end
    end
end

%%

pijs_mode = zeros(N,N);
pijs_se = zeros(N,N);
pijs_mode_mle = zeros(N,N);
pijs_se_mle = zeros(N,N);

m=22000;

pi_sim = zeros(m,N);
qi_sim = zeros(m,N);
for i=1:N
    qi_sim(:,i) = betarnd(numberOfCrisesRecov(i)+1, numberOfCrisesPeriods(i) - numberOfCrisesRecov(i)+1, 1,m);
    pi_sim(:,i) = betarnd(recovered_of_downturn(i)+1, number_of_downturns(i)-recovered_of_downturn(i)+1, 1,m);
end
R0c = zeros(N,m);
R0d = zeros(N,m);
R0c_mle = zeros(N,m);
R0d_mle = zeros(N,m);
pijs_sim = zeros(N,N,m);
pijs_sim_mle = zeros(N,N,m);


qi_sim_prior = betarnd(1, 1, 1,m);
pi_sim_prior = betarnd(1, 1, 1,m);
R0c_prior = zeros(N,m);
R0d_prior = zeros(N,m);
pijs_sim_prior = zeros(N,N,m);



sigma_n  = 0.02;
sigma_z = 0.7;
n_0 = 0.1;
p_0 = zeros(N,N);
for i=1:N
    for j=1:N
        p_0(i,j) = betarnd(Number_of_successful_infections(i,j)+prior_pijs_alpha(i,j), Number_of_failed_infections(i,j)+prior_pijs_beta(i,j));
    end
end

[p_draws, n_draws, alpha_p, alpha_n, R0d_draws, R0c_draws] = MHsampler_Epi_noscale(m, sigma_n, n_0, p_0, Number_of_successful_infections, Number_of_failed_infections, Two_level', prior_pijs_alpha, prior_pijs_beta, pi_sim, qi_sim);

%%
n_postburnin = n_draws(ceil(0.1*m):end);
[f,xi]=ksdensity(n_postburnin, 'Support', [0,1]);
figure(200)
plot(xi,f, 'k', 'LineWidth', 1.5)
xlabel('$n$','Interpreter','latex')
set(gca,'FontSize',14)

%Mean_acc_ratio_n = table(mean(alpha_n,1));
%Mean_acc_ratio_q = table(mean(alpha_q,1));
R0d_mean_mode_var = zeros(3,N);
R0c_mean_mode_var = zeros(3,N);
R0d_mean_mode_var(1,:) = nanmean(R0d_draws);
R0d_mean_mode_var(3,:) = sqrt(nanvar(R0d_draws));
R0c_mean_mode_var(1,:) = nanmean(R0c_draws);
R0c_mean_mode_var(3,:) = sqrt(nanvar(R0c_draws));
n_mean_mode_var = zeros(3,1);
n_mean_mode_var(1,:) = nanmean(n_draws);
n_mean_mode_var(3,:) = sqrt(nanvar(n_draws));
x=0:0.0001:1;

pd = fitdist(n_draws,'Kernel','Kernel','epanechnikov','Width',0.1);
[y,ind]=max(pdf(pd,x));
n_mean_mode_var(2,:) = x(ind);

pijs_mode_alt = zeros(N,N);
x=0:0.0001:1;
pijs_se = zeros(N,N);

save('withpriors_noscale_bothdata.mat')

%%
load('withpriors_noscale_bothdata.mat');
load('mle_bothdata.mat');
%
count = 1;
X = 0.000:0.000001:0.999999;
figure(23);
prior_weight = 1;

for i=1:(N+1)
    for j=1:(N+1)
        if j==1&&i>1
            subplot(N+1,N+1,count);
            text(-0.5, 0.5, Sectors_semi_ab(i-1), 'FontSize',14);
            set(gca,'visible','off')
            set(gca,'xtick',[])
            set(gca,'ytick',[])
        elseif i==1&&j>1
            subplot(N+1,N+1,count);
            text(0, 0.0, Sectors_semi_ab(j-1), 'FontSize',14);
            set(gca,'visible','off')
            set(gca,'xtick',[])
            set(gca,'ytick',[])
        elseif i~=j&&i>1&&j>1
            subplot(N+1,N+1,count);
            p_draws_ij= p_draws(ceil(0.1*m):end,i-1,j-1);
            p_draws_ij_mle= p_draws_mle(ceil(0.1*m):end,i-1,j-1);

            [f,xi]=ksdensity(p_draws_ij(p_draws_ij<=1), 'Kernel','epanechnikov', 'BandWidth',0.3, 'Support', [0,1]);
            plot(xi,f, 'k', 'LineWidth', 1.5)
            hold on;
            y1 = betapdf(X, prior_pijs_alpha(i-1,j-1),  prior_pijs_beta(i-1,j-1));
            y2 = betapdf(X,  prior_pijs_alpha(i-1,j-1)+Number_of_successful_infections(i-1,j-1),  prior_pijs_beta(i-1,j-1)+Number_of_failed_infections(i-1,j-1));
            plot(X,y1,  '--r', 'LineWidth', 1)
            hold on;
            plot(X,y2,  ':b', 'LineWidth', 1)
            %hold on;
            %[f,xi]=ksdensity(p_draws_ij_mle(p_draws_ij_mle<=1), 'Kernel','epanechnikov', 'BandWidth',0.3, 'Support', [0,1]);
            %plot(xi,f, 'g', 'LineWidth', 1.5)
            %title(Sectors_semi_ab(i) +" to "+Sectors_semi_ab(j));
            if i==N+1 & j==N
                legend({'Posterior','Prior','Proposal distr.'});%, 'MLE'});
            end
            %set(gca,'xtick',[])
            set(gca,'fontsize', 12);
            ax = gca;
            outerpos = ax.OuterPosition;
            ti = ax.TightInset; 
            left = outerpos(1) + ti(1);
            bottom = outerpos(2) + ti(2);
            ax_width = outerpos(3) - ti(1) - ti(3);
            ax_height = outerpos(4) - ti(2) - ti(4);
            ax.Position = [left bottom ax_width ax_height];
            hold off;
            
            
            
            pd = fitdist(p_draws_ij,'Kernel','Kernel','epanechnikov','Width',0.0001);
            [y,ind]=max(pdf(pd,x));
            pijs_mode_alt(i,j) = x(ind);
            pijs_se(i,j) = sqrt(var(p_draws_ij));
        else
            
        end
        count=count+1;
          

    end
end




r0x=0:0.001:15;
for i=1:N
    pd = fitdist(R0d_draws(:,i),'Kernel','Kernel','epanechnikov','Width',0.1);
    [y,ind]=max(pdf(pd,r0x));
    R0d_mean_mode_var(2,i) = r0x(ind);
    pd = fitdist(R0c_draws(:,i),'Kernel','Kernel','epanechnikov','Width',0.1);
    [y,ind]=max(pdf(pd,r0x));
    R0c_mean_mode_var(2,i) = r0x(ind);
end
%%
figure(12)
count=1;

for i=1:N
    subplot(N,1,count);
    y1 = betapdf(X, 1, 1);
    y2 = betapdf(X, recovered_of_downturn(i)+1, number_of_downturns(i)-recovered_of_downturn(i)+1);
    plot(X,y1,'Color','r', 'LineWidth', 1.5)
    hold on;
    plot(X,y2,'Color','b','LineWidth', 1.5)
    title(Sectors(i));
    hold off;
    count=count+1;
end

figure(13);
count=1;
for i=1:N
    subplot(N,1,count);
    y3 = betapdf(X, numberOfCrisesRecov(i)+1, numberOfCrisesPeriods(i) - numberOfCrisesRecov(i)+1);
    plot(X,y1,'Color', 'r', 'LineWidth', 1.5)
    hold on;
    plot(X,y3, 'Color', 'b', 'LineWidth', 1.5)
    title(Sectors(i));
    hold off;
    count=count+1;
end

save('withpriors_noscale_bothdata.mat')


%% 

load('withpriors_noscale_bothdata.mat');
figure(18);
i = 5;
j = 4;
subplot(3,1,1);
p_draws_ij= p_draws(ceil(0.1*m):end,i,j);
[f,xi]=ksdensity(p_draws_ij(p_draws_ij<=1), 'Kernel','epanechnikov', 'BandWidth',0.3, 'Support', [0,1]);
plot(xi,f, 'k', 'LineWidth', 1.5)
hold on;
y1 = betapdf(X, prior_pijs_alpha(i,j),  prior_pijs_beta(i-1,j-1));
y2 = betapdf(X,  prior_pijs_alpha(i,j)+Number_of_successful_infections(i,j),  prior_pijs_beta(i,j)+Number_of_failed_infections(i,j));
plot(X,y1,  '--r', 'LineWidth', 1)
hold on;
plot(X,y2,  ':b', 'LineWidth', 1)
%title(Sectors_semi_ab(i) +" to "+Sectors_semi_ab(j));
%set(gca,'xtick',[])
set(gca,'fontsize', 20);
hold off;

subplot(3,1,2);
i = 5;
j = 7;
p_draws_ij= p_draws(ceil(0.1*m):end,i,j);
[f,xi]=ksdensity(p_draws_ij(p_draws_ij<=1), 'Kernel','epanechnikov', 'BandWidth',0.3, 'Support', [0,1]);
plot(xi,f, 'k', 'LineWidth', 1.5)
hold on;
y1 = betapdf(X, prior_pijs_alpha(i,j),  prior_pijs_beta(i-1,j-1));
y2 = betapdf(X,  prior_pijs_alpha(i,j)+Number_of_successful_infections(i,j),  prior_pijs_beta(i,j)+Number_of_failed_infections(i,j));
plot(X,y1,  '--r', 'LineWidth', 1)
hold on;
plot(X,y2,  ':b', 'LineWidth', 1)
%title(Sectors_semi_ab(i) +" to "+Sectors_semi_ab(j));
%set(gca,'xtick',[])
set(gca,'fontsize', 20);

hold off;


subplot(3,1,3);
i = 4;
j = 3;
p_draws_ij= p_draws(ceil(0.1*m):end,i,j);
[f,xi]=ksdensity(p_draws_ij(p_draws_ij<=1), 'Kernel','epanechnikov', 'BandWidth',0.3, 'Support', [0,1]);
plot(xi,f, 'k', 'LineWidth', 1.5)
hold on;
y1 = betapdf(X, prior_pijs_alpha(i,j),  prior_pijs_beta(i-1,j-1));
y2 = betapdf(X,  prior_pijs_alpha(i,j)+Number_of_successful_infections(i,j),  prior_pijs_beta(i,j)+Number_of_failed_infections(i,j));
plot(X,y1,  '--r', 'LineWidth', 1)
hold on;
plot(X,y2,  ':b', 'LineWidth', 1)
%title(Sectors_semi_ab(i) +" to "+Sectors_semi_ab(j));
%set(gca,'xtick',[])
set(gca,'fontsize', 20);

hold off;


save('withpriors_noscale_bothdata.mat')

%%

