This term paper is mandatory for this module. Please see the course outline and postgraduate handbook for details. An electronic copy must be submitted via Faser no later than 12pm (noon) on the date of the submission deadline.

You are required answer the questions, and also print your code and any graphs that you are asked to save and attach them to your submission.

You must answer ALL of the following questions. There are a total of 100 marks.

1    Term paper problem:  Simulating a pandemic

In this exercise you will code the basic infectious disease model, known as the “SIR Model”. The letters stand for Susceptible, Infected, Recovered. This model predicts how many people will become infected over time during a pandemic. The model is very simple, and very easy to code. The model can be used to provide a basic description of how the COVID pandemic would spread through a country without any government intervention. We will then extend the model to add more economic features.

In the text below I explain to you how the model works, and then set the questions you must complete. You will write code to simulate and plot the course of the epidemic under various scenarios. After we finish our Python programming course, you will be provided with half-finished codes which you will complete in order to answer the questions. You may alternatively code your answers in another language if you desire.

1.1    The basic SIR model

Consider a country with an initial population of N people at the beginning of the pandemic.   Let t = 0, 1, 2,...,T index time.  We will take a time period to be one week, and define t = 0 as the first week of the pandemic, t = 1 as the second week, and so on. At time t = 0 suppose that one person is exogenously infected with the virus (e.g.  we could imagine they just came back from holiday).  From time t = 1 onwards, the virus will spread throughout the population. We stop the simulation at period T.

People can be in one of four states:

1. Let It  denote the total number of people currently Infected with the virus at date t.

2. Let Rt denote the total number of people who have had the virus in the past and now Recovered.

3. Let St denote the total number of people who have not had the virus yet, and so are Susceptible to catching it.

4. Let Dt  denote the total number of people who have Died from the virus up to date t.

Since an individual person can only be one of these things at a given time, the total number of people in the four states must always add up to the total population:

St + It + Rt + Dt = N                                                          (1)

We make some simplifying assumptions:  1) We suppose that immunity from the virus last forever, so that once you have had the virus and are Recovered, you never become Susceptible again, 2) We suppose that there are no births of new children in the model, or deaths from other causes, so the population is always N. In fact, since Dt people have died from the disease, the total population alive at any time t is N − Dt .

At time 0, when the pandemic starts, we have the following values of the four states. Let’s suppose that the population of the country is one million people:  N = 1,000,000.  At time 0 we have I0  = 1 person who is initially infected with the disease. At time 0, no one has died or recovered from the disease yet, so we have R0  = D0  = 0.  This leaves the rest of the non-infected population to be Susceptible: S0 = N − I0 = 1,000,000 − 1 = 999,999. In summary:

S0 = 999, 999         I0 = 1         R0 = 0         D0 = 0                                   (2)

Over time, the one infected person will infect more people, and It  will grow.  The pandemic evolves because of the natural forces of people becoming infected,  recovering,  and potentially dying.   For simplicity, we assume the law of large numbers holds to remove randomness from the model, so that we know exactly how many people will transition between states each period. We assume:

1. A fraction δ of infected people die every period (δIt transition from It to Dt+1)

2. A fraction γ of infected people recover every period (γIt transition from It to Rt+1)

3. If there are St  susceptible people and It  infected people this period, this will lead to αN(St)It   new infections next period. We call α the contact rate. (αN(St)It   transition from St to It+1)

The forces δ and γ are intuitive, but the process of new infections is a little more mysterious: why is it αN(St)It ? This is derived from a random matching assumption, but I will not discuss it further here. All you need to know is that α controls how strongly the virus spreads: if α is large, then infected people infect many more people, and if α is small then infected people do not infect many people. Think of α as capturing two ideas:  1) the more social contact people have on average (e.g.  going to restaurants, working, seeing friends) the more chances they have to pass on the virus, and the higher is α, 2) the more infectious the virus is, the more likely it is to be passed on to a new person whenever people meet, and the higher is α .

Given the forces above, the following four equations describe the evolution of the pandemic:

It+1 = (1 − γ − δ)It +                                                               (3)

Rt+1 = Rt+ γIt                                                                                      (4)

Dt+1 = Dt + δIt                                                                                      (5)

St+1 = N − It+1 − Rt+1 − Dt+1                                                                      (6)

Equations (3) to (6) are the whole model. It is simple, as promised! The intuition for each equation is as follows.  In (3) we see that the number of infected people at t + 1 is simply the remaining infected people at t who did not recover or die (the fraction 1 − γ − δ) plus the new infected people αStIt/N. In (4), the new total number of recovered people, Rt+1, is the old number of recovered, Rt, plus the people who just recovered this period, γIt. Similarly, in (5), the new total number of dead people, Dt+1 is the old number, Dt, plus the new deaths δIt. Finally, in (6) we can compute the number of people who are still susceptible, St+1 as simple the total population N, less the people in the other three states (this is just equation (1) rearranged).1

The model is very easy to use. If we know the number of people in each state this week (St , It , Rt , and Dt) we can use equations (3) to (6) to compute the number of people in each state next week (St+1 , It+1 , Rt+1, and Dt+1).  This means that we can start with the number in each state at time t = 0 as given in (2), and simulate the model forward by next computing the solution at time t = 1, then t = 2, and so on up to period T.

1.2    The R number (a.k.a. the Reproduction number)

Now, let’s define the famous “R number” that is always discussed on the news. The R number is short for the “Reproduction number”.  At any time t, it is defined as the number of new people that each infected person will infect on average.  Since we already use Rt  to denote the number of recovered people, I will use a scripted letter, Rt, to denote the R number. The definition is:

