Financial Contagion (Epidemic spread) Model Meets Problem
Recently I want to rebuild the financial contagion model in the paper: Contagion in Financial Networks by Prasanna Gai. Now I am stuck in the first figure:
(figure 3 actually).
What I Did
I used Python and networkx.
First, build the ER network with 1000 nodes, and the probability depends on which average degree I want to simulate. For example, if I want to simulate the average degree at 3, then the probability of generating ER network is 3/(1000-1) where 1000 is the network size.
Then for each node, I find how many nodes point towards it and count, to calculate the AiIB (weights). If node 1 has 3 nodes point towards it, then the weights on these edges are AiIB (0.2 in paper) / 3 (number of neighbours).
To simulate contagion, first randomly chose one node to remove all its assets. then it cannot pay back the liabilities to its neighbours if the liability is more than the capital buffer (Ki, 0.04 in the paper). For those banks received liabilities from more than one banks, even each link's weight is less than the Ki, if the sum of these liabilities is more than Ki it still is considered as bankrupt. The model is just like epidemic spread where the new bankrupt banks will influence the new batch of banks, ended at no more banks are bankrupted in this system.
As contagion is defined as more than 5% of the banks bankrupt in this network ( 50 in this case).
To plot the figure, each average degree needs to test 100 times:
probability = number of contagion happened / simulation times 100 here
extent = [under the case of happening contagion] sum of the fraction of bankrupt banks/number of contagion happened.
The raw code is available on GitHub. By running er_100.py you can get the figure of mine like this:
If you have any trouble of the code please let me know. (the code needs at least 1 hours to run on the GCP with 8 vCPU...)
I also tried the network with 60 nodes here is it looks like: 
it has a little similar shape with figure 1. But this still not good and the small network is not what I want.
I don't know what problems with my code. In my view, I have considered everywhere and it should get similar results. I even started to question the authority of the paper...
If you have any ideas please help me with this.
This is a "hard" question to answer.
I still do not find any clue to the code. Then I rewrite the code in R and run it, here is the draft figure I got:
As you can see now the figure is just the one in the paper. But the algorithm and the structure were totally the same with I wrote it in Python.
Maybe this is a case that shows Python just can't do. If anyone feels interesting of this question, want to divide more among the difference between Python and R, this is a great example. And I am very happy to offer any helps.
By the way the model code in R is available in GitHub and it still under update.
For those who spent time on reading my description, Thanks for your time.
update:
I am also can't believe this because in my view the code is doing count and calculate that's simple. I print out very thing at every step and check every node, from 10 nodes small network to 1000 nodes network, the log file reached over 50G. All looks normal and the number (bankrupt one) just not reach the threshold. Not like in R, with totally the same structure, the results is just same with the paper.
I really don't know why and have no idea.