In [ ]:
fig, ax = plt.subplots()
ax.plot(t,k_array)
Out[ ]:
[<matplotlib.lines.Line2D at 0x7f7ef91ff370>]
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5) #set default figure size
import numpy as np
In your course(s) in macroeconomics, you should be familiar with the Solow growth model. A feature of the solow growth model is that the dynamics of capital growth determines (short term) economic growth. We can write the growth of capital, k, with the following difference equation:
$$k_{t+1} = sAk_t^{\alpha} + (1-\delta)k_t$$The above is derived from a aggregate production function:
$$y = Ak^{\alpha}$$The difference equation then describes the process of capital change.
On the positive side, you add capital by investing, which in macroeconomics is the same as saving. So we add capital by the share of saving,s, times the total productivity (you can think: which share of production do we not consume).
On the negative side, every year, a certain share, $\delta$ of capital disappears in wear-and-tear (depreciation).
We can assign some numbers to our parameters. Here we have a 30% savings rate, a 2 as total factor productivity and 40% depreciation. These numbers are absurd.
A, s, alpha, delta = 2, 0.3, 0.3, 0.4
We can then define our function that will iterate the capital stock, k, forward by one period
def g(k):
return s*A * k**alpha + (1 - delta) * k
We can then iterate numerically setting the initial capital stock to .1
k_array = [.1]
for i in range(19):
k_array.append(g(k_array[i]))
k_array=np.array(k_array)
t = np.arange(0,20,1)
k_array
array([0.1 , 0.36071234, 0.65830253, 0.9242535 , 1.14053991,
1.30846737, 1.4354798 , 1.53001722, 1.5996583 , 1.65060626,
1.6877033 , 1.71462678, 1.73412194, 1.74821545, 1.75839224,
1.76573479, 1.77102933, 1.77484549, 1.77759524, 1.77957616])
fig, ax = plt.subplots()
ax.plot(t,k_array)
[<matplotlib.lines.Line2D at 0x7f7ef91ff370>]
Notice that the capital stock approaches a constant point - this is what we call a steady-state, where the system is in dynamic equilibrium: the amount being added in capital (through saving) is exactly equal to what is being lost to depreciation.
Using our production function, we can also see the time trend of production
y_array = A*k_array**alpha
fig, ax = plt.subplots(2)
ax[0].plot(t,k_array, label="k")
ax[1].plot(t,y_array, color="green", label="y")
fig.legend()
<matplotlib.legend.Legend at 0x7f7ef92828e0>
So the implication of this model is that economic growth also moves towards a steady-state of no economic growth
If we start with an initial condition of k=3, how does this change the dynamics of k and y. Can you explain why?
We can use some functions for drawing typical graphical figures from quantecon - these should be familiar from your macroeconomic classes. You don't need to spend much time understanding exactly what they do, instead we can just use them to try to understand the underlying dynamics.
def subplots(fs):
"Custom subplots with axes throught the origin"
fig, ax = plt.subplots(figsize=fs)
# Set the axes through the origin
for spine in ['left', 'bottom']:
ax.spines[spine].set_position('zero')
ax.spines[spine].set_color('green')
for spine in ['right', 'top']:
ax.spines[spine].set_color('none')
return fig, ax
def plot45(g, xmin, xmax, x0, num_arrows=6, var='x'):
xgrid = np.linspace(xmin, xmax, 200)
fig, ax = subplots((6.5, 6))
ax.set_xlim(xmin, xmax)
ax.set_ylim(xmin, xmax)
hw = (xmax - xmin) * 0.01
hl = 2 * hw
arrow_args = dict(fc="k", ec="k", head_width=hw,
length_includes_head=True, lw=1,
alpha=0.6, head_length=hl)
ax.plot(xgrid, g(xgrid), 'b-', lw=2, alpha=0.6, label='g')
ax.plot(xgrid, xgrid, 'k-', lw=1, alpha=0.7, label='45')
x = x0
xticks = [xmin]
xtick_labels = [xmin]
for i in range(num_arrows):
if i == 0:
ax.arrow(x, 0.0, 0.0, g(x), **arrow_args) # x, y, dx, dy
else:
ax.arrow(x, x, 0.0, g(x) - x, **arrow_args)
ax.plot((x, x), (0, x), 'k', ls='dotted')
ax.arrow(x, g(x), g(x) - x, 0, **arrow_args)
xticks.append(x)
xtick_labels.append(r'${}_{}$'.format(var, str(i)))
x = g(x)
xticks.append(x)
xtick_labels.append(r'${}_{}$'.format(var, str(i+1)))
ax.plot((x, x), (0, x), 'k-', ls='dotted')
xticks.append(xmax)
xtick_labels.append(xmax)
ax.set_xticks(xticks)
ax.set_yticks(xticks)
ax.set_xticklabels(xtick_labels)
ax.set_yticklabels(xtick_labels)
bbox = (0., 1.04, 1., .104)
legend_args = {'bbox_to_anchor': bbox, 'loc': 'upper right'}
ax.legend(ncol=2, frameon=False, **legend_args, fontsize=14)
plt.show()
def ts_plot(g, xmin, xmax, x0, ts_length=6, var='x'):
fig, ax = subplots((7, 5.5))
ax.set_ylim(xmin, xmax)
ax.set_xlabel(r'$t$', fontsize=14)
ax.set_ylabel(r'${}_t$'.format(var), fontsize=14)
x = np.empty(ts_length)
x[0] = x0
for t in range(ts_length-1):
x[t+1] = g(x[t])
ax.plot(range(ts_length),
x,
'bo-',
alpha=0.6,
lw=2,
label=r'${}_t$'.format(var))
ax.legend(loc='best', fontsize=14)
ax.set_xticks(range(ts_length))
plt.show()
Below we'll create a chart showing kapital, $k_t$, on the x-axis and the difference equation function, $g(k_t) = k_{t+1}$, on the y-axis, with a 45 degree line indicating the points where $k_t = k_{t+1}$.
xmin, xmax = 0, 4 # Suitable plotting region.
plot45(g, xmin, xmax, 0, num_arrows=0)
Using this chart, we can show the dynamics in the difference equation. Starting out with an initially low level of initial capital, we quickly add new amounts of capital.
However, as the amount of existing capital grows, then the absolute amount of capital we lose each year to depreciation grows. We gradually approach a steady-state, which is shown as the place where the blue line crosses the 45 degree black-line. Here, the amount of capital we add each year from new saving/investment exactly matches how much capital we lose through depreciation.
k0 = 0.25
plot45(g, xmin, xmax, k0, num_arrows=5, var='k')
/var/folders/q_/twndxknx06n5v1nllkjtz0km0000gn/T/ipykernel_69475/2261905364.py:50: UserWarning: linestyle is redundantly defined by the 'linestyle' keyword argument and the fmt string "k-" (-> linestyle='-'). The keyword argument will take precedence. ax.plot((x, x), (0, x), 'k-', ls='dotted')
What point represents the steady-state in the above chart.
Show the above plots when the initial condition for capital is k=3
The above model was non-linear, yet the dynamics were still quite orderly and what we call regular, but we can also have more chaotic-seeming, or irregular dynamics. An example is what we call a quadratic map, where the iteration function is:
$$x_{t+1} = 4*x_t*(1-x_t)$$Where x is between 0 and 1. What type of dynamics can this describe?
xmin, xmax = 0, 1
g2 = lambda x: 4 * x * (1 - x)
x_array = [0.3]
for i in range(19):
x_array.append(g2(x_array[i]))
x_array=np.array(x_array)
t = np.arange(0,20,1)
fig, ax = plt.subplots()
ax.plot(t,x_array)
[<matplotlib.lines.Line2D at 0x7ff44153ec40>]
This almost looks random, but it is actually completly deterministic.
We can plot the trajectory using the quantecon 45\% line
plot45(g, xmin, xmax, x0, num_arrows=6)
/var/folders/fs/kr83xdyx2jq_8gww_xfmtxfr0000gn/T/ipykernel_1632/2261905364.py:50: UserWarning: linestyle is redundantly defined by the 'linestyle' keyword argument and the fmt string "k-" (-> linestyle='-'). The keyword argument will take precedence. ax.plot((x, x), (0, x), 'k-', ls='dotted')
Even with just one variable, you can generate some pretty complex dynamics!
During the Covid 19 pandemic you may have seen presentations of models of the spread of the virus under different scenarios. These types of models are generally referred to as SIRS models.
A description of a SIR model (with economic implications) can be found in this working paper here, which I base my description here on.
The following code is based on this post on the website r-bloggers.
In this model you divide the population (N) into 4 categories (N is usually set to one, so you are looking at proportion of population, where the categories sum to 1).
S: Number susceptible - that is, those that have not been previously infected and have no immunity
E: The number who are exposed
I: The number who are infected
R: The number who recover or die (and are no longer susceptible)
The model is based on 5 relatively simple difference (or differential) equations which describe in a simple manner the spread of the disease.
We'll start with the simplest relationship - the change in the number of recovered/dead is a function of the number of infected:
$$dR/dt = \gamma I$$As a difference equation we can write:
$$R_t - R_{t-1} = \gamma I_t$$$$R_t = R_{t-1} + \gamma I_t$$Here $\gamma$ is the rate at which people either recover/die from the disease once they have it. So if it takes on average 18 days from infection to get over the disease, then $\gamma=1/18$.
The rate of infection is determined by the number of people exposed to the disease and the rate of recovery:
$$dI/dt = \sigma E - \gamma I$$Again, as a difference equation:
$$I_t - I_{t-1} = \sigma E_t - \gamma I_t$$So here $\sigma$ is interpreted as the rate which people who are exposed to the disease become infected. This is the inverse of the incubation period--how long from exposure to infection. So if that is 5.2 days, then $\sigma=1/5.2$.
The rate of recovery/death counts negatively towards the total rate of infection
The rate of exposure is in turn a function of the proportion of the population who are susceptible and the number who are infected:
$$dE/dt = \beta_t \frac{S}{N}I - \sigma E$$or
$$E_t - E_{t-1} = \beta_t \frac{S_t}{N_t}I_t - \sigma E_t $$$\beta_t$ is probably the most important parameter in terms of understanding the spread of the disease and the public-health measures put in place. This parameter represents how easily people "bump-in" to one another, and spread the virus. By putting in "social-distancing" measures, the government is trying to reduce the amoung of "bumping in". We'll see the effect of this in our simulation below.
We also have a rate $\rho_t$ which is the ratio of this "bump-in" rate, $\beta_t$, to the recovery/death rate, $\gamma$. This ratio has gotten a lot of attention when the media discusses these models.
$$\rho_t = \frac{\beta_t}{\gamma}$$(I spent a lot of time being confused about this, because in the paper I referenced, they write $\rho_t$ as $R_t$, which is completly distinct from the number of people who are recovered, R. This seems to be a gross violation of good notation.)
Finally, the rate of change of those who are susceptible is negatively proportional to the bump-in rate, the existing proportion who already are susceptible and the proportion who are infected already.
$$dS/dt = - \beta_t \frac{S}{N}I$$or
$$S_t - S_{t-1} = - \beta_t \frac{S_t}{N_t}I_t$$We can program in a simplified simulation of this model in R. Below we basically ignore the nuance between exposure and infection, ignoring the dynamics that come in the period between being exposed to the infection and showing symptoms. We also fix $\beta_t$
First, we create our main function, called SIR. As input, we enter the population N, and the initial number of susceptible (S0), infected (I0), our bump-in parameter (beta), and our recovery rate (gamma), and the number of periods t.
import pandas as pd
def SIR(N, S0, I0, beta, gamma, t):
S=[S0]
I=[I0]
for i in range(1,t):
S.append(S[i-1] - beta*S[i-1]/N*I[i-1]) #1
I.append(I[i-1] + beta*S[i-1]/N*I[i-1] - gamma * 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)
So, walking through, we first initialize our vectors of S and I and set the initial conditions
Then we begin a loop through the time periods.
For each period (following 1) the number susceptible is equal to the previous periods susceptible minus those that have newly become infected.
The number infected is equal to the previous number plus our rate of new infections, less those that have recovered
The model exits early if the the number of infected or susceptible goes below 1 (basically, the case where the virus is eliminated from the population, either because everyone has already recovered or died, or because it has been completly eliminated through public health measures.)
We collect the series into a dataframe and then limit the values to those where S and I are above 0.
We add together the number infected and recovered/died at all time periods - these are those that are no longer susceptible
Now we'll run our model with:
SIR_mod1 = SIR(N=1e6, S0=1e6-1, I0=1, beta=0.4, gamma=1/18, t=300)
SIR_mod1
| time | S | I | R | AR | rate | |
|---|---|---|---|---|---|---|
| 0 | 0 | 999999.000000 | 1.000000 | 0.000000 | 0.000001 | 2.688888 |
| 1 | 1 | 999998.600000 | 1.344444 | 0.055556 | 0.000001 | 2.688888 |
| 2 | 2 | 999998.062224 | 1.807530 | 0.130247 | 0.000002 | 2.688888 |
| 3 | 3 | 999997.339213 | 2.430122 | 0.230665 | 0.000003 | 2.688888 |
| 4 | 4 | 999996.367167 | 3.267161 | 0.365672 | 0.000004 | 2.688887 |
| ... | ... | ... | ... | ... | ... | ... |
| 286 | 286 | 405.151017 | 1.224354 | 999593.624629 | 0.999595 | 2.289051 |
| 287 | 287 | 405.150819 | 1.156532 | 999593.692649 | 0.999595 | 2.289051 |
| 288 | 288 | 405.150632 | 1.092468 | 999593.756900 | 0.999595 | 2.289051 |
| 289 | 289 | 405.150455 | 1.031952 | 999593.817593 | 0.999595 | 2.289051 |
| 290 | 290 | 405.150287 | 0.974789 | 999593.874924 | 0.999595 | 2.289051 |
291 rows × 6 columns
We can use the plotting function within pandas, first we can look at the number of infections:
SIR_mod1.I.plot()
<Axes: >
Above we get the typical bell shape of infection, which should seem familiar after going through a pandemic.
Here we get the rate of infection over time:
SIR_mod1.rate.plot()
<Axes: >
Now we move on to a set of dynamics where we also include a random process. One of the simplist such models is an autoregressive 1 (AR1) model. This type of modelling is seen often in time series statistical modelling. The general formulation we write for such a model is:
$$X_{t+1} = aX_t + b + cW_{t+1}$$where $cW_{t+1}$ describes the (random) process.
You may have seen a simpler form of this equation:
$$X_{t+1} = aX_t + b +\epsilon_{t+1}$$Here we have set $cW_{t+1} = \epsilon_{t+1}$
Where the $\epsilon_{t+1}$ represents a random "shock" from a normal distribution - though this is just one of many potential probability distributions and procesees you could formulate
Iterating our AR model backwards from a certain t back to our initial condition we could write our equation in what is called moving average form as:
$$X_t = a^t X_0 + b \sum_{j=0}^{t_1} a^j + c \sum_{j=0}^{t-1} a^j \epsilon_{t-j}$$What we can see from the above equation is that an AR system is defined in terms of the sum of its shocks.
Consider the system above at t=4. Iterate backwards to the initial condition $X_0$ to show the system in MA form. If we set a = .8, what will the equation look like?
Let's try simulating with some real numbers, we'll start with a simple AR(1) model with normal random (iid) process from a N(0,1) distribution. And we'll set a=.8 and the initial value at 0. We'll set b=0. So we have
$$x_t = .8 x_{t-1} + \epsilon_t$$T=50
time = np.arange(T)
x = np.empty(T)
a=.8
#set initial value to be random N(0,1)
x[0] = np.random.normal(0,1)
for t in range(1,50):
x[t] = a*x[t-1] + np.random.normal()
x
array([-1.19324639, -1.77607128, 0.44942104, -0.29886564, -0.58993938,
0.18200874, -0.49797719, -0.34901524, -1.88034084, -1.83860015,
-1.84196321, -1.0822786 , -2.10133894, -1.27150065, -2.72898305,
-2.27448231, -0.94931534, -0.47609661, -2.03949451, -1.25212438,
-1.9687932 , -2.03991867, -1.11808446, -2.60357297, -1.61888514,
-2.57823399, -2.35821763, -2.96125083, -1.93235013, -0.55234647,
-1.38111428, 1.45481469, 2.35188872, 1.06482383, 1.4716766 ,
-0.61038256, -2.17456238, -2.00650836, -1.27615799, -1.16879317,
-1.71741774, -2.39219604, -2.75252326, -3.01706648, -3.21048712,
-4.54777702, -3.7022227 , -2.31717874, -1.88690766, -0.01479741])
We can chart these
fig, ax = plt.subplots()
ax.plot(time, x)
[<matplotlib.lines.Line2D at 0x7f7ed8e2c2e0>]
Now let's do the same exercise 100 times and show all of the series at the same time.
T=50
S=100
time = np.arange(T)
X=np.empty([T,S])
#initial values for each simulation, s
X[0,:] = np.random.normal(0,1,S)
for t in range(1, T):
X[t,:] = a*X[t-1,:] + np.random.normal(0,1, S);
Now let's plot our series
fig, ax = plt.subplots()
for s in range(S):
ax.plot(time, X[:,s], color="blue", alpha=.1)
We can notice a few things about the above series:
After a shock, the series seems to come back to a certain average value. Another way of saying this is that we could take the average of different portions of this series and get about the same answer.
The variance seems also to be roughly constant throughout. If we look at the different sections of the series, and pay attention to the different simulation values, we find that roughly speaking the variance appears about constant.
What we have described can be called covariance stationarity which is a version of a more general idea of stationarity. Basically this means that if we measure the covariance between any two points in a series, then that covariance goes towards zero as we move those two points farther away from each other. In other words, shocks dissipate.
Quantecon describes stationarity a bit more broadly in terms of convergence towards a fixed probability distribution.
Now let's go on to see what a non-stationary distribution looks like.
Now we will set our a=1. Our equation will then look like:
$$x_t = x_{t-1} + \epsilon_t$$So, in words, the value today is equal to the value yesterday plus some random innovation, $\epsilon$. This is what we call a random walk
Backwards iterate the above function to t=0 from for example t=4. What does this tell you about how shocks integrate into the values of the series?
x = np.cumsum(np.random.normal(0,1,50))
time = np.arange(50)
If you tried the exercise above, you can see that any value at a certain time in a random walk is simply the sum of all the past random shocks, which is the way we created the above random walk series. Another way of saying this is that all the shocks to the system are permanent.
fig, ax = plt.subplots()
ax.plot(time, x)
[<matplotlib.lines.Line2D at 0x7f7ef97b4610>]
And now we'll generate many more series from our random walk
X=np.empty([T,S])
for s in range(S):
X[:,s] = np.cumsum(np.random.normal(0,1,T))
fig, ax = plt.subplots()
for s in range(S):
ax.plot(time,X[:,s], color="blue", alpha=.1)
For our random walk series, we see that it would be difficult to define a mean value for the process. Because shocks are permanent in the random walk, then the series has no tendency to go back to some well defined mean value.
We can also notice from running multiple series that the variance increases over time, so we do not have a well defined constant variance for a random walk.
A random walk is covariance non-stationary in that the shocks are permanent, and therefor the covariance between any two points in a series will not go towards zero no matter how far apart in time we move these points from each other.
Again, quantecon defines non-stationarity in terms of a series not converging to a fixed probability distribution.
A related concept to stationarity is the idea of Ergodicity. QuantEcon gives a mathematical definition, but here we will give some intuition through computation.
First consider our AR(1) process with a=.8, where we generate one time-series and then where we estimate our time series mean and standard deviation
T=5000
time = np.arange(T)
x = np.empty(T)
a=.8
#set initial value to be random N(0,1)
x[0] = np.random.normal(0,1)
for t in range(1,T):
x[t] = a*x[t-1] + np.random.normal()
ts_mean = np.mean(x)
ts_std = np.std(x)
print(ts_mean, ts_std)
-0.0759047801544835 1.6107623872318737
Now we will again generate many series from the same process and this time we will choose a certain time period (arbitrarily the mid-point of time) and then take the mean and standard deviation of the cross-section at that time period:
S=10000
midT = round(T/2)
X=np.empty([T,S])
X[0,:] = np.random.normal(0,1,S)
for t in range(1, T):
X[t,:] = a*X[t-1,:] + np.random.normal(0,1, S);
cs_mean = np.mean(X[midT,:])
cs_std = np.std(X[midT,:])
print(cs_mean, cs_std)
-0.030307788398191315 1.6810131648159379
With time series length of 50 and 100 simulations, chances are that you will get means and standard deviations that are quite different from each other.
But in the above, try putting in T=500 and S=1000 and then T=5000 and S=10000.
What you should notice is that the mean and variance (which if we assume normality) define the distribution converge towards each other. This is what we mean by ergodicity. The distribution of the time series converges to the cross-section distribution as T and S get larger.
If we have a process that is ergodic, this helps us a great deal when we are doing statistics, because it means we can use many of the same tools we use for cross-section statistics (most of what is taught in introductory courses - OLS, etc) for time series as well.
But it is also a warning for using standard statistical tools on data that is non-ergodic. Implicitly when you run a regression of one variable on another and those series are time series, then you are implicitly assuming ergodicity.
Are random walks ergodic? Explain and show computations to support
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.
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 error terms. 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).
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.