diff --git a/in-work/solow.md b/in-work/solow.md index b92b139d2..60b49eb48 100644 --- a/in-work/solow.md +++ b/in-work/solow.md @@ -3,94 +3,98 @@ jupytext: text_representation: extension: .md format_name: myst + format_version: 0.13 + jupytext_version: 1.14.1 kernelspec: - display_name: Python 3 + display_name: Python 3 (ipykernel) language: python name: python3 --- # The Solow-Swan Growth Model -We consider a model due -to Robert Solow (1925--2014) and Trevor Swan (1918--1989). +In this lecture we review a famous model due +to [Robert Solow (1925--2014)](https://en.wikipedia.org/wiki/Robert_Solow) and [Trevor Swan (1918--1989)](https://en.wikipedia.org/wiki/Trevor_Swan). -The Solow--Swan economy contains a collection of identical agents, each of -whom saves the same fixed fraction of their current incomes. +The model is used to study growth over the long run. -Savings sustain or increase the stock of capital. - -Capital is combined with labor to produce -output, which in turn is paid out to workers and owners of capital. - -To keep things simple, we ignore population and productivity growth. +Although the model is simple, it contains some interesting lessons. We will use the following imports -```{code-cell} python3 +```{code-cell} ipython3 import matplotlib.pyplot as plt plt.rcParams["figure.figsize"] = (11, 5) #set default figure size import numpy as np ``` - ## The Model +In a Solow--Swan economy, agents save a fixed fraction of their current +incomes. + +Savings sustain or increase the stock of capital. + +Capital is combined with labor to produce output, which in turn is paid out to +workers and owners of capital. + +To keep things simple, we ignore population and productivity growth. + For each integer $t \geq 0$, output $Y_t$ in period $t$ is given by $Y_t = F(K_t, L_t)$, where $K_t$ is capital, $L_t$ is labor and $F$ is an aggregate -production function. +production function. The function $F$ is assumed to be nonnegative and **homogeneous of degree one**, meaning that -% + $$ F(\lambda K, \lambda L) = \lambda F(K, L) \quad \text{for all } \lambda \geq 0. $$ -Production functions with this property include -% -* the \navy{Cobb-Douglas}\index{Cobb-Douglas} function $F(K, L) = A K^{\alpha} +Production functions with this property include + +* the **Cobb-Douglas** function $F(K, L) = A K^{\alpha} L^{1-\alpha}$ with $0 \leq \alpha \leq 1$ and -* the \navy{CES} function $F(K, L) = \left\{ a K^\rho + b L^\rho \right\}^{1/\rho}$ +* the **CES** function $F(K, L) = \left\{ a K^\rho + b L^\rho \right\}^{1/\rho}$ with $a, b, \rho > 0$. We assume a closed economy, so domestic investment equals aggregate domestic -saving. +saving. -The saving rate is a positive constant $s$, so that aggregate -investment and saving both equal $s Y_t$. +The saving rate is a constant $s$ satisfying $0 < s < 1$, so that aggregate +investment and saving both equal $s Y_t$. Capital depreciates: without replenishing through investment, one unit of capital today -becomes $1-\delta$ units tomorrow. +becomes $1-\delta$ units tomorrow. + +Thus, -Thus, -% $$ K_{t+1} = s F(K_t, L_t) + (1 - \delta) K_t. $$ -% -Without population growth, $L_t$ equals some constant $L$. + +Without population growth, $L_t$ equals some constant $L$. Setting $k_t := K_t / L$ and using homogeneity of degree one now yields -% -\begin{equation*} - k_{t+1} + +$$ + k_{t+1} = s \frac{F(K_t, L)}{L} + (1 - \delta) k_t = s F(k_t, 1) + (1 - \delta) k_t. -\end{equation*} -% +$$ + With $f(k) := F(k, 1)$, the final expression for capital dynamics is -% -\begin{equation} - \label{eq:solow} - k_{t+1} = g(k_t) + +```{math} +:label: solow + k_{t+1} = g(k_t) \text{ where } g(k) := s f(k) + (1 - \delta) k. -\end{equation} -% +``` Our aim is to learn about the evolution of $k_t$ over time, given an exogenous initial capital stock $k_0$. @@ -99,377 +103,518 @@ given an exogenous initial capital stock $k_0$. ## A Graphical Perspective To understand the dynamics of the sequence $(k_t)_{t \geq 0}$ we use a 45 -degree diagram. +degree diagram. To do so, we first need to specify the functional form for $f$ and assign values to the parameters. -We choose the Cobb--Douglas specification $f(k) = A k^\alpha$. We set $A=2.0$, -$\alpha=0.3$, $s=0.3$ and $\delta=0.4$. +We choose the Cobb--Douglas specification $f(k) = A k^\alpha$ and set $A=2.0$, +$\alpha=0.3$, $s=0.3$ and $\delta=0.4$. + +The function $g$ from {eq}`solow` is then plotted, along with the 45 +degree line. -The function $g$ from \eqref{eq:solow} is then plotted, along with the 45 -degree line. ++++ +Let's define the constants. -```{code-cell} python3 +```{code-cell} ipython3 A, s, alpha, delta = 2, 0.3, 0.3, 0.4 x0 = 0.25 -num_arrows = 8 -ts_length = 12 xmin, xmax = 0, 3 -g = lambda k: A * s * k**alpha + (1 - delta) * k +``` -kstar = ((s * A) / delta)**(1/(1 - alpha)) +Now, we define the function $g$ -xgrid = np.linspace(xmin, xmax, 12000) +```{code-cell} ipython3 +def g(A, s, alpha, delta, k): + return A * s * k**alpha + (1 - delta) * k +``` -fig, ax = plt.subplots() +Let's plot the 45 degree diagram of $g$ -#ax.set_xlim(xmin, xmax) -#ax.set_ylim(xmin, xmax) +```{code-cell} ipython3 +def plot45(kstar=None): + xgrid = np.linspace(xmin, xmax, 12000) -lb = r'$g(k) = sAk^{\alpha} + (1 - \delta)k$' -ax.plot(xgrid, g(xgrid), lw=2, alpha=0.6, label=lb) -ax.plot(xgrid, xgrid, 'k-', lw=1, alpha=0.7, label='45') + fig, ax = plt.subplots(figsize=(6, 6)) -fps = (kstar,) + ax.set_xlim(xmin, xmax) -ax.plot(fps, fps, 'go', ms=10, alpha=0.6) + g_values = g(A, s, alpha, delta, xgrid) -ax.annotate(r'$k^* = (sA / \delta)^{(1/(1-\alpha))}$', - xy=(kstar, kstar), - xycoords='data', - xytext=(-40, -60), - textcoords='offset points', - fontsize=14, - arrowprops=dict(arrowstyle="->")) + ymin, ymax = np.min(g_values), np.max(g_values) + ax.set_ylim(ymin, ymax) -ax.legend(loc='upper left', frameon=False, fontsize=12) + lb = r'$g(k) = sAk^{\alpha} + (1 - \delta)k$' + ax.plot(xgrid, g_values, lw=2, alpha=0.6, label=lb) + ax.plot(xgrid, xgrid, 'k-', lw=1, alpha=0.7, label='45') -ax.set_xticks((0, 1, 2, 3)) -ax.set_yticks((0, 1, 2, 3)) + if kstar: + fps = (kstar,) -ax.set_xlabel('$k_t$', fontsize=12) -ax.set_ylabel('$k_{t+1}$', fontsize=12) + ax.plot(fps, fps, 'go', ms=10, alpha=0.6) -plt.show() + ax.annotate(r'$k^* = (sA / \delta)^{(1/(1-\alpha))}$', + xy=(kstar, kstar), + xycoords='data', + xytext=(-40, -60), + textcoords='offset points', + fontsize=14, + arrowprops=dict(arrowstyle="->")) + + ax.legend(loc='upper left', frameon=False, fontsize=12) + + ax.set_xticks((0, 1, 2, 3)) + ax.set_yticks((0, 1, 2, 3)) + + ax.set_xlabel('$k_t$', fontsize=12) + ax.set_ylabel('$k_{t+1}$', fontsize=12) + + plt.show() +``` + +```{code-cell} ipython3 +plot45() ``` Suppose, at some $k_t$, the value $g(k_t)$ lies strictly above the 45 degree line. -Then we have $k_{t+1} = g(k_t) > k_t$ and capital per worker rises. +Then we have $k_{t+1} = g(k_t) > k_t$ and capital per worker rises. -If $g(k_t) < k_t$ then capital per worker falls. +If $g(k_t) < k_t$ then capital per worker falls. -If $g(k_t) = k_t$, then we are at a \navy{steady state} and $k_t$ remainds constant. +If $g(k_t) = k_t$, then we are at a **steady state** and $k_t$ remainds constant. -A steady state of the model is a fixed point of the mapping $g$. +(A steady state of the model is a [fixed point](https://en.wikipedia.org/wiki/Fixed_point_(mathematics)) of the mapping $g$.) From the shape of the function $g$ in the figure, we see that -there is a unique steady state in $(0, \infty)$. +there is a unique steady state in $(0, \infty)$. It solves $k = s Ak^{\alpha} + (1-\delta)k$ and hence is given by -% -\begin{equation} - \label{eq:kstarss} - k^* := \left( \frac{s A}{\delta} \right)^{1/(1 - \alpha)}. -\end{equation} +```{math} +:label: kstarss + k^* := \left( \frac{s A}{\delta} \right)^{1/(1 - \alpha)}. +``` If initial capital is below $k^*$, then capital increases over time. -If initial capital is above this level, then the reverse is true. +If initial capital is above this level, then the reverse is true. + ++++ -Thus, we can -say that $(k_t)$ converges to $k^*$, regardless of initial capital -$k_0$. +Let's plot the 45 degree diagram to show the $k^*$ in the plot + +```{code-cell} ipython3 +kstar = ((s * A) / delta)**(1/(1 - alpha)) +plot45(kstar) +``` + + +From our graphical analysis, it appears that $(k_t)$ converges to $k^*$, regardless of initial capital +$k_0$. This is a form of global stability. -The next figure shows two time paths for capital, from -two distinct initial conditions, under the parameterization listed above. +The next figure shows three time paths for capital, from +three distinct initial conditions, under the parameterization listed above. At this parameterization, $k^* \approx 1.78$. -As expected, the time paths in the figure both converge to this value. ++++ + +Let's define the constants and three distinct intital conditions -```{code-cell} python3 +```{code-cell} ipython3 A, s, alpha, delta = 2, 0.3, 0.3, 0.4 x0 = np.array([.25, 1.25, 3.25]) -ts_length = 20 -ts = np.zeros(ts_length) +ts_length = 20 xmin, xmax = 0, ts_length ymin, ymax = 0, 3.5 +``` -g = lambda k: A * s * k**alpha + (1 - delta) * k -# - - -k_star = (s * A / delta)**(1/(1-alpha)) -k_star +```{code-cell} ipython3 +def simulate_ts(x0_values, ts_length): -# + -fig, ax = plt.subplots() -ax.set_xlim(xmin, xmax) -ax.set_ylim(ymin, ymax) + k_star = (s * A / delta)**(1/(1-alpha)) + fig, ax = plt.subplots() + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + ts = np.zeros(ts_length) -# simulate and plot time series -for x_init, col in zip(x0, colors): - ts[0] = x_init - for t in range(1, ts_length): - ts[t] = g(ts[t-1]) - ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6, label=r'$k_0=%g$' %x_init) - -ax.legend(fontsize=12) + # simulate and plot time series + for x_init in x0_values: + ts[0] = x_init + for t in range(1, ts_length): + ts[t] = g(A, s, alpha, delta, ts[t-1]) + ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6, + label=r'$k_0=%g$' %x_init) + ax.plot(np.arange(ts_length), np.full(ts_length,k_star), + alpha=0.6, color='red', label=r'$k_*$') + ax.legend(fontsize=10) -ax.set_xlabel(r'$t$', fontsize=14) -ax.set_ylabel(r'$k_t$', fontsize=14) + ax.set_xlabel(r'$t$', fontsize=14) + ax.set_ylabel(r'$k_t$', fontsize=14) -plt.show() + plt.show() ``` +```{code-cell} ipython3 +simulate_ts(x0, ts_length) +``` +As expected, the time paths in the figure both converge to this value. ++++ ## Growth in Continuous Time In this section we investigate a continuous time version of the Solow--Swan -growth model. +growth model. We will see how the smoothing provided by continuous time can simplify analysis. Recall that the discrete time dynamics for capital are -given by $k_{t+1} = s f(k_t) + (1 - \delta) k_t$. +given by $k_{t+1} = s f(k_t) + (1 - \delta) k_t$. A simple rearrangement gives the rate of change per unit of time: -% -\begin{equation} + +$$ \Delta k_t = s f(k_t) - \delta k_t - \quad \text{where} \quad + \quad \text{where} \quad \Delta k_t := k_{t+1} - k_t. -\end{equation} -% +$$ Taking the time step to zero gives the continuous time limit -% -\begin{equation} - \label{eq:solowc} + +```{math} +:label: solowc k'_t = s f(k_t) - \delta k_t \qquad \text{with} \qquad - k'_t := \frac{\diff}{\diff t} k_t. -\end{equation} -% + k'_t := \frac{d}{dt} k_t. +``` Our aim is to learn about the evolution of $k_t$ over time, given initial stock $k_0$. -A **steady state** for~\eqref{eq:solowc} is a value $k^*$ +A **steady state** for {eq}`solowc` is a value $k^*$ at which capital is unchanging, meaning $k'_t = 0$ or, equivalently, -$s f(k^*) = \delta k^*$. - -As with the discrete time case, when -$f(0)=0$ there is a trivial steady state at $k^*=0$. - -Hence we restrict the -state space for capital to $(0, \infty)$. +$s f(k^*) = \delta k^*$. -Let's also assume -$f(k) = Ak^\alpha$, so $k^*$ solves -$s A k^\alpha = \delta k$. +We assume +$f(k) = Ak^\alpha$, so $k^*$ solves +$s A k^\alpha = \delta k$. -The solution is the same as the discrete time case---see~\eqref{eq:kstarss}. +The solution is the same as the discrete time case---see {eq}`kstarss`. -The dynamics of the Cobb--Douglas case are represented in +The dynamics are represented in the next figure, maintaining the parameterization we used -above. +above. Writing $k'_t = g(k_t)$ with $g(k) = s Ak^\alpha - \delta k$, values of $k$ with $g(k) > 0$ imply that $k'_t > 0$, so -capital is increasing. +capital is increasing. When $g(k) < 0$, the opposite occurs. Once again, high marginal returns to savings at low levels of capital combined with low rates of return at high levels of capital combine to yield global stability. ++++ -```{code-cell} python3 +To see this in a figure, let's define the constants + +```{code-cell} ipython3 A, s, alpha, delta = 2, 0.3, 0.3, 0.4 +``` +Next we define the function $g$ for growth in continuous time -k_grid = np.linspace(0, 2.8, 10000) +```{code-cell} ipython3 +def g_con(A, s, alpha, delta, k): + return A * s * k**alpha - delta * k +``` -g = lambda k: A * s * k**alpha + - delta * k +```{code-cell} ipython3 +def plot_gcon(kstar=None): -kstar = ((s * A) / delta)**(1/(1 - alpha)) + k_grid = np.linspace(0, 2.8, 10000) -# + -fig, ax = plt.subplots() + fig, ax = plt.subplots() + ax.plot(k_grid, g_con(A, s, alpha, delta, k_grid), label='$g(k)$') + ax.plot(k_grid, 0 * k_grid, label="$k'=0$") -ax.plot(k_grid, g(k_grid), label='$g(k)$') -ax.plot(k_grid, 0 * k_grid, label="$k'=0$") + if kstar: + fps = (kstar,) -fps = (kstar,) + ax.plot(fps, 0, 'go', ms=10, alpha=0.6) -ax.plot(fps, 0, 'go', ms=10, alpha=0.6) + ax.annotate(r'$k^* = (sA / \delta)^{(1/(1-\alpha))}$', + xy=(kstar, 0), + xycoords='data', + xytext=(0, 60), + textcoords='offset points', + fontsize=12, + arrowprops=dict(arrowstyle="->")) -ax.annotate(r'$k^* = (sA / \delta)^{(1/(1-\alpha))}$', - xy=(kstar, 0), - xycoords='data', - xytext=(0, 60), - textcoords='offset points', - fontsize=12, - arrowprops=dict(arrowstyle="->")) + ax.legend(loc='lower left', fontsize=12) -ax.legend(loc='lower left', fontsize=12) + ax.set_xlabel("$k$",fontsize=10) + ax.set_ylabel("$k'$", fontsize=10) -ax.set_xlabel("$k$",fontsize=12) -ax.set_ylabel("$k'$", fontsize=12) + ax.set_xticks((0, 1, 2, 3)) + ax.set_yticks((-0.3, 0, 0.3)) -ax.set_xticks((0, 1, 2, 3)) -ax.set_yticks((-0.3, 0, 0.3)) + plt.show() +``` -plt.show() +```{code-cell} ipython3 +kstar = ((s * A) / delta)**(1/(1 - alpha)) +plot_gcon(kstar) ``` This shows global stability heuristically for a fixed parameterization, but how would we show the same thing formally for a continuum of plausible parameters? -In the discrete time case, -a neat expression for $k_t$ is hard to obtain. +In the discrete time case, a neat expression for $k_t$ is hard to obtain. In continuous time the process is easier: we can obtain a relatively simple -expression for $k_t$ that specifies the entire path. +expression for $k_t$ that specifies the entire path. The first step is to set $x_t := k_t^{1-\alpha}$, so that $x'_t = (1-\alpha) k_t^{-\alpha} -k'_t$. +k'_t$. Substituting into $k'_t = sAk_t^\alpha - \delta k_t$ leads to the linear differential equation -% -\begin{equation}\label{eq:xsolow} + +```{math} +:label: xsolow x'_t = (1-\alpha) (sA - \delta x_t). -\end{equation} -% +``` This equation has the exact solution -% -\begin{equation*} - x_t + +$$ + x_t = \left( - k_0^{1-\alpha} - \frac{sA}{\delta} + k_0^{1-\alpha} - \frac{sA}{\delta} \right) - \me^{-\delta (1-\alpha) t} + + \mathrm{e}^{-\delta (1-\alpha) t} + \frac{sA}{\delta}. -\end{equation*} -% +$$ -(You can confirm that this function $x_t$ satisfies~\eqref{eq:xsolow} by -differentiating it.) +(You can confirm that this function $x_t$ satisfies {eq}`xsolow` by +differentiating it with respect to $t$.) Converting back to $k_t$ yields -% -\begin{equation}\label{eq:ssivs} - k_t - = + +```{math} +:label: ssivs + k_t + = \left[ \left( - k_0^{1-\alpha} - \frac{sA}{\delta} + k_0^{1-\alpha} - \frac{sA}{\delta} \right) - \me^{-\delta (1-\alpha) t} + + \mathrm{e}^{-\delta (1-\alpha) t} + \frac{sA}{\delta} \right]^{1/(1-\alpha)}. -\end{equation} -% +``` Since $\delta > 0$ and $\alpha \in (0, 1)$, we see immediately that $k_t \to k^*$ as $t \to \infty$ independent of $k_0$. Thus, global stability holds. ++++ + +## Exercises + +```{exercise-start} +:label: ex1 +``` +### Exercise 1 + +plot per capita consumption $c$ at the steady state, as a function of the savings rate $s$, where $0 \leq s \leq 1$. + +Use the Cobb--Douglas specification $f(k) = A k^\alpha$. + +Set $A=2.0, \alpha=0.3,$ and $\delta=0.5$ + +Also, find the approximate value of $s$ that maximizes the $C^*(s)$ and show it in the plot. + +```{exercise-end} +``` + +```{solution-start} ex1 +:class: dropdown +``` + ++++ + +Steady state consumption at savings rate $s$ is given by + + +$$ + c^*(s) = (1-s)f(k^*) = (1-s)A(k^*)^\alpha +$$ + +```{code-cell} ipython3 +A = 2.0 +alpha = 0.3 +delta = 0.5 +``` + +```{code-cell} ipython3 +s_grid = np.linspace(0, 1, 1000) +k_star = ((s_grid * A) / delta)**(1/(1 - alpha)) +c_star = (1 - s_grid) * A * k_star ** alpha +``` + +Let's find the value of $s$ that maximizes $c^*$. + +```{code-cell} ipython3 +c_max_index = np.argmax(c_star) +s_star_max = s_grid[c_max_index] + +x_s_max = np.array([s_star_max, s_star_max]) +y_s_max = np.array([0, c_star[c_max_index]]) + +fig, ax = plt.subplots() + +fps = (c_star[c_max_index],) + +# Highlight the maximum point with a marker +ax.plot((s_star_max, ), (c_star[c_max_index],), 'go', ms=8, alpha=0.6) +ax.annotate(r'$s^*$', + xy=(s_star_max, c_star[c_max_index]), + xycoords='data', + xytext=(20, -50), + textcoords='offset points', + fontsize=12, + arrowprops=dict(arrowstyle="->")) +ax.plot(s_grid, c_star, label=r'$C*(s)$') +ax.plot(x_s_max, y_s_max, alpha=0.5, ls='dotted') +ax.set_xlabel(r'$s$') +ax.set_ylabel(r'$C^*(s)$') +ax.legend() + +plt.show() +``` + +Incidentally, the rate of savings which maximizes steady state level of per capita consumption is called the [Golden Rule savings rate](https://en.wikipedia.org/wiki/Golden_Rule_savings_rate). + +```{solution-end} +``` + ++++ -## Stochastic Productivity +```{exercise-start} +:label: ex2 +``` +### Exercise 2 + +#### Stochastic Productivity To bring the Solow--Swan model closer to data, we need to think about handling -random fluctuations in aggregate quantities. +random fluctuations in aggregate quantities. Among other things, this will eliminate the unrealistic prediction that per-capita output $y_t = A -k^\alpha_t$ converges to a constant $y^* := A (k^*)^\alpha$. +k^\alpha_t$ converges to a constant $y^* := A (k^*)^\alpha$. We shift to discrete time for the following discussion. One approach is to replace constant productivity with some -stochastic sequence $(A_t)_{t \geq 1}$. +stochastic sequence $(A_t)_{t \geq 1}$. Dynamics are now -% -\begin{equation}\label{eq:solowran} + +```{math} +:label: solowran k_{t+1} = s A_{t+1} f(k_t) + (1 - \delta) k_t. -\end{equation} -% +``` -The next figure shows some time series generated by this -model when $f$ is Cobb--Douglas, $(A_t)$ is {\sc iid} and lognormal, and other -parameters are as above. +We suppose $f$ is Cobb--Douglas and $(A_t)$ is IID and lognormal. -Now the long run convergence obtained in in the deterministic case breaks +Now the long run convergence obtained in the deterministic case breaks down, since the system is hit with new shocks at each point in time. +Consider $A=2.0, s=0.6, \alpha=0.3,$ and $\delta=0.5$ + +Generate and plot the time series $k_t$. + +```{exercise-end} +``` -```{code-cell} python3 ++++ + +```{solution-start} ex2 +:class: dropdown +``` + +Let's define the constants for lognormal distribution and initial values used for simulation + +```{code-cell} ipython3 +# Define the constants sig = 0.2 mu = np.log(2) - sig**2 / 2 +A = 2.0 +s = 0.6 +alpha = 0.3 +delta = 0.5 +x0 = [.25, 3.25] # list of initial values used for simulation +``` +Let's define the function *k_next* to find the next value of $k$ + +```{code-cell} ipython3 def lgnorm(): return np.exp(mu + sig * np.random.randn()) -def G(k): - return lgnorm() * s * k**alpha + (1 - delta) * k - -fig, ax = plt.subplots() - +def k_next(s, alpha, delta, k): + return lgnorm() * s * k**alpha + (1 - delta) * k +``` -colors = ('g', 'b') -x0 = np.array([.25, 3.25]) -ts_length = 50 -ts = np.zeros(ts_length) +```{code-cell} ipython3 +def ts_plot(x_values, ts_length): + fig, ax = plt.subplots() + ts = np.zeros(ts_length) -xmin, xmax = 0, ts_length -ymin, ymax = 0, 3.5 + # simulate and plot time series + for x_init in x_values: + ts[0] = x_init + for t in range(1, ts_length): + ts[t] = k_next(s, alpha, delta, ts[t-1]) + ax.plot(np.arange(ts_length), ts, '-o', ms=4, + alpha=0.6, label=r'$k_0=%g$' %x_init) -ax.set_xlim(xmin, xmax) -ax.set_ylim(ymin, ymax) + ax.legend(loc='best', fontsize=10) -# simulate and plot time series -for x_init, col in zip(x0, colors): - ts[0] = x_init - for t in range(1, ts_length): - ts[t] = G(ts[t-1]) - ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6, label=r'$k_0=%g$' %x_init) - -ax.legend(fontsize=12) + ax.set_xlabel(r'$t$', fontsize=12) + ax.set_ylabel(r'$k_t$', fontsize=12) -ax.set_xlabel(r'$t$', fontsize=14) -ax.set_ylabel(r'$k_t$', fontsize=14) + plt.show() +``` -plt.show() +```{code-cell} ipython3 +ts_plot(x0, 50) ``` -At the same time, if we look at the figure carefully, -we observe that the influence of initial conditions seems to die out, and the -two time series eventually fluctuate in similar ranges. -This hints at the fact that, for this model, stability is not lost after all. -There is a higher-level notion of stability at work here, sometimes called -stochastic stability. +```{solution-end} +``` + +```{code-cell} ipython3 -These notions are discussed in more advanced lectures on Markov processes. +```