diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 9517b4113..9a1fd3161 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -3,10 +3,12 @@ jupytext: text_representation: extension: .md format_name: myst + format_version: 0.13 + jupytext_version: 1.17.2 kernelspec: - display_name: Python 3 - language: python name: python3 + display_name: Python 3 (ipykernel) + language: python --- ```{raw} jupyter @@ -25,10 +27,9 @@ kernelspec: In addition to what's in Anaconda, this lecture will need the following libraries: -```{code-cell} ipython ---- -tags: [hide-output] ---- +```{code-cell} ipython3 +:tags: [hide-output] + !pip install quantecon ``` @@ -46,15 +47,14 @@ Our objectives are to Let's start with some standard imports: -```{code-cell} ipython +```{code-cell} ipython3 import matplotlib.pyplot as plt -plt.rcParams["figure.figsize"] = (11, 5) #set default figure size import numpy as np ``` We'll also use the following for various tasks described below: -```{code-cell} ipython +```{code-cell} ipython3 from quantecon import LinearStateSpace import cmath import math @@ -63,7 +63,7 @@ from sympy import Symbol, init_printing from cmath import sqrt ``` -### Samuelson's Model +### Samuelson's model Samuelson used a *second-order linear difference equation* to represent a model of national output based on three components: @@ -201,7 +201,7 @@ no random shocks hit aggregate demand --- has only transient fluctuations. We can convert the model to one that has persistent irregular fluctuations by adding a random shock to aggregate demand. -### Stochastic Version of the Model +### Stochastic version of the model We create a **random** or **stochastic** version of the model by adding a random process of **shocks** or **disturbances** @@ -212,10 +212,10 @@ equation**: ```{math} :label: second_stochastic -Y_t = G_t + a (1-b) Y_{t-1} - a b Y_{t-2} + \sigma \epsilon_{t} +Y_t = (a+b) Y_{t-1} - b Y_{t-2} + (\gamma + G_t) + \sigma \epsilon_t ``` -### Mathematical Analysis of the Model +### Mathematical analysis of the model To get started, let's set $G_t \equiv 0$, $\sigma = 0$, and $\gamma = 0$. @@ -305,10 +305,10 @@ Notice that $$ \begin{aligned} - Y_t & = & c_1 (r e^{i \omega})^t + c_2 (r e^{-i \omega})^t \\ - & = & c_1 r^t e^{i\omega t} + c_2 r^t e^{-i \omega t} \\ - & = & c_1 r^t [\cos(\omega t) + i \sin(\omega t) ] + c_2 r^t [\cos(\omega t) - i \sin(\omega t) ] \\ - & = & (c_1 + c_2) r^t \cos(\omega t) + i (c_1 - c_2) r^t \sin(\omega t) + Y_t & = c_1 (r e^{i \omega})^t + c_2 (r e^{-i \omega})^t \\ + & = c_1 r^t e^{i\omega t} + c_2 r^t e^{-i \omega t} \\ + & = c_1 r^t [\cos(\omega t) + i \sin(\omega t) ] + c_2 r^t [\cos(\omega t) - i \sin(\omega t) ] \\ + & = (c_1 + c_2) r^t \cos(\omega t) + i (c_1 - c_2) r^t \sin(\omega t) \end{aligned} $$ @@ -354,7 +354,7 @@ absolute values strictly less than one, the absolute value of the larger one governs the rate of convergence to the steady state of the non stochastic version of the model. -### Things This Lecture Does +### Things this lecture does We write a function to generate simulations of a $\{Y_t\}$ sequence as a function of time. @@ -402,10 +402,9 @@ $$ We'll start by drawing an informative graph from page 189 of {cite}`Sargent1987` -```{code-cell} python3 ---- -tags: [output_scroll] ---- +```{code-cell} ipython3 +:tags: [hide-input, output_scroll] + def param_plot(): """This function creates the graph on page 189 of @@ -495,9 +494,9 @@ difference equation parameter pairs in the Samuelson model are such that: Later we'll present the graph with a red mark showing the particular point implied by the setting of $(a,b)$. -### Function to Describe Implications of Characteristic Polynomial +### Function to describe implications of characteristic polynomial -```{code-cell} python3 +```{code-cell} ipython3 def categorize_solution(ρ1, ρ2): """This function takes values of ρ1 and ρ2 and uses them @@ -517,17 +516,17 @@ therefore damped oscillations') therefore get smooth convergence to a steady state') ``` -```{code-cell} python3 -### Test the categorize_solution function +To test the categorize_solution function, +```{code-cell} ipython3 categorize_solution(1.3, -.4) ``` -### Function for Plotting Paths +### Function for plotting paths A useful function for our work below is -```{code-cell} python3 +```{code-cell} ipython3 def plot_y(function=None): """Function plots path of Y_t""" @@ -540,7 +539,7 @@ def plot_y(function=None): plt.show() ``` -### Manual or "by hand" Root Calculations +### Manual or "by hand" root calculations The following function calculates roots of the characteristic polynomial using high school algebra. @@ -550,10 +549,10 @@ using high school algebra. The function also plots a $Y_t$ starting from initial conditions that we set -```{code-cell} python3 +```{code-cell} ipython3 # This is a 'manual' method -def y_nonstochastic(y_0=100, y_1=80, α=.92, β=.5, γ=10, n=80): +def y_nonstochastic(y_0=100, y_1=80, a=.92, b=.5, γ=10, n=80): """Takes values of parameters and computes the roots of characteristic polynomial. It tells whether they are real or complex and whether they @@ -564,28 +563,28 @@ def y_nonstochastic(y_0=100, y_1=80, α=.92, β=.5, γ=10, n=80): roots = [] - ρ1 = α + β - ρ2 = -β + ρ1 = a + b + ρ2 = -b - print(f'ρ_1 is {ρ1}') - print(f'ρ_2 is {ρ2}') + print(f'ρ_1 is {ρ1:.2f}') + print(f'ρ_2 is {ρ2:.2f}') discriminant = ρ1 ** 2 + 4 * ρ2 if discriminant == 0: roots.append(-ρ1 / 2) print('Single real root: ') - print(''.join(str(roots))) + print(''.join(f"{r:.2f}" for r in roots)) elif discriminant > 0: roots.append((-ρ1 + sqrt(discriminant).real) / 2) roots.append((-ρ1 - sqrt(discriminant).real) / 2) print('Two real roots: ') - print(''.join(str(roots))) + print(' '.join(f"{r:.2f}" for r in roots)) else: roots.append((-ρ1 + sqrt(discriminant)) / 2) roots.append((-ρ1 - sqrt(discriminant)) / 2) print('Two complex roots: ') - print(''.join(str(roots))) + print(' '.join(f"{r.real:.2f}{r.imag:+.2f}j" for r in roots)) if all(abs(root) < 1 for root in roots): print('Absolute values of roots are less than one') @@ -604,7 +603,7 @@ def y_nonstochastic(y_0=100, y_1=80, α=.92, β=.5, γ=10, n=80): plot_y(y_nonstochastic()) ``` -### Reverse-Engineering Parameters to Generate Damped Cycles +### Reverse-engineering parameters to generate damped cycles The next cell writes code that takes as inputs the modulus $r$ and phase $\phi$ of a conjugate pair of complex numbers in polar form @@ -618,11 +617,7 @@ $$ - It then reverse-engineers $(a,b)$ and $(\rho_1, \rho_2)$, pairs that would generate those roots -```{code-cell} python3 -### code to reverse-engineer a cycle -### y_t = r^t (c_1 cos(ϕ t) + c2 sin(ϕ t)) -### - +```{code-cell} ipython3 def f(r, ϕ): """ Takes modulus r and angle ϕ of complex number r exp(j ϕ) @@ -639,51 +634,51 @@ def f(r, ϕ): b = -ρ2 # Reverse-engineer a and b that validate these a = ρ1 - b return ρ1, ρ2, a, b +``` -## Now let's use the function in an example -## Here are the example parameters +Now let's use the function in an example. Here are the example parameters: +```{code-cell} ipython3 r = .95 -period = 10 # Length of cycle in units of time +period = 10 # Length of cycle in units of time ϕ = 2 * math.pi/period ## Apply the function - ρ1, ρ2, a, b = f(r, ϕ) -print(f"a, b = {a}, {b}") -print(f"ρ1, ρ2 = {ρ1}, {ρ2}") +print(f"a, b = {a:.2f}, {b:.2f}") +print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ``` -```{code-cell} python3 +```{code-cell} ipython3 ## Print the real components of ρ1 and ρ2 ρ1 = ρ1.real ρ2 = ρ2.real -ρ1, ρ2 +print(f"ρ1 = {ρ1:.2f}, ρ2 = {ρ2:.2f}") ``` -### Root Finding Using Numpy +### Root finding using numpy Here we'll use numpy to compute the roots of the characteristic polynomial -```{code-cell} python3 +```{code-cell} ipython3 r1, r2 = np.roots([1, -ρ1, -ρ2]) p1 = cmath.polar(r1) p2 = cmath.polar(r2) -print(f"r, ϕ = {r}, {ϕ}") -print(f"p1, p2 = {p1}, {p2}") -# print(f"g1, g2 = {g1}, {g2}") +print(f"r, ϕ = {r:.2f}, {ϕ:.2f}") +print(f"p1, p2 = ({p1[0]:.2f}, {p1[1]:.2f}), ({p2[0]:.2f}, {p2[1]:.2f})") +# print(f"g1, g2 = {g1:.2f}, {g2:.2f}") -print(f"a, b = {a}, {b}") -print(f"ρ1, ρ2 = {ρ1}, {ρ2}") +print(f"a, b = {a:.2f}, {b:.2f}") +print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ``` -```{code-cell} python3 +```{code-cell} ipython3 ##=== This method uses numpy to calculate roots ===# @@ -731,38 +726,37 @@ def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): plot_y(y_nonstochastic()) ``` -### Reverse-Engineered Complex Roots: Example +### Reverse-engineered complex roots: Example The next cell studies the implications of reverse-engineered complex roots. We'll generate an **undamped** cycle of period 10 -```{code-cell} python3 +```{code-cell} ipython3 r = 1 # Generates undamped, nonexplosive cycles period = 10 # Length of cycle in units of time ϕ = 2 * math.pi/period ## Apply the reverse-engineering function f - ρ1, ρ2, a, b = f(r, ϕ) # Drop the imaginary part so that it is a valid input into y_nonstochastic a = a.real b = b.real -print(f"a, b = {a}, {b}") +print(f"a, b = {a:.2f}, {b:.2f}") ytemp = y_nonstochastic(α=a, β=b, y_0=20, y_1=30) plot_y(ytemp) ``` -### Digression: Using Sympy to Find Roots +### Digression: Using Sympy to find roots We can also use sympy to compute analytic formulas for the roots -```{code-cell} python3 +```{code-cell} ipython3 init_printing() r1 = Symbol("ρ_1") @@ -772,58 +766,58 @@ z = Symbol("z") sympy.solve(z**2 - r1*z - r2, z) ``` -```{code-cell} python3 -a = Symbol("α") -b = Symbol("β") +```{code-cell} ipython3 +a = Symbol("a") +b = Symbol("b") r1 = a + b r2 = -b sympy.solve(z**2 - r1*z - r2, z) ``` -## Stochastic Shocks +## Stochastic shocks Now we'll construct some code to simulate the stochastic version of the model that emerges when we add a random shock process to aggregate demand -```{code-cell} python3 -def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5): - +```{code-cell} ipython3 +def y_stochastic(y_0=0, y_1=0, a=0.8, b=0.2, γ=10, n=100, σ=5): + """This function takes parameters of a stochastic version of the model and proceeds to analyze the roots of the characteristic polynomial and also generate a simulation. """ # Useful constants - ρ1 = α + β - ρ2 = -β + ρ1 = a + b + ρ2 = -b # Categorize solution categorize_solution(ρ1, ρ2) # Find roots of polynomial roots = np.roots([1, -ρ1, -ρ2]) - print(roots) + print(f"Roots are {[f'{root:.2f}' for root in roots]}") # Check if real or complex if all(isinstance(root, complex) for root in roots): - print('Roots are complex') + print("Roots are complex") else: - print('Roots are real') + print("Roots are real") # Check if roots are less than one if all(abs(root) < 1 for root in roots): - print('Roots are less than one') + print("Roots are less than one") else: - print('Roots are not less than one') + print("Roots are not less than one") # Generate shocks ϵ = np.random.normal(0, 1, n) # Define transition equation - def transition(x, t): return ρ1 * \ - x[t - 1] + ρ2 * x[t - 2] + γ + σ * ϵ[t] + def transition(x, t): + return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + σ * ϵ[t] # Set initial conditions y_t = [y_0, y_1] @@ -839,34 +833,33 @@ plot_y(y_stochastic()) Let's do a simulation in which there are shocks and the characteristic polynomial has complex roots -```{code-cell} python3 +```{code-cell} ipython3 r = .97 -period = 10 # Length of cycle in units of time +period = 10 # Length of cycle in units of time ϕ = 2 * math.pi/period -### Apply the reverse-engineering function f - +# Apply the reverse-engineering function f ρ1, ρ2, a, b = f(r, ϕ) # Drop the imaginary part so that it is a valid input into y_nonstochastic a = a.real b = b.real -print(f"a, b = {a}, {b}") -plot_y(y_stochastic(y_0=40, y_1 = 42, α=a, β=b, σ=2, n=100)) +print(f"a, b = {a:.2f}, {b:.2f}") +plot_y(y_stochastic(y_0=40, y_1=42, a=a, b=b, σ=2, n=100)) ``` -## Government Spending +## Government spending This function computes a response to either a permanent or one-off increase in government expenditures -```{code-cell} python3 +```{code-cell} ipython3 def y_stochastic_g(y_0=20, y_1=20, - α=0.8, - β=0.2, + a=0.8, + b=0.2, γ=10, n=100, σ=2, @@ -879,8 +872,8 @@ def y_stochastic_g(y_0=20, """ # Useful constants - ρ1 = α + β - ρ2 = -β + ρ1 = a + b + ρ2 = -b # Categorize solution categorize_solution(ρ1, ρ2) @@ -948,24 +941,24 @@ def y_stochastic_g(y_0=20, A permanent government spending shock can be simulated as follows -```{code-cell} python3 +```{code-cell} ipython3 plot_y(y_stochastic_g(g=10, g_t=20, duration='permanent')) ``` We can also see the response to a one time jump in government expenditures -```{code-cell} python3 +```{code-cell} ipython3 plot_y(y_stochastic_g(g=500, g_t=50, duration='one-off')) ``` -## Wrapping Everything Into a Class +## Wrapping everything into a class Up to now, we have written functions to do the work. Now we'll roll up our sleeves and write a Python class called `Samuelson` for the Samuelson model -```{code-cell} python3 +```{code-cell} ipython3 class Samuelson(): """This class represents the Samuelson model, otherwise known as the @@ -976,7 +969,7 @@ class Samuelson(): .. math:: - Y_t = + \alpha (1 + \beta) Y_{t-1} - \alpha \beta Y_{t-2} + Y_t = a (1 + \beta) Y_{t-1} - a \beta Y_{t-2} Parameters ---------- @@ -984,9 +977,9 @@ class Samuelson(): Initial condition for Y_0 y_1 : scalar Initial condition for Y_1 - α : scalar + a : scalar Marginal propensity to consume - β : scalar + b : scalar Accelerator coefficient n : int Number of iterations @@ -1007,8 +1000,8 @@ class Samuelson(): def __init__(self, y_0=100, y_1=50, - α=1.3, - β=0.2, + a=1.3, + b=0.2, γ=10, n=100, σ=0, @@ -1016,11 +1009,11 @@ class Samuelson(): g_t=0, duration=None): - self.y_0, self.y_1, self.α, self.β = y_0, y_1, α, β + self.y_0, self.y_1, self.a, self.b = y_0, y_1, a, b self.n, self.g, self.g_t, self.duration = n, g, g_t, duration self.γ, self.σ = γ, σ - self.ρ1 = α + β - self.ρ2 = -β + self.ρ1 = a + b + self.ρ2 = -b self.roots = np.roots([1, -self.ρ1, -self.ρ2]) def root_type(self): @@ -1122,7 +1115,7 @@ class Samuelson(): ax.grid() # Add parameter values to plot - paramstr = f'$\\alpha={self.α:.2f}$ \n $\\beta={self.β:.2f}$ \n \ + paramstr = f'$a={self.a:.2f}$ \n $b={self.b:.2f}$ \n \ $\\gamma={self.γ:.2f}$ \n $\\sigma={self.σ:.2f}$ \n \ $\\rho_1={self.ρ1:.2f}$ \n $\\rho_2={self.ρ2:.2f}$' props = dict(fc='white', pad=10, alpha=0.5) @@ -1158,33 +1151,33 @@ class Samuelson(): return fig ``` -### Illustration of Samuelson Class +### Illustration of Samuelson class Now we'll put our Samuelson class to work on an example -```{code-cell} python3 -sam = Samuelson(α=0.8, β=0.5, σ=2, g=10, g_t=20, duration='permanent') +```{code-cell} ipython3 +sam = Samuelson(a=0.8, b=0.5, σ=2, g=10, g_t=20, duration='permanent') sam.summary() ``` -```{code-cell} python3 +```{code-cell} ipython3 sam.plot() plt.show() ``` -### Using the Graph +### Using the graph We'll use our graph to show where the roots lie and how their location is consistent with the behavior of the path just graphed. The red $+$ sign shows the location of the roots -```{code-cell} python3 +```{code-cell} ipython3 sam.param_plot() plt.show() ``` -## Using the LinearStateSpace Class +## Using the LinearStateSpace class It turns out that we can use the [QuantEcon.py](https://quantecon.org/quantecon-py) [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to do @@ -1193,14 +1186,14 @@ much of the work that we have done from scratch above. Here is how we map the Samuelson model into an instance of a `LinearStateSpace` class -```{code-cell} python3 +```{code-cell} ipython3 """This script maps the Samuelson model in the the ``LinearStateSpace`` class """ -α = 0.8 -β = 0.9 -ρ1 = α + β -ρ2 = -β +a = 0.8 +b = 0.9 +ρ1 = a + b +ρ2 = -b γ = 10 σ = 1 g = 10 @@ -1211,8 +1204,8 @@ A = [[1, 0, 0], [0, 1, 0]] G = [[γ + g, ρ1, ρ2], # this is Y_{t+1} - [γ, α, 0], # this is C_{t+1} - [0, β, -β]] # this is I_{t+1} + [γ, a, 0], # this is C_{t+1} + [0, b, -b]] # this is I_{t+1} μ_0 = [1, 100, 50] C = np.zeros((3,1)) @@ -1235,12 +1228,12 @@ axes[-1].set_xlabel('Iteration') plt.show() ``` -### Other Methods in the `LinearStateSpace` Class +### Other methods in the `LinearStateSpace` class Let's plot **impulse response functions** for the instance of the Samuelson model using a method in the `LinearStateSpace` class -```{code-cell} python3 +```{code-cell} ipython3 imres = sam_t.impulse_response() imres = np.asarray(imres) y1 = imres[:, :, 0] @@ -1251,18 +1244,18 @@ y1.shape Now let's compute the zeros of the characteristic polynomial by simply calculating the eigenvalues of $A$ -```{code-cell} python3 +```{code-cell} ipython3 A = np.asarray(A) w, v = np.linalg.eig(A) -print(w) +print(np.round(w, 2)) ``` -### Inheriting Methods from `LinearStateSpace` +### Inheriting methods from `LinearStateSpace` We could also create a subclass of `LinearStateSpace` (inheriting all its methods and attributes) to add more functions to use -```{code-cell} python3 +```{code-cell} ipython3 class SamuelsonLSS(LinearStateSpace): """ @@ -1272,21 +1265,21 @@ class SamuelsonLSS(LinearStateSpace): def __init__(self, y_0=100, y_1=50, - α=0.8, - β=0.9, + a=0.8, + b=0.9, γ=10, σ=1, g=10): - self.α, self.β = α, β + self.a, self.b = a, b self.y_0, self.y_1, self.g = y_0, y_1, g self.γ, self.σ = γ, σ # Define intial conditions self.μ_0 = [1, y_0, y_1] - self.ρ1 = α + β - self.ρ2 = -β + self.ρ1 = a + b + self.ρ2 = -b # Define transition matrix self.A = [[1, 0, 0], @@ -1295,8 +1288,8 @@ class SamuelsonLSS(LinearStateSpace): # Define output matrix self.G = [[γ + g, self.ρ1, self.ρ2], # this is Y_{t+1} - [γ, α, 0], # this is C_{t+1} - [0, β, -β]] # this is I_{t+1} + [γ, a, 0], # this is C_{t+1} + [0, b, -b]] # this is I_{t+1} self.C = np.zeros((3, 1)) self.C[1] = σ # stochastic @@ -1371,30 +1364,30 @@ class SamuelsonLSS(LinearStateSpace): Let's show how we can use the `SamuelsonLSS` -```{code-cell} python3 +```{code-cell} ipython3 samlss = SamuelsonLSS() ``` -```{code-cell} python3 +```{code-cell} ipython3 samlss.plot_simulation(100, stationary=False) plt.show() ``` -```{code-cell} python3 +```{code-cell} ipython3 samlss.plot_simulation(100, stationary=True) plt.show() ``` -```{code-cell} python3 +```{code-cell} ipython3 samlss.plot_irf(100) plt.show() ``` -```{code-cell} python3 +```{code-cell} ipython3 samlss.multipliers() ``` -## Pure Multiplier Model +## Pure multiplier model Let's shut down the accelerator by setting $b=0$ to get a pure multiplier model @@ -1402,23 +1395,23 @@ multiplier model - the absence of cycles gives an idea about why Samuelson included the accelerator -```{code-cell} python3 -pure_multiplier = SamuelsonLSS(α=0.95, β=0) +```{code-cell} ipython3 +pure_multiplier = SamuelsonLSS(a=0.95, b=0) ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier.plot_simulation() ``` -```{code-cell} python3 -pure_multiplier = SamuelsonLSS(α=0.8, β=0) +```{code-cell} ipython3 +pure_multiplier = SamuelsonLSS(a=0.8, b=0) ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier.plot_simulation() ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier.plot_irf(100) ```