From e4fdac9aff210bb6f0e56a7b0ecd1e0ef23d0d68 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Tue, 12 Aug 2025 12:03:47 +0530 Subject: [PATCH 01/17] Fixing the headings to match Quantecon style guidelines and fixing the mathematical equation --- lectures/samuelson.md | 136 +++++++++++++++++++++--------------------- 1 file changed, 68 insertions(+), 68 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 9517b4113..3fae9fcf7 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,7 +47,7 @@ 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 @@ -54,7 +55,7 @@ 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 +64,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 +202,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 +213,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$. @@ -354,7 +355,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 +403,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: [output_scroll] + def param_plot(): """This function creates the graph on page 189 of @@ -495,9 +495,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 +517,17 @@ therefore damped oscillations') therefore get smooth convergence to a steady state') ``` -```{code-cell} python3 +```{code-cell} ipython3 ### Test the categorize_solution function 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 +540,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,7 +550,7 @@ 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): @@ -604,7 +604,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,7 +618,7 @@ $$ - It then reverse-engineers $(a,b)$ and $(\rho_1, \rho_2)$, pairs that would generate those roots -```{code-cell} python3 +```{code-cell} ipython3 ### code to reverse-engineer a cycle ### y_t = r^t (c_1 cos(ϕ t) + c2 sin(ϕ t)) ### @@ -655,7 +655,7 @@ print(f"a, b = {a}, {b}") print(f"ρ1, ρ2 = {ρ1}, {ρ2}") ``` -```{code-cell} python3 +```{code-cell} ipython3 ## Print the real components of ρ1 and ρ2 ρ1 = ρ1.real @@ -664,12 +664,12 @@ print(f"ρ1, ρ2 = {ρ1}, {ρ2}") ρ1, ρ2 ``` -### 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) @@ -683,7 +683,7 @@ print(f"a, b = {a}, {b}") print(f"ρ1, ρ2 = {ρ1}, {ρ2}") ``` -```{code-cell} python3 +```{code-cell} ipython3 ##=== This method uses numpy to calculate roots ===# @@ -731,14 +731,14 @@ 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 @@ -758,11 +758,11 @@ 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,7 +772,7 @@ z = Symbol("z") sympy.solve(z**2 - r1*z - r2, z) ``` -```{code-cell} python3 +```{code-cell} ipython3 a = Symbol("α") b = Symbol("β") r1 = a + b @@ -781,13 +781,13 @@ 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 +```{code-cell} ipython3 def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5): """This function takes parameters of a stochastic version of @@ -839,7 +839,7 @@ 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 @@ -857,12 +857,12 @@ print(f"a, b = {a}, {b}") plot_y(y_stochastic(y_0=40, y_1 = 42, α=a, β=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, @@ -948,24 +948,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 +976,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 ---------- @@ -1158,33 +1158,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 +```{code-cell} ipython3 sam = Samuelson(α=0.8, β=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,7 +1193,7 @@ 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 """ @@ -1235,12 +1235,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 +1251,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) ``` -### 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): """ @@ -1371,30 +1371,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 +1402,23 @@ multiplier model - the absence of cycles gives an idea about why Samuelson included the accelerator -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier = SamuelsonLSS(α=0.95, β=0) ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier.plot_simulation() ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier = SamuelsonLSS(α=0.8, β=0) ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier.plot_simulation() ``` -```{code-cell} python3 +```{code-cell} ipython3 pure_multiplier.plot_irf(100) ``` From bfb0b9f911a053732a47b274bc921248cdc8e906 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Tue, 12 Aug 2025 14:40:23 +0530 Subject: [PATCH 02/17] Rectifying the notations in the lecture This commit rectifies the notations and replaces the greek letters (alpha and beta) with (a and b) --- lectures/samuelson.md | 82 +++++++++++++++++++++---------------------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 3fae9fcf7..074fda74d 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -553,7 +553,7 @@ that we set ```{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,8 +564,8 @@ 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}') @@ -687,7 +687,7 @@ print(f"ρ1, ρ2 = {ρ1}, {ρ2}") ##=== This method uses numpy to calculate roots ===# -def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): +def y_nonstochastic(y_0=100, y_1=80, a=.9, b=.8, γ=10, n=80): """ Rather than computing the roots of the characteristic polynomial by hand as we did earlier, this function @@ -695,8 +695,8 @@ def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): """ # Useful constants - ρ1 = α + β - ρ2 = -β + ρ1 = a + b + ρ2 = -b categorize_solution(ρ1, ρ2) @@ -754,7 +754,7 @@ b = b.real print(f"a, b = {a}, {b}") -ytemp = y_nonstochastic(α=a, β=b, y_0=20, y_1=30) +ytemp = y_nonstochastic(a=a, b=b, y_0=20, y_1=30) plot_y(ytemp) ``` @@ -773,8 +773,8 @@ sympy.solve(z**2 - r1*z - r2, z) ``` ```{code-cell} ipython3 -a = Symbol("α") -b = Symbol("β") +a = Symbol("a") +b = Symbol("b") r1 = a + b r2 = -b @@ -788,7 +788,7 @@ model that emerges when we add a random shock process to aggregate demand ```{code-cell} ipython3 -def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5): +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 @@ -796,8 +796,8 @@ def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5): """ # Useful constants - ρ1 = α + β - ρ2 = -β + ρ1 = a + b + ρ2 = -b # Categorize solution categorize_solution(ρ1, ρ2) @@ -854,7 +854,7 @@ 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)) +plot_y(y_stochastic(y_0=40, y_1 = 42, a=a, b=b, σ=2, n=100)) ``` ## Government spending @@ -865,8 +865,8 @@ in government expenditures ```{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 +879,8 @@ def y_stochastic_g(y_0=20, """ # Useful constants - ρ1 = α + β - ρ2 = -β + ρ1 = a + b + ρ2 = -b # Categorize solution categorize_solution(ρ1, ρ2) @@ -984,9 +984,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 +1007,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 +1016,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 +1122,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) @@ -1163,7 +1163,7 @@ class Samuelson(): Now we'll put our Samuelson class to work on an example ```{code-cell} ipython3 -sam = Samuelson(α=0.8, β=0.5, σ=2, g=10, g_t=20, duration='permanent') +sam = Samuelson(a=0.8, b=0.5, σ=2, g=10, g_t=20, duration='permanent') sam.summary() ``` @@ -1197,10 +1197,10 @@ Here is how we map the Samuelson model into an instance of a """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 +1211,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)) @@ -1272,21 +1272,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 +1295,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 @@ -1403,7 +1403,7 @@ multiplier model accelerator ```{code-cell} ipython3 -pure_multiplier = SamuelsonLSS(α=0.95, β=0) +pure_multiplier = SamuelsonLSS(a=0.95, b=0) ``` ```{code-cell} ipython3 @@ -1411,7 +1411,7 @@ pure_multiplier.plot_simulation() ``` ```{code-cell} ipython3 -pure_multiplier = SamuelsonLSS(α=0.8, β=0) +pure_multiplier = SamuelsonLSS(a=0.8, b=0) ``` ```{code-cell} ipython3 From 515848339e35c6d736d1cd311bad36ecf4ab2eb8 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Mon, 18 Aug 2025 20:12:23 +0530 Subject: [PATCH 03/17] Update the math equation of the lecture --- lectures/samuelson.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 074fda74d..fe6a2f176 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -213,7 +213,7 @@ equation**: ```{math} :label: second_stochastic -Y_t = (a+b)Y_{t-1} - b Y_{t-2} + (\gamma + G_t) + \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 From a84faeb3afe290f42ff86362208751c6a3082862 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Wed, 20 Aug 2025 10:59:48 +0530 Subject: [PATCH 04/17] Fixing the math alignments after 'Notice that' --- lectures/samuelson.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index fe6a2f176..202fc9087 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -306,10 +306,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} $$ From 4c78e09b80e581fcc2abfc56b4d07c1fe97d67d5 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Wed, 20 Aug 2025 11:13:03 +0530 Subject: [PATCH 05/17] remove global figure size setting --- lectures/samuelson.md | 1 - 1 file changed, 1 deletion(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 202fc9087..2853c972c 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -49,7 +49,6 @@ Let's start with some standard imports: ```{code-cell} ipython3 import matplotlib.pyplot as plt -plt.rcParams["figure.figsize"] = (11, 5) #set default figure size import numpy as np ``` From 15bf4c7c2b231bf5f8a75fb407242ef8b8a09366 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Thu, 21 Aug 2025 11:10:23 +0530 Subject: [PATCH 06/17] Rectifying the comments on the samuelson lectures --- lectures/samuelson.md | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 2853c972c..71ecbc8be 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -516,9 +516,9 @@ therefore damped oscillations') therefore get smooth convergence to a steady state') ``` -```{code-cell} ipython3 -### Test the categorize_solution function +To test the categorize_solution function, +```{code-cell} ipython3 categorize_solution(1.3, -.4) ``` @@ -618,10 +618,6 @@ $$ pairs that would generate those roots ```{code-cell} ipython3 -### code to reverse-engineer a cycle -### y_t = r^t (c_1 cos(ϕ t) + c2 sin(ϕ t)) -### - def f(r, ϕ): """ Takes modulus r and angle ϕ of complex number r exp(j ϕ) @@ -638,16 +634,16 @@ 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 ϕ = 2 * math.pi/period ## Apply the function - ρ1, ρ2, a, b = f(r, ϕ) print(f"a, b = {a}, {b}") @@ -743,8 +739,7 @@ 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 - +# 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 From 410a575e2d6d2fc09e3dc3010029dba0bb9093c0 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Thu, 21 Aug 2025 14:39:31 +0530 Subject: [PATCH 07/17] Truncating most of the numbers to two values after decimal point --- lectures/samuelson.md | 67 +++++++++++++++++++++---------------------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 71ecbc8be..7883bfd90 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -566,25 +566,25 @@ def y_nonstochastic(y_0=100, y_1=80, a=.92, b=.5, γ=10, n=80): ρ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') @@ -640,14 +640,14 @@ 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} ipython3 @@ -656,7 +656,7 @@ print(f"ρ1, ρ2 = {ρ1}, {ρ2}") ρ1 = ρ1.real ρ2 = ρ2.real -ρ1, ρ2 +print(f"ρ1 = {ρ1:.2f}, ρ2 = {ρ2:.2f}") ``` ### Root finding using numpy @@ -670,19 +670,19 @@ 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} ipython3 ##=== This method uses numpy to calculate roots ===# -def y_nonstochastic(y_0=100, y_1=80, a=.9, b=.8, γ=10, n=80): +def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): """ Rather than computing the roots of the characteristic polynomial by hand as we did earlier, this function @@ -690,8 +690,8 @@ def y_nonstochastic(y_0=100, y_1=80, a=.9, b=.8, γ=10, n=80): """ # Useful constants - ρ1 = a + b - ρ2 = -b + ρ1 = α + β + ρ2 = -β categorize_solution(ρ1, ρ2) @@ -739,16 +739,16 @@ 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 +## 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=a, b=b, y_0=20, y_1=30) +ytemp = y_nonstochastic(α=a, β=b, y_0=20, y_1=30) plot_y(ytemp) ``` @@ -783,7 +783,7 @@ demand ```{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. @@ -798,26 +798,26 @@ def y_stochastic(y_0=0, y_1=0, a=0.8, b=0.2, γ=10, n=100, σ=5): # 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] @@ -836,19 +836,18 @@ Let's do a simulation in which there are shocks and the characteristic polynomia ```{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=a, b=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 @@ -1248,7 +1247,7 @@ calculating the eigenvalues of $A$ ```{code-cell} ipython3 A = np.asarray(A) w, v = np.linalg.eig(A) -print(w) +print(np.round(w, 2)) ``` ### Inheriting methods from `LinearStateSpace` From 4246118e665a4adbe53f127a73bb042b094602a6 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Sat, 23 Aug 2025 16:06:16 +0530 Subject: [PATCH 08/17] Hiding the def param_plots() code input which generates the graph --- lectures/samuelson.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 7883bfd90..9a1fd3161 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -403,7 +403,7 @@ $$ We'll start by drawing an informative graph from page 189 of {cite}`Sargent1987` ```{code-cell} ipython3 -:tags: [output_scroll] +:tags: [hide-input, output_scroll] def param_plot(): From a997e1f2af7858a3f6432325e31829147b2fb7ce Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Fri, 29 Aug 2025 14:35:42 +0530 Subject: [PATCH 09/17] Replace a and b with unicode characters alpha and beta --- lectures/samuelson.md | 148 +++++++++++++++++++++--------------------- 1 file changed, 74 insertions(+), 74 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 9a1fd3161..a4ad092fb 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -118,10 +118,10 @@ Let's assume that * $\{Y_t\}$ is a sequence of levels of national income, yet another endogenous variable. -- $a$ is the marginal propensity to consume in the Keynesian - consumption function $C_t = a Y_{t-1} + \gamma$. -- $b$ is the "accelerator coefficient" in the "investment - accelerator" $I_t = b (Y_{t-1} - Y_{t-2})$. +- $⍺$ is the marginal propensity to consume in the Keynesian + consumption function $C_t = ⍺ Y_{t-1} + \gamma$. +- $β$ is the "accelerator coefficient" in the "investment + accelerator" $I_t = β (Y_{t-1} - Y_{t-2})$. - $\{\epsilon_{t}\}$ is an IID sequence standard normal random variables. - $\sigma \geq 0$ is a "volatility" parameter --- setting $\sigma = 0$ recovers the non-stochastic case @@ -132,7 +132,7 @@ The model combines the consumption function ```{math} :label: consumption -C_t = a Y_{t-1} + \gamma +C_t = ⍺ Y_{t-1} + \gamma ``` with the investment accelerator @@ -140,7 +140,7 @@ with the investment accelerator ```{math} :label: accelerator -I_t = b (Y_{t-1} - Y_{t-2}) +I_t = β (Y_{t-1} - Y_{t-2}) ``` and the national income identity @@ -151,10 +151,10 @@ and the national income identity Y_t = C_t + I_t + G_t ``` -- The parameter $a$ is peoples' *marginal propensity to consume* +- The parameter $⍺$ is peoples' *marginal propensity to consume* out of income - equation {eq}`consumption` asserts that people consume a fraction of - $a \in (0,1)$ of each additional dollar of income. -- The parameter $b > 0$ is the investment accelerator coefficient - equation + $⍺ \in (0,1)$ of each additional dollar of income. +- The parameter $β > 0$ is the investment accelerator coefficient - equation {eq}`accelerator` asserts that people invest in physical capital when income is increasing and disinvest when it is decreasing. @@ -162,7 +162,7 @@ Equations {eq}`consumption`, {eq}`accelerator`, and {eq}`income_identity` imply the following second-order linear difference equation for national income: $$ -Y_t = (a+b) Y_{t-1} - b Y_{t-2} + (\gamma + G_t) +Y_t = (⍺+β) Y_{t-1} - β Y_{t-2} + (\gamma + G_t) $$ or @@ -173,7 +173,7 @@ or Y_t = \rho_1 Y_{t-1} + \rho_2 Y_{t-2} + (\gamma + G_t) ``` -where $\rho_1 = (a+b)$ and $\rho_2 = -b$. +where $\rho_1 = (⍺+β)$ and $\rho_2 = -β$. To complete the model, we require two **initial conditions**. @@ -184,7 +184,7 @@ $$ Y_{-1} = \bar Y_{-1}, \quad Y_{-2} = \bar Y_{-2} $$ -We'll ordinarily set the parameters $(a,b)$ so that starting from +We'll ordinarily set the parameters $(⍺,β)$ so that starting from an arbitrary pair of initial conditions $(\bar Y_{-1}, \bar Y_{-2})$, national income $Y_t$ converges to a constant value as $t$ becomes large. @@ -212,7 +212,7 @@ equation**: ```{math} :label: second_stochastic -Y_t = (a+b) Y_{t-1} - b Y_{t-2} + (\gamma + G_t) + \sigma \epsilon_t +Y_t = (⍺+β) Y_{t-1} - β Y_{t-2} + (\gamma + G_t) + \sigma \epsilon_t ``` ### Mathematical analysis of the model @@ -341,7 +341,7 @@ We say that $\check p$ is the **period** because in that amount of time the cosi (Draw a cosine function to convince yourself of this please) **Remark:** Following {cite}`Samuelson1939`, we want to choose the parameters -$a, b$ of the model so that the absolute values (of the possibly +$⍺, β$ of the model so that the absolute values (of the possibly complex) roots $\lambda_1, \lambda_2$ of the characteristic polynomial are both strictly less than one: @@ -360,7 +360,7 @@ We write a function to generate simulations of a $\{Y_t\}$ sequence as a functio The function requires that we put in initial conditions for $Y_{-1}, Y_{-2}$. -The function checks that $a, b$ are set so that $\lambda_1, \lambda_2$ are less than +The function checks that $⍺, β$ are set so that $\lambda_1, \lambda_2$ are less than unity in absolute value (also called "modulus"). The function also tells us whether the roots are complex, and, if they are complex, returns both their real and complex parts. @@ -477,7 +477,7 @@ plt.show() ``` The graph portrays regions in which the $(\lambda_1, \lambda_2)$ -root pairs implied by the $(\rho_1 = (a+b), \rho_2 = - b)$ +root pairs implied by the $(\rho_1 = (⍺+β), \rho_2 = - β)$ difference equation parameter pairs in the Samuelson model are such that: - $(\lambda_1, \lambda_2)$ are complex with modulus less than @@ -492,7 +492,7 @@ difference equation parameter pairs in the Samuelson model are such that: convergence to the steady state without damped cycles. Later we'll present the graph with a red mark showing the particular -point implied by the setting of $(a,b)$. +point implied by the setting of $(⍺,β)$. ### Function to describe implications of characteristic polynomial @@ -552,7 +552,7 @@ that we set ```{code-cell} ipython3 # This is a 'manual' method -def y_nonstochastic(y_0=100, y_1=80, a=.92, b=.5, γ=10, n=80): +def y_nonstochastic(y_0=100, y_1=80, α=.92, β=.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 @@ -563,8 +563,8 @@ def y_nonstochastic(y_0=100, y_1=80, a=.92, b=.5, γ=10, n=80): roots = [] - ρ1 = a + b - ρ2 = -b + ρ1 = α + β + ρ2 = -β print(f'ρ_1 is {ρ1:.2f}') print(f'ρ_2 is {ρ2:.2f}') @@ -614,7 +614,7 @@ $$ - The code assumes that these two complex numbers are the roots of the characteristic polynomial -- It then reverse-engineers $(a,b)$ and $(\rho_1, \rho_2)$, +- It then reverse-engineers $(⍺,β)$ and $(\rho_1, \rho_2)$, pairs that would generate those roots ```{code-cell} ipython3 @@ -624,16 +624,16 @@ def f(r, ϕ): and creates ρ1 and ρ2 of characteristic polynomial for which r exp(j ϕ) and r exp(- j ϕ) are complex roots. - Returns the multiplier coefficient a and the accelerator coefficient b + Returns the multiplier coefficient ⍺ and the accelerator coefficient β that verifies those roots. """ g1 = cmath.rect(r, ϕ) # Generate two complex roots g2 = cmath.rect(r, -ϕ) ρ1 = g1 + g2 # Implied ρ1, ρ2 ρ2 = -g1 * g2 - b = -ρ2 # Reverse-engineer a and b that validate these - a = ρ1 - b - return ρ1, ρ2, a, b + β = -ρ2 # Reverse-engineer a and b that validate these + α = ρ1 - β + return ρ1, ρ2, α, β ``` Now let's use the function in an example. Here are the example parameters: @@ -644,9 +644,9 @@ period = 10 # Length of cycle in units of time ϕ = 2 * math.pi/period ## Apply the function -ρ1, ρ2, a, b = f(r, ϕ) +ρ1, ρ2, α, β = f(r, ϕ) -print(f"a, b = {a:.2f}, {b:.2f}") +print(f"α, β = {α:.2f}, {β:.2f}") print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ``` @@ -674,7 +674,7 @@ 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:.2f}, {b:.2f}") +print(f"α, β = {α:.2f}, {β:.2f}") print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ``` @@ -740,15 +740,15 @@ 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, ϕ) +ρ1, ρ2, α, β = f(r, ϕ) # Drop the imaginary part so that it is a valid input into y_nonstochastic -a = a.real -b = b.real +α = α.real +β = β.real -print(f"a, b = {a:.2f}, {b:.2f}") +print(f"a, b = {α:.2f}, {β:.2f}") -ytemp = y_nonstochastic(α=a, β=b, y_0=20, y_1=30) +ytemp = y_nonstochastic(α=α, β=β, y_0=20, y_1=30) plot_y(ytemp) ``` @@ -767,10 +767,10 @@ sympy.solve(z**2 - r1*z - r2, z) ``` ```{code-cell} ipython3 -a = Symbol("a") -b = Symbol("b") -r1 = a + b -r2 = -b +α = Symbol("α") +β = Symbol("β") +r1 = α + β +r2 = -β sympy.solve(z**2 - r1*z - r2, z) ``` @@ -782,7 +782,7 @@ model that emerges when we add a random shock process to aggregate demand ```{code-cell} ipython3 -def y_stochastic(y_0=0, y_1=0, a=0.8, b=0.2, γ=10, n=100, σ=5): +def y_stochastic(y_0=0, y_1=0, α=0.8, β=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 @@ -790,8 +790,8 @@ def y_stochastic(y_0=0, y_1=0, a=0.8, b=0.2, γ=10, n=100, σ=5): """ # Useful constants - ρ1 = a + b - ρ2 = -b + ρ1 = α + β + ρ2 = -β # Categorize solution categorize_solution(ρ1, ρ2) @@ -840,14 +840,14 @@ 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, ϕ) +ρ1, ρ2, α, β = f(r, ϕ) # Drop the imaginary part so that it is a valid input into y_nonstochastic -a = a.real -b = b.real +α = α.real +β = β.real -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)) +print(f"a, b = {α:.2f}, {β:.2f}") +plot_y(y_stochastic(y_0=40, y_1=42, α=α, β=β, σ=2, n=100)) ``` ## Government spending @@ -858,8 +858,8 @@ in government expenditures ```{code-cell} ipython3 def y_stochastic_g(y_0=20, y_1=20, - a=0.8, - b=0.2, + α=0.8, + β=0.2, γ=10, n=100, σ=2, @@ -872,8 +872,8 @@ def y_stochastic_g(y_0=20, """ # Useful constants - ρ1 = a + b - ρ2 = -b + ρ1 = α + β + ρ2 = -β # Categorize solution categorize_solution(ρ1, ρ2) @@ -969,7 +969,7 @@ class Samuelson(): .. math:: - Y_t = a (1 + \beta) Y_{t-1} - a \beta Y_{t-2} + Y_t = α (1 + β) Y_{t-1} - α β Y_{t-2} Parameters ---------- @@ -977,9 +977,9 @@ class Samuelson(): Initial condition for Y_0 y_1 : scalar Initial condition for Y_1 - a : scalar + α : scalar Marginal propensity to consume - b : scalar + β : scalar Accelerator coefficient n : int Number of iterations @@ -1000,8 +1000,8 @@ class Samuelson(): def __init__(self, y_0=100, y_1=50, - a=1.3, - b=0.2, + α=1.3, + β=0.2, γ=10, n=100, σ=0, @@ -1009,11 +1009,11 @@ class Samuelson(): g_t=0, duration=None): - self.y_0, self.y_1, self.a, self.b = y_0, y_1, a, b + self.y_0, self.y_1, self.α, self.β = y_0, y_1, α, β self.n, self.g, self.g_t, self.duration = n, g, g_t, duration self.γ, self.σ = γ, σ - self.ρ1 = a + b - self.ρ2 = -b + self.ρ1 = α + β + self.ρ2 = -β self.roots = np.roots([1, -self.ρ1, -self.ρ2]) def root_type(self): @@ -1115,7 +1115,7 @@ class Samuelson(): ax.grid() # Add parameter values to plot - paramstr = f'$a={self.a:.2f}$ \n $b={self.b:.2f}$ \n \ + paramstr = f'$α={self.α:.2f}$ \n $β={self.β:.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) @@ -1156,7 +1156,7 @@ class Samuelson(): Now we'll put our Samuelson class to work on an example ```{code-cell} ipython3 -sam = Samuelson(a=0.8, b=0.5, σ=2, g=10, g_t=20, duration='permanent') +sam = Samuelson(α=0.8, β=0.5, σ=2, g=10, g_t=20, duration='permanent') sam.summary() ``` @@ -1190,10 +1190,10 @@ Here is how we map the Samuelson model into an instance of a """This script maps the Samuelson model in the the ``LinearStateSpace`` class """ -a = 0.8 -b = 0.9 -ρ1 = a + b -ρ2 = -b +α = 0.8 +β = 0.9 +ρ1 = α + β +ρ2 = -β γ = 10 σ = 1 g = 10 @@ -1204,8 +1204,8 @@ A = [[1, 0, 0], [0, 1, 0]] G = [[γ + g, ρ1, ρ2], # this is Y_{t+1} - [γ, a, 0], # this is C_{t+1} - [0, b, -b]] # this is I_{t+1} + [γ, α, 0], # this is C_{t+1} + [0, β, -β]] # this is I_{t+1} μ_0 = [1, 100, 50] C = np.zeros((3,1)) @@ -1265,21 +1265,21 @@ class SamuelsonLSS(LinearStateSpace): def __init__(self, y_0=100, y_1=50, - a=0.8, - b=0.9, + α=0.8, + β=0.9, γ=10, σ=1, g=10): - self.a, self.b = a, b + self.α, self.β = α, β 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 = a + b - self.ρ2 = -b + self.ρ1 = α + β + self.ρ2 = -β # Define transition matrix self.A = [[1, 0, 0], @@ -1288,8 +1288,8 @@ class SamuelsonLSS(LinearStateSpace): # Define output matrix self.G = [[γ + g, self.ρ1, self.ρ2], # this is Y_{t+1} - [γ, a, 0], # this is C_{t+1} - [0, b, -b]] # this is I_{t+1} + [γ, α, 0], # this is C_{t+1} + [0, β, -β]] # this is I_{t+1} self.C = np.zeros((3, 1)) self.C[1] = σ # stochastic @@ -1396,7 +1396,7 @@ multiplier model accelerator ```{code-cell} ipython3 -pure_multiplier = SamuelsonLSS(a=0.95, b=0) +pure_multiplier = SamuelsonLSS(α=0.95, β=0) ``` ```{code-cell} ipython3 @@ -1404,7 +1404,7 @@ pure_multiplier.plot_simulation() ``` ```{code-cell} ipython3 -pure_multiplier = SamuelsonLSS(a=0.8, b=0) +pure_multiplier = SamuelsonLSS(α=0.8, β=0) ``` ```{code-cell} ipython3 From 290d43fe48136cdc7738dcdf6ea167503bd053e3 Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Mon, 8 Sep 2025 08:28:43 +0530 Subject: [PATCH 10/17] replace unicode characters in text with latex characters --- lectures/samuelson.md | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index a4ad092fb..2be1d2c71 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -118,9 +118,9 @@ Let's assume that * $\{Y_t\}$ is a sequence of levels of national income, yet another endogenous variable. -- $⍺$ is the marginal propensity to consume in the Keynesian +- $\alpha$ is the marginal propensity to consume in the Keynesian consumption function $C_t = ⍺ Y_{t-1} + \gamma$. -- $β$ is the "accelerator coefficient" in the "investment +- $\beta$ is the "accelerator coefficient" in the "investment accelerator" $I_t = β (Y_{t-1} - Y_{t-2})$. - $\{\epsilon_{t}\}$ is an IID sequence standard normal random variables. - $\sigma \geq 0$ is a "volatility" @@ -151,10 +151,10 @@ and the national income identity Y_t = C_t + I_t + G_t ``` -- The parameter $⍺$ is peoples' *marginal propensity to consume* +- The parameter $\alpha$ is peoples' *marginal propensity to consume* out of income - equation {eq}`consumption` asserts that people consume a fraction of - $⍺ \in (0,1)$ of each additional dollar of income. -- The parameter $β > 0$ is the investment accelerator coefficient - equation + $\alpha \in (0,1)$ of each additional dollar of income. +- The parameter $\beta > 0$ is the investment accelerator coefficient - equation {eq}`accelerator` asserts that people invest in physical capital when income is increasing and disinvest when it is decreasing. @@ -173,7 +173,7 @@ or Y_t = \rho_1 Y_{t-1} + \rho_2 Y_{t-2} + (\gamma + G_t) ``` -where $\rho_1 = (⍺+β)$ and $\rho_2 = -β$. +where $\rho_1 = (\alpha+\beta)$ and $\rho_2 = -\beta$. To complete the model, we require two **initial conditions**. @@ -184,7 +184,7 @@ $$ Y_{-1} = \bar Y_{-1}, \quad Y_{-2} = \bar Y_{-2} $$ -We'll ordinarily set the parameters $(⍺,β)$ so that starting from +We'll ordinarily set the parameters $(\alpha,\beta)$ so that starting from an arbitrary pair of initial conditions $(\bar Y_{-1}, \bar Y_{-2})$, national income $Y_t$ converges to a constant value as $t$ becomes large. @@ -341,7 +341,7 @@ We say that $\check p$ is the **period** because in that amount of time the cosi (Draw a cosine function to convince yourself of this please) **Remark:** Following {cite}`Samuelson1939`, we want to choose the parameters -$⍺, β$ of the model so that the absolute values (of the possibly +$\alpha, \beta$ of the model so that the absolute values (of the possibly complex) roots $\lambda_1, \lambda_2$ of the characteristic polynomial are both strictly less than one: @@ -360,7 +360,7 @@ We write a function to generate simulations of a $\{Y_t\}$ sequence as a functio The function requires that we put in initial conditions for $Y_{-1}, Y_{-2}$. -The function checks that $⍺, β$ are set so that $\lambda_1, \lambda_2$ are less than +The function checks that $\alpha, \beta$ are set so that $\lambda_1, \lambda_2$ are less than unity in absolute value (also called "modulus"). The function also tells us whether the roots are complex, and, if they are complex, returns both their real and complex parts. @@ -477,7 +477,7 @@ plt.show() ``` The graph portrays regions in which the $(\lambda_1, \lambda_2)$ -root pairs implied by the $(\rho_1 = (⍺+β), \rho_2 = - β)$ +root pairs implied by the $(\rho_1 = (\alpha+\beta), \rho_2 = - \beta)$ difference equation parameter pairs in the Samuelson model are such that: - $(\lambda_1, \lambda_2)$ are complex with modulus less than @@ -492,7 +492,7 @@ difference equation parameter pairs in the Samuelson model are such that: convergence to the steady state without damped cycles. Later we'll present the graph with a red mark showing the particular -point implied by the setting of $(⍺,β)$. +point implied by the setting of $(\alpha,\beta)$. ### Function to describe implications of characteristic polynomial @@ -614,7 +614,7 @@ $$ - The code assumes that these two complex numbers are the roots of the characteristic polynomial -- It then reverse-engineers $(⍺,β)$ and $(\rho_1, \rho_2)$, +- It then reverse-engineers $(\alpha,\beta)$ and $(\rho_1, \rho_2)$, pairs that would generate those roots ```{code-cell} ipython3 @@ -624,7 +624,7 @@ def f(r, ϕ): and creates ρ1 and ρ2 of characteristic polynomial for which r exp(j ϕ) and r exp(- j ϕ) are complex roots. - Returns the multiplier coefficient ⍺ and the accelerator coefficient β + Returns the multiplier coefficient $\alpha$ and the accelerator coefficient $\beta$ that verifies those roots. """ g1 = cmath.rect(r, ϕ) # Generate two complex roots @@ -969,7 +969,7 @@ class Samuelson(): .. math:: - Y_t = α (1 + β) Y_{t-1} - α β Y_{t-2} + Y_t = $\alpha$ (1 + $\beta$) Y_{t-1} - $\alpha$ $\beta$ Y_{t-2} Parameters ---------- From 52b6d0165c24c2060b8d0e3617c282b16dc41d5a Mon Sep 17 00:00:00 2001 From: bishmaybarik Date: Mon, 8 Sep 2025 08:39:19 +0530 Subject: [PATCH 11/17] update math blocks unicode characters with latex notations for consistency --- lectures/samuelson.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 2be1d2c71..410096567 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -119,9 +119,9 @@ Let's assume that another endogenous variable. - $\alpha$ is the marginal propensity to consume in the Keynesian - consumption function $C_t = ⍺ Y_{t-1} + \gamma$. + consumption function $C_t = \alpha Y_{t-1} + \gamma$. - $\beta$ is the "accelerator coefficient" in the "investment - accelerator" $I_t = β (Y_{t-1} - Y_{t-2})$. + accelerator" $I_t = \beta (Y_{t-1} - Y_{t-2})$. - $\{\epsilon_{t}\}$ is an IID sequence standard normal random variables. - $\sigma \geq 0$ is a "volatility" parameter --- setting $\sigma = 0$ recovers the non-stochastic case @@ -132,7 +132,7 @@ The model combines the consumption function ```{math} :label: consumption -C_t = ⍺ Y_{t-1} + \gamma +C_t = \alpha Y_{t-1} + \gamma ``` with the investment accelerator @@ -140,7 +140,7 @@ with the investment accelerator ```{math} :label: accelerator -I_t = β (Y_{t-1} - Y_{t-2}) +I_t = \beta (Y_{t-1} - Y_{t-2}) ``` and the national income identity @@ -162,7 +162,7 @@ Equations {eq}`consumption`, {eq}`accelerator`, and {eq}`income_identity` imply the following second-order linear difference equation for national income: $$ -Y_t = (⍺+β) Y_{t-1} - β Y_{t-2} + (\gamma + G_t) +Y_t = (\alpha+\beta) Y_{t-1} - \beta Y_{t-2} + (\gamma + G_t) $$ or @@ -212,7 +212,7 @@ equation**: ```{math} :label: second_stochastic -Y_t = (⍺+β) Y_{t-1} - β Y_{t-2} + (\gamma + G_t) + \sigma \epsilon_t +Y_t = (\alpha+\beta) Y_{t-1} - \beta Y_{t-2} + (\gamma + G_t) + \sigma \epsilon_t ``` ### Mathematical analysis of the model From b5846617dec6198160622c82349e9b3ff89c4860 Mon Sep 17 00:00:00 2001 From: mmcky Date: Thu, 11 Sep 2025 11:32:10 +1000 Subject: [PATCH 12/17] FIX: matplotlib deprecation warning --- lectures/samuelson.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 410096567..7fb21c592 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -4,11 +4,11 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.2 + jupytext_version: 1.16.7 kernelspec: - name: python3 display_name: Python 3 (ipykernel) language: python + name: python3 --- ```{raw} jupyter @@ -1140,10 +1140,10 @@ class Samuelson(): label = rf'$\lambda_{i+1} = {sam.roots[i].real:.2f} {operator[i]} {sam.roots[i].imag:.2f}i$' else: label = rf'$\lambda_{i+1} = {sam.roots[i].real:.2f}$' - ax.scatter(0, 0, 0, label=label) # dummy to add to legend + ax.scatter(0, 0, s=0, label=label) # dummy to add to legend # Add ρ pair to plot - ax.scatter(self.ρ1, self.ρ2, 100, 'red', '+', + ax.scatter(self.ρ1, self.ρ2, s=100, c='red', marker='+', label=r'$(\ \rho_1, \ \rho_2 \ )$', zorder=5) plt.legend(fontsize=12, loc=3) From fac1303350a897af342108333c1e08b8a4b97485 Mon Sep 17 00:00:00 2001 From: mmcky Date: Fri, 12 Sep 2025 09:07:04 +1000 Subject: [PATCH 13/17] update emphasis from style-guide --- lectures/samuelson.md | 37 ++++++++++++++++++------------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 7fb21c592..17bad3493 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -175,7 +175,7 @@ Y_t = \rho_1 Y_{t-1} + \rho_2 Y_{t-2} + (\gamma + G_t) where $\rho_1 = (\alpha+\beta)$ and $\rho_2 = -\beta$. -To complete the model, we require two **initial conditions**. +To complete the model, we require two *initial conditions*. If the model is to generate time series for $t=0, \ldots, T$, we require initial values @@ -192,8 +192,8 @@ a constant value as $t$ becomes large. We are interested in studying - the transient fluctuations in $Y_t$ as it converges to its - **steady state** level -- the **rate** at which it converges to a steady state level + *steady state* level +- the *rate* at which it converges to a steady state level The deterministic version of the model described so far --- meaning that no random shocks hit aggregate demand --- has only transient fluctuations. @@ -203,11 +203,10 @@ fluctuations by adding a random shock to aggregate demand. ### Stochastic version of the model -We create a **random** or **stochastic** version of the model by adding -a random process of **shocks** or **disturbances** +We create a *random* or *stochastic* version of the model by adding +a random process of *shocks* or *disturbances* $\{\sigma \epsilon_t \}$ to the right side of equation {eq}`second_order`, -leading to the **second-order scalar linear stochastic difference -equation**: +leading to the *second-order scalar linear stochastic difference equation*: ```{math} :label: second_stochastic @@ -235,7 +234,7 @@ Y_{t+2} - \rho_1 Y_{t+1} - \rho_2 Y_t = 0 ``` To discover the properties of the solution of {eq}`second_stochastic2`, -it is useful first to form the **characteristic polynomial** +it is useful first to form the *characteristic polynomial* for {eq}`second_stochastic2`: ```{math} @@ -246,7 +245,7 @@ z^2 - \rho_1 z - \rho_2 where $z$ is possibly a complex number. -We want to find the two **zeros** (a.k.a. **roots**) -- namely +We want to find the two *zeros* (a.k.a. *roots*) -- namely $\lambda_1, \lambda_2$ -- of the characteristic polynomial. These are two special values of $z$, say $z= \lambda_1$ and @@ -260,7 +259,7 @@ the characteristic polynomial {eq}`polynomial` equals zero: z^2 - \rho_1 z - \rho_2 = (z- \lambda_1 ) (z -\lambda_2) = 0 ``` -Equation {eq}`polynomial_sol` is said to **factor** the characteristic polynomial. +Equation {eq}`polynomial_sol` is said to *factor* the characteristic polynomial. When the roots are complex, they will occur as a complex conjugate pair. @@ -287,8 +286,8 @@ $$ (To read about the polar form, see [here](https://www.khanacademy.org/math/precalculus/x9e81a4f98389efdf:complex/x9e81a4f98389efdf:complex-mul-div-polar/a/complex-number-polar-form-review)) -Given **initial conditions** $Y_{-1}, Y_{-2}$, we want to generate -a **solution** of the difference equation {eq}`second_stochastic2`. +Given *initial conditions* $Y_{-1}, Y_{-2}$, we want to generate +a *solution* of the difference equation {eq}`second_stochastic2`. It can be represented as @@ -333,14 +332,14 @@ $$ where $v$ and $\theta$ are constants that must be chosen to satisfy initial conditions for $Y_{-1}, Y_{-2}$. This formula shows that when the roots are complex, $Y_t$ displays -oscillations with **period** $\check p = -\frac{2 \pi}{\omega}$ and **damping factor** $r$. +oscillations with *period* $\check p = +\frac{2 \pi}{\omega}$ and *damping factor* $r$. -We say that $\check p$ is the **period** because in that amount of time the cosine wave $\cos(\omega t + \theta)$ goes through exactly one complete cycles. +We say that $\check p$ is the *period* because in that amount of time the cosine wave $\cos(\omega t + \theta)$ goes through exactly one complete cycles. (Draw a cosine function to convince yourself of this please) -**Remark:** Following {cite}`Samuelson1939`, we want to choose the parameters +*Remark:* Following {cite}`Samuelson1939`, we want to choose the parameters $\alpha, \beta$ of the model so that the absolute values (of the possibly complex) roots $\lambda_1, \lambda_2$ of the characteristic polynomial are both strictly less than one: @@ -349,7 +348,7 @@ $$ | \lambda_j | < 1 \quad \quad \text{for } j = 1, 2 $$ -**Remark:** When both roots $\lambda_1, \lambda_2$ of the characteristic polynomial have +*Remark:* When both roots $\lambda_1, \lambda_2$ of the characteristic polynomial have 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. @@ -731,7 +730,7 @@ plot_y(y_nonstochastic()) The next cell studies the implications of reverse-engineered complex roots. -We'll generate an **undamped** cycle of period 10 +We'll generate an *undamped* cycle of period 10 ```{code-cell} ipython3 r = 1 # Generates undamped, nonexplosive cycles @@ -1230,7 +1229,7 @@ plt.show() ### Other methods in the `LinearStateSpace` class -Let's plot **impulse response functions** for the instance of the +Let's plot *impulse response functions* for the instance of the Samuelson model using a method in the `LinearStateSpace` class ```{code-cell} ipython3 From 509a5defc2d12585b9b7bf6d9179210f0fc62e7d Mon Sep 17 00:00:00 2001 From: mmcky Date: Fri, 12 Sep 2025 09:50:48 +1000 Subject: [PATCH 14/17] run black (line-length=80) over notebook --- lectures/samuelson.md | 419 +++++++++++++++++++++++------------------- 1 file changed, 227 insertions(+), 192 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 17bad3493..842109c24 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -27,6 +27,7 @@ kernelspec: In addition to what's in Anaconda, this lecture will need the following libraries: + ```{code-cell} ipython3 :tags: [hide-output] @@ -405,13 +406,12 @@ We'll start by drawing an informative graph from page 189 of {cite}`Sargent1987` :tags: [hide-input, output_scroll] def param_plot(): - """This function creates the graph on page 189 of Sargent Macroeconomic Theory, second edition, 1987. """ fig, ax = plt.subplots(figsize=(10, 6)) - ax.set_aspect('equal') + ax.set_aspect("equal") # Set axis xmin, ymin = -3, -2 @@ -420,57 +420,93 @@ def param_plot(): # Set axis labels ax.set(xticks=[], yticks=[]) - ax.set_xlabel(r'$\rho_2$', fontsize=16) - ax.xaxis.set_label_position('top') - ax.set_ylabel(r'$\rho_1$', rotation=0, fontsize=16) - ax.yaxis.set_label_position('right') + ax.set_xlabel(r"$\rho_2$", fontsize=16) + ax.xaxis.set_label_position("top") + ax.set_ylabel(r"$\rho_1$", rotation=0, fontsize=16) + ax.yaxis.set_label_position("right") # Draw (t1, t2) points ρ1 = np.linspace(-2, 2, 100) - ax.plot(ρ1, -abs(ρ1) + 1, c='black') - ax.plot(ρ1, np.full_like(ρ1, -1), c='black') - ax.plot(ρ1, -(ρ1**2 / 4), c='black') + ax.plot(ρ1, -abs(ρ1) + 1, c="black") + ax.plot(ρ1, np.full_like(ρ1, -1), c="black") + ax.plot(ρ1, -(ρ1**2 / 4), c="black") # Turn normal axes off - for spine in ['left', 'bottom', 'top', 'right']: + for spine in ["left", "bottom", "top", "right"]: ax.spines[spine].set_visible(False) # Add arrows to represent axes - axes_arrows = {'arrowstyle': '<|-|>', 'lw': 1.3} - ax.annotate('', xy=(xmin, 0), xytext=(xmax, 0), arrowprops=axes_arrows) - ax.annotate('', xy=(0, ymin), xytext=(0, ymax), arrowprops=axes_arrows) + axes_arrows = {"arrowstyle": "<|-|>", "lw": 1.3} + ax.annotate("", xy=(xmin, 0), xytext=(xmax, 0), arrowprops=axes_arrows) + ax.annotate("", xy=(0, ymin), xytext=(0, ymax), arrowprops=axes_arrows) # Annotate the plot with equations - plot_arrowsl = {'arrowstyle': '-|>', 'connectionstyle': "arc3, rad=-0.2"} - plot_arrowsr = {'arrowstyle': '-|>', 'connectionstyle': "arc3, rad=0.2"} - ax.annotate(r'$\rho_1 + \rho_2 < 1$', xy=(0.5, 0.3), xytext=(0.8, 0.6), - arrowprops=plot_arrowsr, fontsize='12') - ax.annotate(r'$\rho_1 + \rho_2 = 1$', xy=(0.38, 0.6), xytext=(0.6, 0.8), - arrowprops=plot_arrowsr, fontsize='12') - ax.annotate(r'$\rho_2 < 1 + \rho_1$', xy=(-0.5, 0.3), xytext=(-1.3, 0.6), - arrowprops=plot_arrowsl, fontsize='12') - ax.annotate(r'$\rho_2 = 1 + \rho_1$', xy=(-0.38, 0.6), xytext=(-1, 0.8), - arrowprops=plot_arrowsl, fontsize='12') - ax.annotate(r'$\rho_2 = -1$', xy=(1.5, -1), xytext=(1.8, -1.3), - arrowprops=plot_arrowsl, fontsize='12') - ax.annotate(r'${\rho_1}^2 + 4\rho_2 = 0$', xy=(1.15, -0.35), - xytext=(1.5, -0.3), arrowprops=plot_arrowsr, fontsize='12') - ax.annotate(r'${\rho_1}^2 + 4\rho_2 < 0$', xy=(1.4, -0.7), - xytext=(1.8, -0.6), arrowprops=plot_arrowsr, fontsize='12') + plot_arrowsl = {"arrowstyle": "-|>", "connectionstyle": "arc3, rad=-0.2"} + plot_arrowsr = {"arrowstyle": "-|>", "connectionstyle": "arc3, rad=0.2"} + ax.annotate( + r"$\rho_1 + \rho_2 < 1$", + xy=(0.5, 0.3), + xytext=(0.8, 0.6), + arrowprops=plot_arrowsr, + fontsize="12", + ) + ax.annotate( + r"$\rho_1 + \rho_2 = 1$", + xy=(0.38, 0.6), + xytext=(0.6, 0.8), + arrowprops=plot_arrowsr, + fontsize="12", + ) + ax.annotate( + r"$\rho_2 < 1 + \rho_1$", + xy=(-0.5, 0.3), + xytext=(-1.3, 0.6), + arrowprops=plot_arrowsl, + fontsize="12", + ) + ax.annotate( + r"$\rho_2 = 1 + \rho_1$", + xy=(-0.38, 0.6), + xytext=(-1, 0.8), + arrowprops=plot_arrowsl, + fontsize="12", + ) + ax.annotate( + r"$\rho_2 = -1$", + xy=(1.5, -1), + xytext=(1.8, -1.3), + arrowprops=plot_arrowsl, + fontsize="12", + ) + ax.annotate( + r"${\rho_1}^2 + 4\rho_2 = 0$", + xy=(1.15, -0.35), + xytext=(1.5, -0.3), + arrowprops=plot_arrowsr, + fontsize="12", + ) + ax.annotate( + r"${\rho_1}^2 + 4\rho_2 < 0$", + xy=(1.4, -0.7), + xytext=(1.8, -0.6), + arrowprops=plot_arrowsr, + fontsize="12", + ) # Label categories of solutions - ax.text(1.5, 1, 'Explosive\n growth', ha='center', fontsize=16) - ax.text(-1.5, 1, 'Explosive\n oscillations', ha='center', fontsize=16) - ax.text(0.05, -1.5, 'Explosive oscillations', ha='center', fontsize=16) - ax.text(0.09, -0.5, 'Damped oscillations', ha='center', fontsize=16) + ax.text(1.5, 1, "Explosive\n growth", ha="center", fontsize=16) + ax.text(-1.5, 1, "Explosive\n oscillations", ha="center", fontsize=16) + ax.text(0.05, -1.5, "Explosive oscillations", ha="center", fontsize=16) + ax.text(0.09, -0.5, "Damped oscillations", ha="center", fontsize=16) # Add small marker to y-axis - ax.axhline(y=1.005, xmin=0.495, xmax=0.505, c='black') - ax.text(-0.12, -1.12, '-1', fontsize=10) - ax.text(-0.12, 0.98, '1', fontsize=10) + ax.axhline(y=1.005, xmin=0.495, xmax=0.505, c="black") + ax.text(-0.12, -1.12, "-1", fontsize=10) + ax.text(-0.12, 0.98, "1", fontsize=10) return fig + param_plot() plt.show() ``` @@ -497,28 +533,31 @@ point implied by the setting of $(\alpha,\beta)$. ```{code-cell} ipython3 def categorize_solution(ρ1, ρ2): - """This function takes values of ρ1 and ρ2 and uses them to classify the type of solution """ - discriminant = ρ1 ** 2 + 4 * ρ2 + discriminant = ρ1**2 + 4 * ρ2 if ρ2 > 1 + ρ1 or ρ2 < -1: - print('Explosive oscillations') + print("Explosive oscillations") elif ρ1 + ρ2 > 1: - print('Explosive growth') + print("Explosive growth") elif discriminant < 0: - print('Roots are complex with modulus less than one; \ -therefore damped oscillations') + print( + "Roots are complex with modulus less than one; \ +therefore damped oscillations" + ) else: - print('Roots are real and absolute values are less than one; \ -therefore get smooth convergence to a steady state') + print( + "Roots are real and absolute values are less than one; \ +therefore get smooth convergence to a steady state" + ) ``` To test the categorize_solution function, ```{code-cell} ipython3 -categorize_solution(1.3, -.4) +categorize_solution(1.3, -0.4) ``` ### Function for plotting paths @@ -527,13 +566,12 @@ A useful function for our work below is ```{code-cell} ipython3 def plot_y(function=None): - """Function plots path of Y_t""" plt.subplots(figsize=(10, 6)) plt.plot(function) - plt.xlabel('Time $t$') - plt.ylabel('$Y_t$', rotation=0) + plt.xlabel("Time $t$") + plt.ylabel("$Y_t$", rotation=0) plt.grid() plt.show() ``` @@ -551,8 +589,8 @@ that we set ```{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, α=0.92, β=0.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 are less than unity in absolute value.It also computes a simulation of @@ -565,32 +603,33 @@ def y_nonstochastic(y_0=100, y_1=80, α=.92, β=.5, γ=10, n=80): ρ1 = α + β ρ2 = -β - print(f'ρ_1 is {ρ1:.2f}') - print(f'ρ_2 is {ρ2:.2f}') + print(f"ρ_1 is {ρ1:.2f}") + print(f"ρ_2 is {ρ2:.2f}") - discriminant = ρ1 ** 2 + 4 * ρ2 + discriminant = ρ1**2 + 4 * ρ2 if discriminant == 0: roots.append(-ρ1 / 2) - print('Single real root: ') - print(''.join(f"{r:.2f}" for r in roots)) + print("Single real root: ") + 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(f"{r:.2f}" for r in roots)) + print("Two real 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(f"{r.real:.2f}{r.imag:+.2f}j" for r in roots)) + print("Two complex 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') + print("Absolute values of roots are less than one") else: - print('Absolute values of roots are not less than one') + print("Absolute values of roots are not less than one") - def transition(x, t): return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + def transition(x, t): + return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ y_t = [y_0, y_1] @@ -599,6 +638,7 @@ def y_nonstochastic(y_0=100, y_1=80, α=.92, β=.5, γ=10, n=80): return y_t + plot_y(y_nonstochastic()) ``` @@ -628,9 +668,9 @@ def f(r, ϕ): """ g1 = cmath.rect(r, ϕ) # Generate two complex roots g2 = cmath.rect(r, -ϕ) - ρ1 = g1 + g2 # Implied ρ1, ρ2 + ρ1 = g1 + g2 # Implied ρ1, ρ2 ρ2 = -g1 * g2 - β = -ρ2 # Reverse-engineer a and b that validate these + β = -ρ2 # Reverse-engineer a and b that validate these α = ρ1 - β return ρ1, ρ2, α, β ``` @@ -638,9 +678,9 @@ def f(r, ϕ): 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 -ϕ = 2 * math.pi/period +r = 0.95 +period = 10 # Length of cycle in units of time +ϕ = 2 * math.pi / period ## Apply the function ρ1, ρ2, α, β = f(r, ϕ) @@ -681,9 +721,8 @@ print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ##=== This method uses numpy to calculate roots ===# -def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): - - """ Rather than computing the roots of the characteristic +def y_nonstochastic(y_0=100, y_1=80, α=0.9, β=0.8, γ=10, n=80): + """Rather than computing the roots of the characteristic polynomial by hand as we did earlier, this function enlists numpy to do the work for us """ @@ -696,22 +735,23 @@ def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): # Find roots of polynomial roots = np.roots([1, -ρ1, -ρ2]) - print(f'Roots are {roots}') + print(f"Roots are {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") # Define transition equation - def transition(x, t): return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + def transition(x, t): + return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ # Set initial conditions y_t = [y_0, y_1] @@ -722,6 +762,7 @@ def y_nonstochastic(y_0=100, y_1=80, α=.9, β=.8, γ=10, n=80): return y_t + plot_y(y_nonstochastic()) ``` @@ -733,10 +774,10 @@ roots. We'll generate an *undamped* cycle of period 10 ```{code-cell} ipython3 -r = 1 # Generates undamped, nonexplosive cycles +r = 1 # Generates undamped, nonexplosive cycles -period = 10 # Length of cycle in units of time -ϕ = 2 * math.pi/period +period = 10 # Length of cycle in units of time +ϕ = 2 * math.pi / period ## Apply the reverse-engineering function f ρ1, ρ2, α, β = f(r, ϕ) @@ -762,7 +803,7 @@ r1 = Symbol("ρ_1") r2 = Symbol("ρ_2") z = Symbol("z") -sympy.solve(z**2 - r1*z - r2, z) +sympy.solve(z**2 - r1 * z - r2, z) ``` ```{code-cell} ipython3 @@ -771,7 +812,7 @@ sympy.solve(z**2 - r1*z - r2, z) r1 = α + β r2 = -β -sympy.solve(z**2 - r1*z - r2, z) +sympy.solve(z**2 - r1 * z - r2, z) ``` ## Stochastic shocks @@ -782,7 +823,6 @@ demand ```{code-cell} ipython3 def y_stochastic(y_0=0, y_1=0, α=0.8, β=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. @@ -815,7 +855,7 @@ def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5): ϵ = np.random.normal(0, 1, n) # Define transition equation - def transition(x, t): + def transition(x, t): return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + σ * ϵ[t] # Set initial conditions @@ -827,16 +867,17 @@ def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5): return y_t + plot_y(y_stochastic()) ``` Let's do a simulation in which there are shocks and the characteristic polynomial has complex roots ```{code-cell} ipython3 -r = .97 +r = 0.97 -period = 10 # Length of cycle in units of time -ϕ = 2 * math.pi/period +period = 10 # Length of cycle in units of time +ϕ = 2 * math.pi / period # Apply the reverse-engineering function f ρ1, ρ2, α, β = f(r, ϕ) @@ -855,17 +896,9 @@ This function computes a response to either a permanent or one-off increase in government expenditures ```{code-cell} ipython3 -def y_stochastic_g(y_0=20, - y_1=20, - α=0.8, - β=0.2, - γ=10, - n=100, - σ=2, - g=0, - g_t=0, - duration='permanent'): - +def y_stochastic_g( + y_0=20, y_1=20, α=0.8, β=0.2, γ=10, n=100, σ=2, g=0, g_t=0, duration="permanent" +): """This program computes a response to a permanent increase in government expenditures that occurs at time 20 """ @@ -883,15 +916,15 @@ def y_stochastic_g(y_0=20, # 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) @@ -923,14 +956,14 @@ def y_stochastic_g(y_0=20, y_t.append(transition(y_t, t)) # Permanent government spending shock - elif duration == 'permanent': + elif duration == "permanent": if t < g_t: y_t.append(transition(y_t, t, g=0)) else: y_t.append(transition(y_t, t, g=g)) # One-off government spending shock - elif duration == 'one-off': + elif duration == "one-off": if t == g_t: y_t.append(transition(y_t, t, g=g)) else: @@ -941,13 +974,13 @@ def y_stochastic_g(y_0=20, A permanent government spending shock can be simulated as follows ```{code-cell} ipython3 -plot_y(y_stochastic_g(g=10, g_t=20, duration='permanent')) +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} ipython3 -plot_y(y_stochastic_g(g=500, g_t=50, duration='one-off')) +plot_y(y_stochastic_g(g=500, g_t=50, duration="one-off")) ``` ## Wrapping everything into a class @@ -958,8 +991,7 @@ Now we'll roll up our sleeves and write a Python class called `Samuelson` for the Samuelson model ```{code-cell} ipython3 -class Samuelson(): - +class Samuelson: """This class represents the Samuelson model, otherwise known as the multiple-accelerator model. The model combines the Keynesian multiplier with the accelerator theory of investment. @@ -996,17 +1028,9 @@ class Samuelson(): """ - def __init__(self, - y_0=100, - y_1=50, - α=1.3, - β=0.2, - γ=10, - n=100, - σ=0, - g=0, - g_t=0, - duration=None): + def __init__( + self, y_0=100, y_1=50, α=1.3, β=0.2, γ=10, n=100, σ=0, g=0, g_t=0, duration=None + ): self.y_0, self.y_1, self.α, self.β = y_0, y_1, α, β self.n, self.g, self.g_t, self.duration = n, g, g_t, duration @@ -1017,11 +1041,11 @@ class Samuelson(): def root_type(self): if all(isinstance(root, complex) for root in self.roots): - return 'Complex conjugate' + return "Complex conjugate" elif len(self.roots) > 1: - return 'Double real' + return "Double real" else: - return 'Single real' + return "Single real" def root_less_than_one(self): if all(abs(root) < 1 for root in self.roots): @@ -1029,15 +1053,15 @@ class Samuelson(): def solution_type(self): ρ1, ρ2 = self.ρ1, self.ρ2 - discriminant = ρ1 ** 2 + 4 * ρ2 + discriminant = ρ1**2 + 4 * ρ2 if ρ2 >= 1 + ρ1 or ρ2 <= -1: - return 'Explosive oscillations' + return "Explosive oscillations" elif ρ1 + ρ2 >= 1: - return 'Explosive growth' + return "Explosive growth" elif discriminant < 0: - return 'Damped oscillations' + return "Damped oscillations" else: - return 'Steady state' + return "Steady state" def _transition(self, x, t, g): @@ -1049,8 +1073,7 @@ class Samuelson(): # Stochastic else: ϵ = np.random.normal(0, 1, self.n) - return self.ρ1 * x[t - 1] + self.ρ2 * x[t - 2] + self.γ + g \ - + self.σ * ϵ[t] + return self.ρ1 * x[t - 1] + self.ρ2 * x[t - 2] + self.γ + g + self.σ * ϵ[t] def generate_series(self): @@ -1069,14 +1092,14 @@ class Samuelson(): y_t.append(self._transition(y_t, t)) # Permanent government spending shock - elif self.duration == 'permanent': + elif self.duration == "permanent": if t < self.g_t: y_t.append(self._transition(y_t, t, g=0)) else: y_t.append(self._transition(y_t, t, g=self.g)) # One-off government spending shock - elif self.duration == 'one-off': + elif self.duration == "one-off": if t == self.g_t: y_t.append(self._transition(y_t, t, g=self.g)) else: @@ -1084,42 +1107,52 @@ class Samuelson(): return y_t def summary(self): - print('Summary\n' + '-' * 50) - print(f'Root type: {self.root_type()}') - print(f'Solution type: {self.solution_type()}') - print(f'Roots: {str(self.roots)}') + print("Summary\n" + "-" * 50) + print(f"Root type: {self.root_type()}") + print(f"Solution type: {self.solution_type()}") + print(f"Roots: {str(self.roots)}") if self.root_less_than_one() == True: - print('Absolute value of roots is less than one') + print("Absolute value of roots is less than one") else: - print('Absolute value of roots is not less than one') + print("Absolute value of roots is not less than one") if self.σ > 0: - print('Stochastic series with σ = ' + str(self.σ)) + print("Stochastic series with σ = " + str(self.σ)) else: - print('Non-stochastic series') + print("Non-stochastic series") if self.g != 0: - print('Government spending equal to ' + str(self.g)) + print("Government spending equal to " + str(self.g)) if self.duration != None: - print(self.duration.capitalize() + - ' government spending shock at t = ' + str(self.g_t)) + print( + self.duration.capitalize() + + " government spending shock at t = " + + str(self.g_t) + ) def plot(self): fig, ax = plt.subplots(figsize=(10, 6)) ax.plot(self.generate_series()) - ax.set(xlabel='Iteration', xlim=(0, self.n)) - ax.set_ylabel('$Y_t$', rotation=0) + ax.set(xlabel="Iteration", xlim=(0, self.n)) + ax.set_ylabel("$Y_t$", rotation=0) ax.grid() # Add parameter values to plot - paramstr = f'$α={self.α:.2f}$ \n $β={self.β:.2f}$ \n \ + paramstr = f"$α={self.α:.2f}$ \n $β={self.β:.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) - ax.text(0.87, 0.05, paramstr, transform=ax.transAxes, - fontsize=12, bbox=props, va='bottom') + $\\rho_1={self.ρ1:.2f}$ \n $\\rho_2={self.ρ2:.2f}$" + props = dict(fc="white", pad=10, alpha=0.5) + ax.text( + 0.87, + 0.05, + paramstr, + transform=ax.transAxes, + fontsize=12, + bbox=props, + va="bottom", + ) return fig @@ -1135,15 +1168,22 @@ class Samuelson(): for i, root in enumerate(self.roots): if isinstance(root, complex): # Need to fill operator for positive as string is split apart - operator = ['+', ''] - label = rf'$\lambda_{i+1} = {sam.roots[i].real:.2f} {operator[i]} {sam.roots[i].imag:.2f}i$' + operator = ["+", ""] + label = rf"$\lambda_{i+1} = {sam.roots[i].real:.2f} {operator[i]} {sam.roots[i].imag:.2f}i$" else: - label = rf'$\lambda_{i+1} = {sam.roots[i].real:.2f}$' - ax.scatter(0, 0, s=0, label=label) # dummy to add to legend + label = rf"$\lambda_{i+1} = {sam.roots[i].real:.2f}$" + ax.scatter(0, 0, s=0, label=label) # dummy to add to legend # Add ρ pair to plot - ax.scatter(self.ρ1, self.ρ2, s=100, c='red', marker='+', - label=r'$(\ \rho_1, \ \rho_2 \ )$', zorder=5) + ax.scatter( + self.ρ1, + self.ρ2, + s=100, + c="red", + marker="+", + label=r"$(\ \rho_1, \ \rho_2 \ )$", + zorder=5, + ) plt.legend(fontsize=12, loc=3) @@ -1155,7 +1195,7 @@ class Samuelson(): Now we'll put our Samuelson class to work on an example ```{code-cell} ipython3 -sam = Samuelson(α=0.8, β=0.5, σ=2, g=10, g_t=20, duration='permanent') +sam = Samuelson(α=0.8, β=0.5, σ=2, g=10, g_t=20, duration="permanent") sam.summary() ``` @@ -1189,6 +1229,7 @@ Here is how we map the Samuelson model into an instance of a """This script maps the Samuelson model in the the ``LinearStateSpace`` class """ + α = 0.8 β = 0.9 ρ1 = α + β @@ -1198,31 +1239,31 @@ Here is how we map the Samuelson model into an instance of a g = 10 n = 100 -A = [[1, 0, 0], - [γ + g, ρ1, ρ2], - [0, 1, 0]] +A = [[1, 0, 0], [γ + g, ρ1, ρ2], [0, 1, 0]] -G = [[γ + g, ρ1, ρ2], # this is Y_{t+1} - [γ, α, 0], # this is C_{t+1} - [0, β, -β]] # this is I_{t+1} +G = [ + [γ + g, ρ1, ρ2], # this is Y_{t+1} + [γ, α, 0], # this is C_{t+1} + [0, β, -β], +] # this is I_{t+1} μ_0 = [1, 100, 50] -C = np.zeros((3,1)) -C[1] = σ # stochastic +C = np.zeros((3, 1)) +C[1] = σ # stochastic sam_t = LinearStateSpace(A, C, G, mu_0=μ_0) x, y = sam_t.simulate(ts_length=n) fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8)) -titles = ['Output ($Y_t$)', 'Consumption ($C_t$)', 'Investment ($I_t$)'] -colors = ['darkblue', 'red', 'purple'] +titles = ["Output ($Y_t$)", "Consumption ($C_t$)", "Investment ($I_t$)"] +colors = ["darkblue", "red", "purple"] for ax, series, title, color in zip(axes, y, titles, colors): ax.plot(series, color=color) ax.set(title=title, xlim=(0, n)) ax.grid() -axes[-1].set_xlabel('Iteration') +axes[-1].set_xlabel("Iteration") plt.show() ``` @@ -1256,19 +1297,12 @@ methods and attributes) to add more functions to use ```{code-cell} ipython3 class SamuelsonLSS(LinearStateSpace): - """ This subclass creates a Samuelson multiplier-accelerator model as a linear state space system. """ - def __init__(self, - y_0=100, - y_1=50, - α=0.8, - β=0.9, - γ=10, - σ=1, - g=10): + + def __init__(self, y_0=100, y_1=50, α=0.8, β=0.9, γ=10, σ=1, g=10): self.α, self.β = α, β self.y_0, self.y_1, self.g = y_0, y_1, g @@ -1281,14 +1315,14 @@ class SamuelsonLSS(LinearStateSpace): self.ρ2 = -β # Define transition matrix - self.A = [[1, 0, 0], - [γ + g, self.ρ1, self.ρ2], - [0, 1, 0]] + self.A = [[1, 0, 0], [γ + g, self.ρ1, self.ρ2], [0, 1, 0]] # 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} + self.G = [ + [γ + g, self.ρ1, self.ρ2], # this is Y_{t+1} + [γ, α, 0], # this is C_{t+1} + [0, β, -β], + ] # this is I_{t+1} self.C = np.zeros((3, 1)) self.C[1] = σ # stochastic @@ -1306,26 +1340,27 @@ class SamuelsonLSS(LinearStateSpace): # values for simulation if stationary == True: try: - self.mu_x, self.mu_y, self.Sigma_x, self.Sigma_y, self.Sigma_yx = \ + self.mu_x, self.mu_y, self.Sigma_x, self.Sigma_y, self.Sigma_yx = ( self.stationary_distributions() + ) self.mu_0 = self.mu_x self.Sigma_0 = self.Sigma_x # Exception where no convergence achieved when - #calculating stationary distributions + # calculating stationary distributions except ValueError: - print('Stationary distribution does not exist') + print("Stationary distribution does not exist") x, y = self.simulate(ts_length) fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8)) - titles = ['Output ($Y_t$)', 'Consumption ($C_t$)', 'Investment ($I_t$)'] - colors = ['darkblue', 'red', 'purple'] + titles = ["Output ($Y_t$)", "Consumption ($C_t$)", "Investment ($I_t$)"] + colors = ["darkblue", "red", "purple"] for ax, series, title, color in zip(axes, y, titles, colors): ax.plot(series, color=color) ax.set(title=title, xlim=(0, n)) ax.grid() - axes[-1].set_xlabel('Iteration') + axes[-1].set_xlabel("Iteration") # Reset distribution parameters to their initial values self.mu_0 = temp_mu @@ -1338,25 +1373,25 @@ class SamuelsonLSS(LinearStateSpace): x, y = self.impulse_response(j) # Reshape into 3 x j matrix for plotting purposes - yimf = np.array(y).flatten().reshape(j+1, 3).T + yimf = np.array(y).flatten().reshape(j + 1, 3).T fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8)) - labels = ['$Y_t$', '$C_t$', '$I_t$'] - colors = ['darkblue', 'red', 'purple'] + labels = ["$Y_t$", "$C_t$", "$I_t$"] + colors = ["darkblue", "red", "purple"] for ax, series, label, color in zip(axes, yimf, labels, colors): ax.plot(series, color=color) ax.set(xlim=(0, j)) ax.set_ylabel(label, rotation=0, fontsize=14, labelpad=10) ax.grid() - axes[0].set_title('Impulse Response Functions') - axes[-1].set_xlabel('Iteration') + axes[0].set_title("Impulse Response Functions") + axes[-1].set_xlabel("Iteration") return fig def multipliers(self, j=5): x, y = self.impulse_response(j) - return np.sum(np.array(y).flatten().reshape(j+1, 3), axis=0) + return np.sum(np.array(y).flatten().reshape(j + 1, 3), axis=0) ``` ### Illustrations From 5c8f323ef9cc2043b69aada6b5a52a9997af36e6 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 15 Sep 2025 12:24:31 +1000 Subject: [PATCH 15/17] rewrite functions for clarity and remove duplicating code and figures --- lectures/samuelson.md | 660 +++++++++++++++++++----------------------- 1 file changed, 300 insertions(+), 360 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 842109c24..8884f4f53 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -533,31 +533,113 @@ point implied by the setting of $(\alpha,\beta)$. ```{code-cell} ipython3 def categorize_solution(ρ1, ρ2): - """This function takes values of ρ1 and ρ2 and uses them - to classify the type of solution """ - + This function takes values of ρ1 and ρ2 and uses them + to classify the type of solution. + """ discriminant = ρ1**2 + 4 * ρ2 if ρ2 > 1 + ρ1 or ρ2 < -1: - print("Explosive oscillations") + return "Explosive oscillations" elif ρ1 + ρ2 > 1: - print("Explosive growth") + return "Explosive growth" elif discriminant < 0: - print( - "Roots are complex with modulus less than one; \ -therefore damped oscillations" - ) + return "Damped oscillations" else: - print( - "Roots are real and absolute values are less than one; \ -therefore get smooth convergence to a steady state" - ) + return "Steady state convergence" + +def analyze_roots(α, β, verbose=True): + """ + Unified function to calculate roots and analyze their properties. + """ + ρ1 = α + β + ρ2 = -β + + # Compute characteristic polynomial roots + roots = np.roots([1, -ρ1, -ρ2]) + + # Classify solution type + solution_type = categorize_solution(ρ1, ρ2) + + # Determine root properties + is_complex = all(isinstance(root, complex) for root in roots) + is_stable = all(abs(root) < 1 for root in roots) + + if verbose: + print(f"ρ1 = {ρ1:.2f}, ρ2 = {ρ2:.2f}") + print(f"Roots: {[f'{root:.2f}' for root in roots]}") + print(f"Root type: {'Complex' if is_complex else 'Real'}") + print(f"Stability: {'Stable' if is_stable else 'Unstable'}") + print(f"Solution type: {solution_type}") + + return { + 'roots': roots, + 'rho1': ρ1, + 'rho2': ρ2, + 'is_complex': is_complex, + 'is_stable': is_stable, + 'solution_type': solution_type + } ``` -To test the categorize_solution function, +We also write a unified simulation function that can handle both +non-stochastic and stochastic versions of the model. + +It allows for government spending paths of a few simple forms which +we specify via a dictionary `g_params` ```{code-cell} ipython3 -categorize_solution(1.3, -0.4) +def simulate_samuelson( + y_0, y_1, α, β, γ=10, n=100, σ=0, g_params=None, seed=0 +): + """ + Unified simulation function for Samuelson model. + + Parameters: + g_params: dict with keys 'g', 'g_t', 'duration' for government spending + seed: random seed for reproducible results + """ + analysis = analyze_roots(α, β, verbose=False) + ρ1, ρ2 = analysis['rho1'], analysis['rho2'] + + # Initialize time series + y_t = [y_0, y_1] + + # Generate shocks if stochastic + if σ > 0: + np.random.seed(seed) + ϵ = np.random.normal(0, 1, n) + + # Simulate forward + for t in range(2, n): + + # Determine government spending + g = 0 + if g_params: + g_val, g_t_val = g_params.get('g', 0), g_params.get('g_t', 0) + duration = g_params.get('duration', None) + if duration == 'permanent' and t >= g_t_val: + g = g_val + elif duration == 'one-off' and t == g_t_val: + g = g_val + elif duration is None: + g = g_val + + # Calculate next value + y_next = ρ1 * y_t[t-1] + ρ2 * y_t[t-2] + γ + g + if σ > 0: + y_next += σ * ϵ[t] + + y_t.append(y_next) + + return y_t, analysis +``` + +We will use this function to run simulations of the model. + +But before doing that, let's test the analysis function + +```{code-cell} ipython3 +analysis = analyze_roots(α=1.3, β=0.4) ``` ### Function for plotting paths @@ -570,9 +652,8 @@ def plot_y(function=None): plt.subplots(figsize=(10, 6)) plt.plot(function) - plt.xlabel("Time $t$") - plt.ylabel("$Y_t$", rotation=0) - plt.grid() + plt.xlabel("$t$") + plt.ylabel("$Y_t$") plt.show() ``` @@ -587,56 +668,34 @@ The function also plots a $Y_t$ starting from initial conditions that we set ```{code-cell} ipython3 -# This is a 'manual' method - - def y_nonstochastic(y_0=100, y_1=80, α=0.92, β=0.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 - are less than unity in absolute value.It also computes a simulation of - length n starting from the two given initial conditions for national - income + """ + This function calculates the roots of the characteristic polynomial + by hand and returns a path of y_t starting from initial conditions """ - roots = [] - - ρ1 = α + β - ρ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) + y_series, analysis = simulate_samuelson(y_0, y_1, α, β, γ, n, 0, None, 42) + + print(f"ρ_1 is {analysis['rho1']:.2f}") + print(f"ρ_2 is {analysis['rho2']:.2f}") + + roots = analysis['roots'] + if len(roots) == 1: print("Single real root: ") - 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(f"{r:.2f}" for r in roots)) - else: - roots.append((-ρ1 + sqrt(discriminant)) / 2) - roots.append((-ρ1 - sqrt(discriminant)) / 2) + print(f"{roots[0]:.2f}") + elif analysis['is_complex']: print("Two complex roots: ") print(" ".join(f"{r.real:.2f}{r.imag:+.2f}j" for r in roots)) - - if all(abs(root) < 1 for root in roots): + else: + print("Two real roots: ") + print(" ".join(f"{r:.2f}" for r in roots)) + + if analysis['is_stable']: print("Absolute values of roots are less than one") else: print("Absolute values of roots are not less than one") - - def transition(x, t): - return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ - - y_t = [y_0, y_1] - - for t in range(2, n): - y_t.append(transition(y_t, t)) - - return y_t + + return y_series plot_y(y_nonstochastic()) @@ -663,39 +722,46 @@ def f(r, ϕ): and creates ρ1 and ρ2 of characteristic polynomial for which r exp(j ϕ) and r exp(- j ϕ) are complex roots. - Returns the multiplier coefficient $\alpha$ and the accelerator coefficient $\beta$ + Returns the multiplier coefficient $\alpha$ + and the accelerator coefficient $\beta$ that verifies those roots. """ - g1 = cmath.rect(r, ϕ) # Generate two complex roots + # Create complex conjugate pair from polar coordinates + g1 = cmath.rect(r, ϕ) g2 = cmath.rect(r, -ϕ) - ρ1 = g1 + g2 # Implied ρ1, ρ2 + + # Calculate corresponding ρ1, ρ2 parameters + ρ1 = g1 + g2 ρ2 = -g1 * g2 - β = -ρ2 # Reverse-engineer a and b that validate these + + # Derive α and β coefficients from ρ parameters + β = -ρ2 α = ρ1 - β return ρ1, ρ2, α, β ``` -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 = 0.95 -period = 10 # Length of cycle in units of time + +# Cycle period in time units +period = 10 ϕ = 2 * math.pi / period -## Apply the function +# Apply the reverse-engineering function ρ1, ρ2, α, β = f(r, ϕ) print(f"α, β = {α:.2f}, {β:.2f}") print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ``` -```{code-cell} ipython3 -## Print the real components of ρ1 and ρ2 - -ρ1 = ρ1.real -ρ2 = ρ2.real +The real parts of the roots are -print(f"ρ1 = {ρ1:.2f}, ρ2 = {ρ2:.2f}") +```{code-cell} ipython3 +print(f"ρ1 = {ρ1.real:.2f}, ρ2 = {ρ2.real:.2f}") ``` ### Root finding using numpy @@ -711,56 +777,26 @@ p2 = cmath.polar(r2) 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"α, β = {α:.2f}, {β:.2f}") print(f"ρ1, ρ2 = {ρ1:.2f}, {ρ2:.2f}") ``` ```{code-cell} ipython3 -##=== This method uses numpy to calculate roots ===# - - def y_nonstochastic(y_0=100, y_1=80, α=0.9, β=0.8, γ=10, n=80): - """Rather than computing the roots of the characteristic - polynomial by hand as we did earlier, this function - enlists numpy to do the work for us + """ + This function enlists numpy to calculate the roots of the characteristic + polynomial. """ - # Useful constants - ρ1 = α + β - ρ2 = -β - - categorize_solution(ρ1, ρ2) - - # Find roots of polynomial - roots = np.roots([1, -ρ1, -ρ2]) - print(f"Roots are {roots}") - - # Check if real or complex - if all(isinstance(root, complex) for root in roots): - print("Roots are complex") - else: - 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") - else: - print("Roots are not less than one") - - # Define transition equation - def transition(x, t): - return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ - - # Set initial conditions - y_t = [y_0, y_1] - - # Generate y_t series - for t in range(2, n): - y_t.append(transition(y_t, t)) - - return y_t + y_series, analysis = simulate_samuelson(y_0, y_1, α, β, γ, n, 0, None, 42) + + print(f"Solution type: {analysis['solution_type']}") + print(f"Roots are {analysis['roots']}") + print(f"Root type: {'Complex' if analysis['is_complex'] else 'Real'}") + print(f"Stability: {'Stable' if analysis['is_stable'] else 'Unstable'}") + + return y_series plot_y(y_nonstochastic()) @@ -779,14 +815,14 @@ 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 +# Apply the reverse-engineering function f ρ1, ρ2, α, β = f(r, ϕ) -# Drop the imaginary part so that it is a valid input into y_nonstochastic +# Extract real parts for numerical computation α = α.real β = β.real -print(f"a, b = {α:.2f}, {β:.2f}") +print(f"α, β = {α:.2f}, {β:.2f}") ytemp = y_nonstochastic(α=α, β=β, y_0=20, y_1=30) plot_y(ytemp) @@ -823,55 +859,27 @@ demand ```{code-cell} ipython3 def y_stochastic(y_0=0, y_1=0, α=0.8, β=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. + """ + This function takes parameters of a stochastic version of + the model, analyzes the roots of the characteristic + polynomial and generates a simulation. """ - # Useful constants - ρ1 = α + β - ρ2 = -β - - # Categorize solution - categorize_solution(ρ1, ρ2) - - # Find roots of polynomial - roots = np.roots([1, -ρ1, -ρ2]) - 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") - else: - 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") - else: - 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] - - # Set initial conditions - y_t = [y_0, y_1] - - # Generate y_t series - for t in range(2, n): - y_t.append(transition(y_t, t)) - - return y_t + y_series, analysis = simulate_samuelson(y_0, y_1, α, β, γ, n, σ, None, 42) + + print(f"Solution type: {analysis['solution_type']}") + print(f"Roots are {[f'{root:.2f}' for root in analysis['roots']]}") + print(f"Root type: {'Complex' if analysis['is_complex'] else 'Real'}") + print(f"Stability: {'Stable' if analysis['is_stable'] else 'Unstable'}") + + return y_series plot_y(y_stochastic()) ``` -Let's do a simulation in which there are shocks and the characteristic polynomial has complex roots +Let's do a simulation in which there are shocks and the characteristic +polynomial has complex roots ```{code-cell} ipython3 r = 0.97 @@ -882,11 +890,11 @@ period = 10 # Length of cycle in units of time # Apply the reverse-engineering function f ρ1, ρ2, α, β = f(r, ϕ) -# Drop the imaginary part so that it is a valid input into y_nonstochastic +# Extract real parts for numerical computation α = α.real β = β.real -print(f"a, b = {α:.2f}, {β:.2f}") +print(f"α, β = {α:.2f}, {β:.2f}") plot_y(y_stochastic(y_0=40, y_1=42, α=α, β=β, σ=2, n=100)) ``` @@ -897,78 +905,27 @@ in government expenditures ```{code-cell} ipython3 def y_stochastic_g( - y_0=20, y_1=20, α=0.8, β=0.2, γ=10, n=100, σ=2, g=0, g_t=0, duration="permanent" + y_0=20, y_1=20, α=0.8, β=0.2, γ=10, + n=100, σ=2, g=0, g_t=0, duration="permanent" ): - """This program computes a response to a permanent increase + """ + This program computes a response to a permanent increase in government expenditures that occurs at time 20 """ - # Useful constants - ρ1 = α + β - ρ2 = -β - - # Categorize solution - categorize_solution(ρ1, ρ2) - - # Find roots of polynomial - roots = np.roots([1, -ρ1, -ρ2]) - print(roots) - - # Check if real or complex - if all(isinstance(root, complex) for root in roots): - print("Roots are complex") - else: - 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") - else: - print("Roots are not less than one") - - # Generate shocks - ϵ = np.random.normal(0, 1, n) - - def transition(x, t, g): - - # Non-stochastic - separated to avoid generating random series - # when not needed - if σ == 0: - return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + g - - # Stochastic - else: - ϵ = np.random.normal(0, 1, n) - return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + g + σ * ϵ[t] - - # Create list and set initial conditions - y_t = [y_0, y_1] - - # Generate y_t series - for t in range(2, n): - - # No government spending - if g == 0: - y_t.append(transition(y_t, t)) - - # Government spending (no shock) - elif g != 0 and duration == None: - y_t.append(transition(y_t, t)) - - # Permanent government spending shock - elif duration == "permanent": - if t < g_t: - y_t.append(transition(y_t, t, g=0)) - else: - y_t.append(transition(y_t, t, g=g)) - - # One-off government spending shock - elif duration == "one-off": - if t == g_t: - y_t.append(transition(y_t, t, g=g)) - else: - y_t.append(transition(y_t, t, g=0)) - return y_t + g_params = ( + {'g': g, 'g_t': g_t, 'duration': duration} if g != 0 else None + ) + y_series, analysis = simulate_samuelson( + y_0, y_1, α, β, γ, n, σ, g_params, 42 + ) + + print(f"Solution type: {analysis['solution_type']}") + print(f"Roots: {analysis['roots']}") + print(f"Root type: {'Complex' if analysis['is_complex'] else 'Real'}") + print(f"Stability: {'Stable' if analysis['is_stable'] else 'Unstable'}") + + return y_series ``` A permanent government spending shock can be simulated as follows @@ -992,15 +949,18 @@ for the Samuelson model ```{code-cell} ipython3 class Samuelson: - """This class represents the Samuelson model, otherwise known as the - multiple-accelerator model. The model combines the Keynesian multiplier + """ + This class represents the Samuelson model, otherwise known as the + multiplier-accelerator model. + The model combines the Keynesian multiplier with the accelerator theory of investment. - The path of output is governed by a linear second-order difference equation + The path of output is governed by a linear + second-order difference equation .. math:: - Y_t = $\alpha$ (1 + $\beta$) Y_{t-1} - $\alpha$ $\beta$ Y_{t-2} + Y_t = \alpha (1 + \beta) Y_{t-1} - \alpha \beta Y_{t-2} Parameters ---------- @@ -1029,82 +989,38 @@ class Samuelson: """ def __init__( - self, y_0=100, y_1=50, α=1.3, β=0.2, γ=10, n=100, σ=0, g=0, g_t=0, duration=None + self, y_0=100, y_1=50, + α=1.3, β=0.2, γ=10, n=100, σ=0, g=0, g_t=0, duration=None ): self.y_0, self.y_1, self.α, self.β = y_0, y_1, α, β self.n, self.g, self.g_t, self.duration = n, g, g_t, duration self.γ, self.σ = γ, σ - self.ρ1 = α + β - self.ρ2 = -β - self.roots = np.roots([1, -self.ρ1, -self.ρ2]) + + # Use unified analysis function + self.analysis = analyze_roots(α, β, verbose=False) + self.ρ1, self.ρ2 = self.analysis['rho1'], self.analysis['rho2'] + self.roots = self.analysis['roots'] def root_type(self): - if all(isinstance(root, complex) for root in self.roots): - return "Complex conjugate" - elif len(self.roots) > 1: - return "Double real" - else: - return "Single real" + return "Complex conjugate" if self.analysis['is_complex'] else "Real" def root_less_than_one(self): - if all(abs(root) < 1 for root in self.roots): - return True + return self.analysis['is_stable'] def solution_type(self): - ρ1, ρ2 = self.ρ1, self.ρ2 - discriminant = ρ1**2 + 4 * ρ2 - if ρ2 >= 1 + ρ1 or ρ2 <= -1: - return "Explosive oscillations" - elif ρ1 + ρ2 >= 1: - return "Explosive growth" - elif discriminant < 0: - return "Damped oscillations" - else: - return "Steady state" - - def _transition(self, x, t, g): - - # Non-stochastic - separated to avoid generating random series - # when not needed - if self.σ == 0: - return self.ρ1 * x[t - 1] + self.ρ2 * x[t - 2] + self.γ + g - - # Stochastic - else: - ϵ = np.random.normal(0, 1, self.n) - return self.ρ1 * x[t - 1] + self.ρ2 * x[t - 2] + self.γ + g + self.σ * ϵ[t] - - def generate_series(self): - - # Create list and set initial conditions - y_t = [self.y_0, self.y_1] - - # Generate y_t series - for t in range(2, self.n): + return self.analysis['solution_type'] - # No government spending - if self.g == 0: - y_t.append(self._transition(y_t, t)) - - # Government spending (no shock) - elif self.g != 0 and self.duration == None: - y_t.append(self._transition(y_t, t)) - - # Permanent government spending shock - elif self.duration == "permanent": - if t < self.g_t: - y_t.append(self._transition(y_t, t, g=0)) - else: - y_t.append(self._transition(y_t, t, g=self.g)) - - # One-off government spending shock - elif self.duration == "one-off": - if t == self.g_t: - y_t.append(self._transition(y_t, t, g=self.g)) - else: - y_t.append(self._transition(y_t, t, g=0)) - return y_t + def generate_series(self, seed=0): + g_params = ( + {'g': self.g, 'g_t': self.g_t, 'duration': self.duration} + if self.g != 0 else None + ) + y_series, _ = simulate_samuelson( + self.y_0, self.y_1, self.α, self.β, self.γ, + self.n, self.σ, g_params, seed + ) + return y_series def summary(self): print("Summary\n" + "-" * 50) @@ -1132,17 +1048,18 @@ class Samuelson: + str(self.g_t) ) - def plot(self): + def plot(self, seed=0): fig, ax = plt.subplots(figsize=(10, 6)) - ax.plot(self.generate_series()) - ax.set(xlabel="Iteration", xlim=(0, self.n)) + ax.plot(self.generate_series(seed)) + ax.set(xlabel="iteration", xlim=(0, self.n)) ax.set_ylabel("$Y_t$", rotation=0) - ax.grid() - # Add parameter values to plot - paramstr = f"$α={self.α:.2f}$ \n $β={self.β:.2f}$ \n \ - $\\gamma={self.γ:.2f}$ \n $\\sigma={self.σ:.2f}$ \n \ - $\\rho_1={self.ρ1:.2f}$ \n $\\rho_2={self.ρ2:.2f}$" + # Display model parameters on the plot + paramstr = ( + f"$α={self.α:.2f}$ \n $β={self.β:.2f}$ \n " + f"$\\gamma={self.γ:.2f}$ \n $\\sigma={self.σ:.2f}$ \n " + f"$\\rho_1={self.ρ1:.2f}$ \n $\\rho_2={self.ρ2:.2f}$" + ) props = dict(fc="white", pad=10, alpha=0.5) ax.text( 0.87, @@ -1158,30 +1075,38 @@ class Samuelson: def param_plot(self): - # Uses the param_plot() function defined earlier (it is then able - # to be used standalone or as part of the model) fig = param_plot() ax = fig.gca() - # Add λ values to legend + # Display eigenvalues in the legend for i, root in enumerate(self.roots): if isinstance(root, complex): - # Need to fill operator for positive as string is split apart + + # Handle sign formatting for complex number display operator = ["+", ""] - label = rf"$\lambda_{i+1} = {sam.roots[i].real:.2f} {operator[i]} {sam.roots[i].imag:.2f}i$" + root_real = self.roots[i].real + root_imag = self.roots[i].imag + label = ( + rf"$\lambda_{i+1} = {root_real:.2f}" + rf"{operator[i]} {root_imag:.2f}i$" + ) else: - label = rf"$\lambda_{i+1} = {sam.roots[i].real:.2f}$" - ax.scatter(0, 0, s=0, label=label) # dummy to add to legend + label = rf"$\lambda_{i+1} = {self.roots[i].real:.2f}$" + + # Add invisible point for legend entry + ax.scatter( + 0, 0, s=0, label=label + ) - # Add ρ pair to plot + # Mark current parameter values on the stability diagram ax.scatter( self.ρ1, self.ρ2, s=100, c="red", marker="+", - label=r"$(\ \rho_1, \ \rho_2 \ )$", + label=r"$(\rho_1, \rho_2)$", zorder=5, ) @@ -1242,14 +1167,14 @@ n = 100 A = [[1, 0, 0], [γ + g, ρ1, ρ2], [0, 1, 0]] G = [ - [γ + g, ρ1, ρ2], # this is Y_{t+1} - [γ, α, 0], # this is C_{t+1} - [0, β, -β], -] # this is I_{t+1} + [γ + g, ρ1, ρ2], # Y_{t+1} + [γ, α, 0], # C_{t+1} + [0, β, -β], # I_{t+1} +] μ_0 = [1, 100, 50] C = np.zeros((3, 1)) -C[1] = σ # stochastic +C[1] = σ # Shock variance sam_t = LinearStateSpace(A, C, G, mu_0=μ_0) @@ -1261,9 +1186,8 @@ colors = ["darkblue", "red", "purple"] for ax, series, title, color in zip(axes, y, titles, colors): ax.plot(series, color=color) ax.set(title=title, xlim=(0, n)) - ax.grid() -axes[-1].set_xlabel("Iteration") +axes[-1].set_xlabel("iteration") plt.show() ``` @@ -1308,71 +1232,89 @@ class SamuelsonLSS(LinearStateSpace): 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] + # Set initial state vector + self.initial_μ = [1, y_0, y_1] self.ρ1 = α + β self.ρ2 = -β - # Define transition matrix + # Construct state transition matrix self.A = [[1, 0, 0], [γ + g, self.ρ1, self.ρ2], [0, 1, 0]] - # Define output matrix + # Construct observation 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} + [γ + g, self.ρ1, self.ρ2], # Y_{t+1} + [γ, α, 0], # C_{t+1} + [0, β, -β], # I_{t+1} + ] self.C = np.zeros((3, 1)) - self.C[1] = σ # stochastic - - # Initialize LSS with parameters from Samuelson model - LinearStateSpace.__init__(self, self.A, self.C, self.G, mu_0=self.μ_0) - - def plot_simulation(self, ts_length=100, stationary=True): + self.C[1] = σ # Shock variance - # Temporarily store original parameters - temp_mu = self.mu_0 - temp_Sigma = self.Sigma_0 + # Initialize the LinearStateSpace instance + LinearStateSpace.__init__( + self, self.A, self.C, self.G, mu_0=self.initial_μ + ) - # Set distribution parameters equal to their stationary - # values for simulation + # Create unicode aliases for mu_0 and Sigma_0 in the parent class + @property + def μ_0(self): + return self.mu_0 + + @μ_0.setter + def μ_0(self, value): + self.mu_0 = value + + @property + def Σ_0(self): + return self.Sigma_0 + + @Σ_0.setter + def Σ_0(self, value): + self.Sigma_0 = value + + def plot_simulation(self, ts_length=100, stationary=True, seed=0): + + # Store original distribution parameters + temp_μ = self.μ_0 + temp_Σ = self.Σ_0 + + # Use stationary distribution for simulation if stationary == True: try: - self.mu_x, self.mu_y, self.Sigma_x, self.Sigma_y, self.Sigma_yx = ( - self.stationary_distributions() - ) - self.mu_0 = self.mu_x - self.Sigma_0 = self.Sigma_x - # Exception where no convergence achieved when - # calculating stationary distributions + (self.μ_x, self.μ_y, self.Σ_x, self.Σ_y, self.Σ_yx + ) = self.stationary_distributions() + self.μ_0 = self.μ_x + self.Σ_0 = self.Σ_x + + # Handle case where stationary distribution doesn't exist except ValueError: print("Stationary distribution does not exist") + np.random.seed(seed) x, y = self.simulate(ts_length) fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8)) - titles = ["Output ($Y_t$)", "Consumption ($C_t$)", "Investment ($I_t$)"] + titles = ["Output ($Y_t$)", + "Consumption ($C_t$)", + "Investment ($I_t$)"] colors = ["darkblue", "red", "purple"] for ax, series, title, color in zip(axes, y, titles, colors): ax.plot(series, color=color) ax.set(title=title, xlim=(0, n)) - ax.grid() - axes[-1].set_xlabel("Iteration") + axes[-1].set_xlabel("iteration") + plt.show() - # Reset distribution parameters to their initial values - self.mu_0 = temp_mu - self.Sigma_0 = temp_Sigma - - return fig + # Restore original distribution parameters + self.μ_0 = temp_μ + self.Σ_0 = temp_Σ def plot_irf(self, j=5): x, y = self.impulse_response(j) - # Reshape into 3 x j matrix for plotting purposes + # Reshape impulse responses for plotting yimf = np.array(y).flatten().reshape(j + 1, 3).T fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8)) @@ -1382,12 +1324,10 @@ class SamuelsonLSS(LinearStateSpace): ax.plot(series, color=color) ax.set(xlim=(0, j)) ax.set_ylabel(label, rotation=0, fontsize=14, labelpad=10) - ax.grid() - axes[0].set_title("Impulse Response Functions") - axes[-1].set_xlabel("Iteration") - - return fig + axes[0].set_title("Impulse response functions") + axes[-1].set_xlabel("iteration") + plt.show() def multipliers(self, j=5): x, y = self.impulse_response(j) From d4ad185bb65f7ae08598ee38144472a54f503610 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 15 Sep 2025 12:29:00 +1000 Subject: [PATCH 16/17] minor updates --- lectures/samuelson.md | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 8884f4f53..36845ff90 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -1151,10 +1151,6 @@ Here is how we map the Samuelson model into an instance of a `LinearStateSpace` class ```{code-cell} ipython3 -"""This script maps the Samuelson model in the the -``LinearStateSpace`` class -""" - α = 0.8 β = 0.9 ρ1 = α + β @@ -1286,7 +1282,7 @@ class SamuelsonLSS(LinearStateSpace): ) = self.stationary_distributions() self.μ_0 = self.μ_x self.Σ_0 = self.Σ_x - + # Handle case where stationary distribution doesn't exist except ValueError: print("Stationary distribution does not exist") From c0e6bc12fca9bd2edad57a228b5a4f7c992026dc Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 15 Sep 2025 12:39:27 +1000 Subject: [PATCH 17/17] revert calculation by hand update --- lectures/samuelson.md | 50 ++++++++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 17 deletions(-) diff --git a/lectures/samuelson.md b/lectures/samuelson.md index 36845ff90..44aaf3b91 100644 --- a/lectures/samuelson.md +++ b/lectures/samuelson.md @@ -662,7 +662,7 @@ def plot_y(function=None): The following function calculates roots of the characteristic polynomial using high school algebra. -(We'll calculate the roots in other ways later) +(We'll calculate the roots in other ways later using `analyze_roots`.) The function also plots a $Y_t$ starting from initial conditions that we set @@ -673,29 +673,45 @@ def y_nonstochastic(y_0=100, y_1=80, α=0.92, β=0.5, γ=10, n=80): This function calculates the roots of the characteristic polynomial by hand and returns a path of y_t starting from initial conditions """ + roots = [] - y_series, analysis = simulate_samuelson(y_0, y_1, α, β, γ, n, 0, None, 42) - - print(f"ρ_1 is {analysis['rho1']:.2f}") - print(f"ρ_2 is {analysis['rho2']:.2f}") - - roots = analysis['roots'] - if len(roots) == 1: + ρ1 = α + β + ρ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(f"{roots[0]:.2f}") - elif analysis['is_complex']: - print("Two complex roots: ") - print(" ".join(f"{r.real:.2f}{r.imag:+.2f}j" for r in roots)) - else: + 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(f"{r:.2f}" for r in roots)) - - if analysis['is_stable']: + else: + roots.append((-ρ1 + sqrt(discriminant)) / 2) + roots.append((-ρ1 - sqrt(discriminant)) / 2) + print("Two complex 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") else: print("Absolute values of roots are not less than one") - - return y_series + + def transition(x, t): + return ρ1 * x[t - 1] + ρ2 * x[t - 2] + γ + + y_t = [y_0, y_1] + + for t in range(2, n): + y_t.append(transition(y_t, t)) + + return y_t plot_y(y_nonstochastic())