Rt =                                                                               (7)

The derivation is given in this footnote, for those interested.2  You will also sometimes hear about the “R0 number” of a virus. This is just the R number, Rt, taken in the first period when there is only one person infected. That is, the R0 is just R0:

R0 =                                                                               (8 

1.3    Algebra questions (30 marks)

This section contains algebraic questions to introduce you to thinking about the model. You must show your algebraic steps in all answers.

1.  (10 marks) Combine equations (3) and (7) to show that

= 1 + (γ + δ)(Rt − 1) It

Hence, argue that if the R number is greater than one (Rt  > ( > 1), and vice versa.

(9)

1), then the pandemic is growing

2. Recall that N = 1,000,000 and S0 = 999,999. On the first day of the pandemic (t = 0), (8) implies

that the R0 number is given by

R0 =                                                                              (10)

We roughly calibrate the model to the COVID virus as follows. Each period of the model is one day.  People are infectious for around one to two weeks, so we set γ = 1/10 = 0.1, so that the average length of infection is ten days. The death rate from COVID appears to be around 1%, so we set δ = 0.01/10 = 0.001.3

(a)  (5 marks) Using (8), and given these values of S0 , N , γ, and δ, show that R0 > 1 as long as α >  × (0.1 + 0.001) ≃ 0.101.

(b)  (5 marks) Initial estimates of the R0 number for COVID suggested a number around 3. Using (8), and given our values of S0 , N , γ, and δ, show that the value of α which implies that R0 = 3 is α = 3 ×  × (0.1 + 0.001) ≃ 0.303.

3.  (10 marks) Suppose that a government policy (such as a test and trace system, or lockdown) reduces the spread of the virus, by reducing α to α = 0.05.  Using (8), calculate the R0  of the virus with this new policy in place. Does the policy increase or decrease the R0?

1.4    Coding exercise 1: the basic SIR model (50 marks)

You will now be asked to simulate the model on the computer. You could code this from scratch yourself if desired.  However, half-completed codes in Python will be made available for you to complete after we finish our coding classes. The half-finished codes are now available on the EC901 Moodle under the name “Coding exercise 1: codes for you to complete” .

For the first coding exercise, consider the basic model from equations (3) to (6).   Assume that α = 0.303, and continue to assume that γ = 0.1 and δ = 0.001.  Recall that S0  = 999,999, I0  = 1, R0 = 0, D0 = 0, and N = 1,000,000. Set T = 250, and hence simulate the pandemic for 250 days.

Hint: You will notice that the pandemic “ends” much before T = 250, since there is no policy in this model, and the population reaches herd immunity.  However, if coded correctly, this implies that over 90% of the population become infected with the virus.

4.  (40 marks) Simulate the path of the pandemic from t = 0 to t = T given the parameter values

above. Plot the values of St , It , Rt, and Dt  over time in well formatted graphs. Save your code and graphs to submit with your answers.

Hint 1: To simulate St , It , Rt, and Dt  over time, you should use a “for” loop.     Hint 2: To make graphs in Python, you will need to use the Matplotlib package.

5. Using your simulated results, answer:

(a)  (2 marks) What fraction of people will have contracted the virus at some point?

Hint: This is everyone who is not still Susceptible at time T, i.e. (N − S(T))/N . (b)  (2 marks) What fraction of people died from the virus?

Hint: This is everyone who is Dead at time T, i.e. D(T)/N .

(c)  (1 mark) At roughly what date is the number of infections, I(t), highest?

6.  (5 marks) Repeat question 4 but now set α = 0.05, saving your answers as before. Does a large pandemic occur in this case? Discuss, making sure to compare the R0  when α = 0.05 to the R0 when α = 0.303.

1.5    Coding exercise 2:  SIR model with optimising agents (20 marks)

For the next questions, we extend the model to incorporate insights from economics.  Specifically, we now allow the contact rate, α, to vary over time. That is, we replace α with a time-varying parameter αt . We suppose that αt  might change over time in response to changes in people’s behaviour. Specifically, the basic SIR model assumes that the contact rate is constant, but this is a strange assumption.  If people are scared of catching the virus, people might voluntarily choose to have less social contact when the virus becomes more prevalent, and they become more worried about catching it.  The idea that people optimise is one of the central lessons of economics, and hence somewhere where economics might be able to help with disease modelling. I provide a simple example in the questions below.

Specifically, let’s suppose that the higher the current number of infections (higher It) the more scared people are of the pandemic, and the more likely they are to stay home or take other safety measures which reduce αt. We represent this via a reduced form behavioural equation:

αt =  (It)−β                                                                                       (11)

Here, β controls the strength of the behavioural response: the larger is β, the more people respond to high infections (high It) by staying home and reducing αt. If β = 0, then the equation becomes αt =  , which is just our basic SIR model from the last section.   controls the average level of the contact

rate. With a time-varying αt, equation (3) is replaced with It+1  = (1 − γ − δ)It + αtN(St)It . Plugging in αt =  (It)−β, this becomes It+1 = (1 − γ − δ)It +  . We use this equation to replace (3) in our

basic SIR model. With this change, the core equations of the model now become:

It+1 = (1 − γ − δ)It +                                                                  (12)

Rt+1 = Rt+ γIt                                                                                    (13)

Dt+1 = Dt + δIt                                                                                    (14)

St+1 = N − It+1 − Rt+1 − Dt+1                                                                   (15)

The model is now described by equations (12) to (15), but otherwise functions identically to the basic model from the earlier questions.