Solution sketch, Numerical Difference Equations

Assignment

1. SIRS Model

Consider the SIRS model above and now assume that the $\beta_t$ is not constant, but instead assume that the $\beta_t$ is a function of the infections, so that when the number of infected goes above 5 \% of the population, social distancing is put in place and the $\beta$ value goes from .2 to .1. Once infections go below 2%, social distancing is relaxed again, and beta increases to .2. Can you change the SIR function to reflect these rules? How does this effect the the path of infections.

2. Seasonality in electricity consumption.

You want to simulate the dynamics of electricity consumption of a country. One aspect of this is including seasonality - that is regular movements over hours, days, weeks, months and years (and potentially more). In particular, if we have daily observations of electricity, we should see weekly seasonality. One way we can model this is through an AR(7) term (monday's are related to monday's, tuesday's to tuesday's, and so forth). Simulate a model with an AR(1) term and an AR(7) term where we have a stochastic process of normal iid's. You can consider the data to be normalized and scaled so that the initial value is zero and the stochastic error terms are distributed as Normal(0,1).

3. Setting up a dynamic economic/business problem as a difference or differential equation that you model numerically. (Part of obligatory Math Project)

You should explain each step of the set-up and computation. It can be helpful to create several scenarios, with various specifications of your model to show sensitivity to assumptions.You should interpret and critically evaluate the result, pointing out weaknesses and ways to improve the model.

solution

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
import numpy as np
import pandas as pd

#1 

def SIR(N, S0, I0, beta, gamma, t):  
    S=[S0]
    I=[I0]
    
    for i in range(1,t):
        perc_infected = I[i-1]/N
        if perc_infected>.05:
            beta = .1
        elif perc_infected<.02:
            beta = .2
        
        I.append(I[i-1] + beta*S[i-1]/N*I[i-1] - gamma * I[i-1]) #1
        S.append(S[i-1] - beta*S[i-1]/N*I[i-1]) #2
        if I[i]<1 or S[i]<1: #3
            break
            
    S = np.array(S)
    I = np.array(I)
    R= N-S-I
    new_t = len(S)
    time = np.arange(new_t)
    
    SIRS_df = pd.DataFrame({"time":time, "S":S, "I":I, "R":R})
    SIRS_df = SIRS_df.loc[(SIRS_df.S>0) & (SIRS_df.I>0),:].copy() #4
    SIRS_df["AR"] = (SIRS_df["I"] + SIRS_df["R"])/N #5
    nr = SIRS_df.shape[0]
    rate = np.array(SIRS_df["I"][1:nr])/np.array(SIRS_df["I"][:(nr-1)])
    rate = rate[0] + rate
    SIRS_df["rate"] =  np.insert(rate, 0, rate[0])
    return(SIRS_df)

SIR_mod1 = SIR(N=1e6, S0=1e6-1, I0=1, beta=0.2, gamma=1/18, t=300)
SIR_mod1

	time	S	I	R	AR	rate
0	0	999999.000000	1.000000	0.000000	0.000001	2.288888
1	1	999998.800000	1.144444	0.055556	0.000001	2.288888
2	2	999998.571112	1.309753	0.119136	0.000001	2.288888
3	3	999998.309161	1.498939	0.191900	0.000002	2.288888
4	4	999998.009374	1.715452	0.275174	0.000002	2.288888
...	...	...	...	...	...	...
295	295	211378.689792	10898.298907	777723.011301	0.788621	2.131258
296	296	210917.956164	10753.571485	778328.472351	0.789082	2.131164
297	297	210464.331900	10609.775111	778925.892989	0.789536	2.131072
298	298	210017.736054	10466.939007	779515.324940	0.789982	2.130982
299	299	209578.087487	10325.090962	780096.821551	0.790422	2.130892

300 rows × 6 columns

SIR_mod1.I.plot()

<AxesSubplot:>

SIR_mod1.rate.plot()

<AxesSubplot:>

2

T=365
time = np.arange(T)
x = np.empty(T)
a1=.5
a7=.3
x[0] = np.random.normal(0,1) #mean 
for t in range(1,5):
    x[t] = a1*x[t-1] + np.random.normal(0,1)

x[5] = a1*x[4] + np.random.normal(0,1) 
x[6] = a1*x[5] + np.random.normal(0,1)

1
2
3
4

for t in range(7,T):
    x[t] = a1*x[t-1] + a7*x[t-7] + np.random.normal(0,1)

fig, ax = plt.subplots()
ax.plot(time, x)

[<matplotlib.lines.Line2D at 0x7fb0e8806df0>]

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

plot_acf(x)

plot_pacf(x)