{ "cells": [ { "cell_type": "markdown", "id": "f93fcd6e", "metadata": {}, "source": [ "$$\n", "\\newcommand{\\argmax}{\\operatorname*{argmax}}\n", "\\newcommand{\\argmin}{\\operatorname*{argmin}}\n", "\\newcommand{\\EE}{\\mathbb{E}}\n", "$$" ] }, { "cell_type": "markdown", "id": "6f5c1e25", "metadata": {}, "source": [ "
\n",
" \n",
" 
\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"id": "dedacb58",
"metadata": {},
"source": [
"# Samuelson Multiplier-Accelerator"
]
},
{
"cell_type": "markdown",
"id": "13e6072d",
"metadata": {},
"source": [
"## Contents\n",
"\n",
"- [Samuelson Multiplier-Accelerator](#Samuelson-Multiplier-Accelerator) \n",
" - [Overview](#Overview) \n",
" - [Details](#Details) \n",
" - [Implementation](#Implementation) \n",
" - [Stochastic shocks](#Stochastic-shocks) \n",
" - [Government spending](#Government-spending) \n",
" - [Wrapping everything into a class](#Wrapping-everything-into-a-class) \n",
" - [Using the LinearStateSpace class](#Using-the-LinearStateSpace-class) \n",
" - [Pure multiplier model](#Pure-multiplier-model) \n",
" - [Summary](#Summary) "
]
},
{
"cell_type": "markdown",
"id": "cadf55ed",
"metadata": {},
"source": [
"In addition to what’s in Anaconda, this lecture will need the following libraries:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13a8b417",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"!pip install quantecon"
]
},
{
"cell_type": "markdown",
"id": "b20190ff",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"This lecture creates non-stochastic and stochastic versions of Paul Samuelson’s celebrated multiplier accelerator model [[Samuelson, 1939](https://python.quantecon.org/zreferences.html#id180)].\n",
"\n",
"In doing so, we extend the example of the Solow model class in [our second OOP lecture](https://python-programming.quantecon.org/python_oop.html).\n",
"\n",
"Our objectives are to\n",
"\n",
"- provide a more detailed example of OOP and classes \n",
"- review a famous model \n",
"- review linear difference equations, both deterministic and stochastic \n",
"\n",
"\n",
"Let’s start with some standard imports:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "269592f0",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"id": "8760c3bd",
"metadata": {},
"source": [
"We’ll also use the following for various tasks described below:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d75cc1c4",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"from quantecon import LinearStateSpace\n",
"import cmath\n",
"import math\n",
"import sympy\n",
"from sympy import Symbol, init_printing\n",
"from cmath import sqrt"
]
},
{
"cell_type": "markdown",
"id": "eca6197b",
"metadata": {},
"source": [
"### Samuelson’s model\n",
"\n",
"Samuelson used a *second-order linear difference equation* to\n",
"represent a model of national output based on three components:\n",
"\n",
"- a *national output identity* asserting that national output or national income is the\n",
" sum of consumption plus investment plus government purchases. \n",
"- a Keynesian *consumption function* asserting that consumption at\n",
" time $ t $ is equal to a constant times national output at time $ t-1 $. \n",
"- an investment *accelerator* asserting that investment at time\n",
" $ t $ equals a constant called the *accelerator coefficient*\n",
" times the difference in output between period $ t-1 $ and\n",
" $ t-2 $. \n",
"\n",
"\n",
"Consumption plus investment plus government purchases\n",
"constitute *aggregate demand,* which automatically calls forth an\n",
"equal amount of *aggregate supply*.\n",
"\n",
"(To read about linear difference equations see [here](https://en.wikipedia.org/wiki/Linear_difference_equation) or chapter IX of [[Sargent, 1987](https://python.quantecon.org/zreferences.html#id308)].)\n",
"\n",
"Samuelson used the model to analyze how particular values of the\n",
"marginal propensity to consume and the accelerator coefficient might\n",
"give rise to transient **business cycles** in national output.\n",
"\n",
"Possible dynamic properties include\n",
"\n",
"- smooth convergence to a constant level of output \n",
"- damped business cycles that eventually converge to a constant level\n",
" of output \n",
"- persistent business cycles that neither dampen nor explode \n",
"\n",
"\n",
"Later we present an extension that\n",
"adds a random shock to the right side of the national income\n",
"identity representing random fluctuations in aggregate demand.\n",
"\n",
"This modification makes national output become governed by a second-order\n",
"**stochastic linear difference equation** that, with appropriate parameter values,\n",
"gives rise to recurrent irregular business cycles.\n",
"\n",
"(To read about stochastic linear difference equations see chapter XI of\n",
"[[Sargent, 1987](https://python.quantecon.org/zreferences.html#id308)].)"
]
},
{
"cell_type": "markdown",
"id": "aacff4f5",
"metadata": {},
"source": [
"## Details\n",
"\n",
"Let’s assume that\n",
"\n",
"- $ \\{G_t\\} $ is a sequence of levels of government expenditures –\n",
" we’ll start by setting $ G_t = G $ for all $ t $. \n",
"- $ \\{C_t\\} $ is a sequence of levels of aggregate consumption\n",
" expenditures, a key endogenous variable in the model. \n",
"- $ \\{I_t\\} $ is a sequence of rates of investment, another key\n",
" endogenous variable. \n",
"- $ \\{Y_t\\} $ is a sequence of levels of national income, yet\n",
" another endogenous variable. \n",
"\n",
"\n",
"- $ \\alpha $ is the marginal propensity to consume in the Keynesian\n",
" consumption function $ C_t = \\alpha Y_{t-1} + \\gamma $. \n",
"- $ \\beta $ is the “accelerator coefficient†in the “investment\n",
" accelerator†$ I_t = \\beta (Y_{t-1} - Y_{t-2}) $. \n",
"- $ \\{\\epsilon_{t}\\} $ is an IID sequence standard normal random variables. \n",
"- $ \\sigma \\geq 0 $ is a “volatilityâ€\n",
" parameter — setting $ \\sigma = 0 $ recovers the non-stochastic case\n",
" that we’ll start with. \n",
"\n",
"\n",
"The model combines the consumption function\n",
"\n",
"\n",
"\n",
"$$\n",
"C_t = \\alpha Y_{t-1} + \\gamma \\tag{41.1}\n",
"$$\n",
"\n",
"with the investment accelerator\n",
"\n",
"\n",
"\n",
"$$\n",
"I_t = \\beta (Y_{t-1} - Y_{t-2}) \\tag{41.2}\n",
"$$\n",
"\n",
"and the national income identity\n",
"\n",
"\n",
"\n",
"$$\n",
"Y_t = C_t + I_t + G_t \\tag{41.3}\n",
"$$\n",
"\n",
"- The parameter $ \\alpha $ is peoples’ **marginal propensity to consume**\n",
" out of income - equation [(41.1)](#equation-consumption) asserts that people consume a fraction of\n",
" $ \\alpha \\in (0,1) $ of each additional dollar of income. \n",
"- The parameter $ \\beta > 0 $ is the investment accelerator coefficient - equation\n",
" [(41.2)](#equation-accelerator) asserts that people invest in physical capital when\n",
" income is increasing and disinvest when it is decreasing. \n",
"\n",
"\n",
"Equations [(41.1)](#equation-consumption), [(41.2)](#equation-accelerator), and [(41.3)](#equation-income-identity)\n",
"imply the following second-order linear difference equation for national income:\n",
"\n",
"$$\n",
"Y_t = (\\alpha+\\beta) Y_{t-1} - \\beta Y_{t-2} + (\\gamma + G_t)\n",
"$$\n",
"\n",
"or\n",
"\n",
"\n",
"\n",
"$$\n",
"Y_t = \\rho_1 Y_{t-1} + \\rho_2 Y_{t-2} + (\\gamma + G_t) \\tag{41.4}\n",
"$$\n",
"\n",
"where $ \\rho_1 = (\\alpha+\\beta) $ and $ \\rho_2 = -\\beta $.\n",
"\n",
"To complete the model, we require two *initial conditions*.\n",
"\n",
"If the model is to generate time series for $ t=0, \\ldots, T $, we\n",
"require initial values\n",
"\n",
"$$\n",
"Y_{-1} = \\bar Y_{-1}, \\quad Y_{-2} = \\bar Y_{-2}\n",
"$$\n",
"\n",
"We’ll ordinarily set the parameters $ (\\alpha,\\beta) $ so that starting from\n",
"an arbitrary pair of initial conditions\n",
"$ (\\bar Y_{-1}, \\bar Y_{-2}) $, national income $ Y_t $ converges to\n",
"a constant value as $ t $ becomes large.\n",
"\n",
"We are interested in studying\n",
"\n",
"- the transient fluctuations in $ Y_t $ as it converges to its\n",
" **steady state** level \n",
"- the *rate* at which it converges to a steady state level \n",
"\n",
"\n",
"The deterministic version of the model described so far — meaning that\n",
"no random shocks hit aggregate demand — has only transient fluctuations.\n",
"\n",
"We can convert the model to one that has persistent irregular\n",
"fluctuations by adding a random shock to aggregate demand."
]
},
{
"cell_type": "markdown",
"id": "cb829013",
"metadata": {},
"source": [
"### Stochastic version of the model\n",
"\n",
"We create a *random* or *stochastic* version of the model by adding\n",
"a random process of *shocks* or *disturbances*\n",
"$ \\{\\sigma \\epsilon_t \\} $ to the right side of equation [(41.4)](#equation-second-order),\n",
"leading to the *second-order scalar linear stochastic difference equation*:\n",
"\n",
"\n",
"\n",
"$$\n",
"Y_t = (\\alpha+\\beta) Y_{t-1} - \\beta Y_{t-2} + (\\gamma + G_t) + \\sigma \\epsilon_t \\tag{41.5}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "2ce9cf0e",
"metadata": {},
"source": [
"### Mathematical analysis of the model\n",
"\n",
"To get started, let’s set $ G_t \\equiv 0 $, $ \\sigma = 0 $, and\n",
"$ \\gamma = 0 $.\n",
"\n",
"Then we can write equation [(41.5)](#equation-second-stochastic) as\n",
"\n",
"$$\n",
"Y_t = \\rho_1 Y_{t-1} + \\rho_2 Y_{t-2}\n",
"$$\n",
"\n",
"or\n",
"\n",
"\n",
"\n",
"$$\n",
"Y_{t+2} - \\rho_1 Y_{t+1} - \\rho_2 Y_t = 0 \\tag{41.6}\n",
"$$\n",
"\n",
"To discover the properties of the solution of [(41.6)](#equation-second-stochastic2),\n",
"it is useful first to form the **characteristic polynomial**\n",
"for [(41.6)](#equation-second-stochastic2):\n",
"\n",
"\n",
"\n",
"$$\n",
"z^2 - \\rho_1 z - \\rho_2 \\tag{41.7}\n",
"$$\n",
"\n",
"where $ z $ is possibly a complex number.\n",
"\n",
"We want to find the two **zeros** (a.k.a. **roots**) – namely\n",
"$ \\lambda_1, \\lambda_2 $ – of the characteristic polynomial.\n",
"\n",
"These are two special values of $ z $, say $ z= \\lambda_1 $ and\n",
"$ z= \\lambda_2 $, such that if we set $ z $ equal to one of\n",
"these values in expression [(41.7)](#equation-polynomial),\n",
"the characteristic polynomial [(41.7)](#equation-polynomial) equals zero:\n",
"\n",
"\n",
"\n",
"$$\n",
"z^2 - \\rho_1 z - \\rho_2 = (z- \\lambda_1 ) (z -\\lambda_2) = 0 \\tag{41.8}\n",
"$$\n",
"\n",
"Equation [(41.8)](#equation-polynomial-sol) is said to *factor* the characteristic polynomial.\n",
"\n",
"When the roots are complex, they will occur as a complex conjugate pair.\n",
"\n",
"When the roots are complex, it is convenient to represent them in the\n",
"polar form\n",
"\n",
"$$\n",
"\\lambda_1 = r e^{i \\omega}, \\ \\lambda_2 = r e^{-i \\omega}\n",
"$$\n",
"\n",
"where $ r $ is the *amplitude* of the complex number and\n",
"$ \\omega $ is its *angle* or *phase*.\n",
"\n",
"These can also be represented as\n",
"\n",
"$$\n",
"\\lambda_1 = r (cos (\\omega) + i \\sin (\\omega))\n",
"$$\n",
"\n",
"$$\n",
"\\lambda_2 = r (cos (\\omega) - i \\sin(\\omega))\n",
"$$\n",
"\n",
"(To read about the polar form, see\n",
"[here](https://www.khanacademy.org/math/precalculus/x9e81a4f98389efdf:complex/x9e81a4f98389efdf:complex-mul-div-polar/a/complex-number-polar-form-review))\n",
"\n",
"Given *initial conditions* $ Y_{-1}, Y_{-2} $, we want to generate\n",
"a *solution* of the difference equation [(41.6)](#equation-second-stochastic2).\n",
"\n",
"It can be represented as\n",
"\n",
"$$\n",
"Y_t = \\lambda_1^t c_1 + \\lambda_2^t c_2\n",
"$$\n",
"\n",
"where $ c_1 $ and $ c_2 $ are constants that depend on the two\n",
"initial conditions and on $ \\rho_1, \\rho_2 $.\n",
"\n",
"When the roots are complex, it is useful to pursue the following calculations.\n",
"\n",
"Notice that\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" Y_t & = c_1 (r e^{i \\omega})^t + c_2 (r e^{-i \\omega})^t \\\\\n",
" & = c_1 r^t e^{i\\omega t} + c_2 r^t e^{-i \\omega t} \\\\\n",
" & = c_1 r^t [\\cos(\\omega t) + i \\sin(\\omega t) ] + c_2 r^t [\\cos(\\omega t) - i \\sin(\\omega t) ] \\\\\n",
" & = (c_1 + c_2) r^t \\cos(\\omega t) + i (c_1 - c_2) r^t \\sin(\\omega t)\n",
" \\end{aligned}\n",
"$$\n",
"\n",
"The only way that $ Y_t $ can be a real number for each $ t $ is if $ c_1 + c_2 $ is a real number and $ c_1 - c_2 $ is an imaginary number.\n",
"\n",
"This happens only when $ c_1 $ and $ c_2 $ are complex conjugates, in which case they can be written in the polar forms\n",
"\n",
"$$\n",
"c_1 = v e^{i \\theta}, \\ \\ c_2 = v e^{- i \\theta}\n",
"$$\n",
"\n",
"So we can write\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" Y_t & = & v e^{i \\theta} r^t e^{i \\omega t} + v e ^{- i \\theta} r^t e^{-i \\omega t} \\\\\n",
" & = & v r^t [ e^{i(\\omega t + \\theta)} + e^{-i (\\omega t +\\theta)}] \\\\\n",
" & = & 2 v r^t \\cos (\\omega t + \\theta)\n",
" \\end{aligned}\n",
"$$\n",
"\n",
"where $ v $ and $ \\theta $ are constants that must be chosen to satisfy initial conditions for $ Y_{-1}, Y_{-2} $.\n",
"\n",
"This formula shows that when the roots are complex, $ Y_t $ displays\n",
"oscillations with *period* $ \\check p =\n",
"\\frac{2 \\pi}{\\omega} $ and *damping factor* $ r $.\n",
"\n",
"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.\n",
"\n",
"(Draw a cosine function to convince yourself of this please)\n",
"\n",
"*Remark:* Following [[Samuelson, 1939](https://python.quantecon.org/zreferences.html#id180)], we want to choose the parameters\n",
"$ \\alpha, \\beta $ of the model so that the absolute values (of the possibly\n",
"complex) roots $ \\lambda_1, \\lambda_2 $ of the characteristic\n",
"polynomial are both strictly less than one:\n",
"\n",
"$$\n",
"| \\lambda_j | < 1 \\quad \\quad \\text{for } j = 1, 2\n",
"$$\n",
"\n",
"*Remark:* When both roots $ \\lambda_1, \\lambda_2 $ of the characteristic polynomial have\n",
"absolute values strictly less than one, the absolute value of the larger\n",
"one governs the rate of convergence to the steady state of the non\n",
"stochastic version of the model."
]
},
{
"cell_type": "markdown",
"id": "97ed7fd4",
"metadata": {},
"source": [
"### Things this lecture does\n",
"\n",
"We write a function to generate simulations of a $ \\{Y_t\\} $ sequence as a function of time.\n",
"\n",
"The function requires that we put in initial conditions for $ Y_{-1}, Y_{-2} $.\n",
"\n",
"The function checks that $ \\alpha, \\beta $ are set so that $ \\lambda_1, \\lambda_2 $ are less than\n",
"unity in absolute value (also called “modulusâ€).\n",
"\n",
"The function also tells us whether the roots are complex, and, if they are complex, returns both their real and complex parts.\n",
"\n",
"If the roots are both real, the function returns their values.\n",
"\n",
"We use our function written to simulate paths that are stochastic (when $ \\sigma >0 $).\n",
"\n",
"We have written the function in a way that allows us to input $ \\{G_t\\} $ paths of a few simple forms, e.g.,\n",
"\n",
"- one time jumps in $ G $ at some time \n",
"- a permanent jump in $ G $ that occurs at some time \n",
"\n",
"\n",
"We proceed to use the Samuelson multiplier-accelerator model as a laboratory to make a simple OOP example.\n",
"\n",
"The “state†that determines next period’s $ Y_{t+1} $ is now not just the current value $ Y_t $ but also the once lagged value $ Y_{t-1} $.\n",
"\n",
"This involves a little more bookkeeping than is required in the Solow model class definition.\n",
"\n",
"We use the Samuelson multiplier-accelerator model as a vehicle for teaching how we can gradually add more features to the class.\n",
"\n",
"We want to have a method in the class that automatically generates a simulation, either non-stochastic ($ \\sigma=0 $) or stochastic ($ \\sigma > 0 $).\n",
"\n",
"We also show how to map the Samuelson model into a simple instance of the `LinearStateSpace` class described [here](https://python.quantecon.org/linear_models.html).\n",
"\n",
"We can use a `LinearStateSpace` instance to do various things that we did above with our homemade function and class.\n",
"\n",
"Among other things, we show by example that the eigenvalues of the matrix $ A $ that we use to form the instance of the `LinearStateSpace` class for the Samuelson model equal the roots of the characteristic polynomial [(41.7)](#equation-polynomial) for the Samuelson multiplier accelerator model.\n",
"\n",
"Here is the formula for the matrix $ A $ in the linear state space system in the case that government expenditures are a constant $ G $:\n",
"\n",
"$$\n",
"A = \\begin{bmatrix} 1 & 0 & 0 \\cr\n",
" \\gamma + G & \\rho_1 & \\rho_2 \\cr\n",
" 0 & 1 & 0 \\end{bmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "0819ac49",
"metadata": {},
"source": [
"## Implementation\n",
"\n",
"We’ll start by drawing an informative graph from page 189 of [[Sargent, 1987](https://python.quantecon.org/zreferences.html#id308)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "35d68641",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def param_plot():\n",
" \"\"\"This function creates the graph on page 189 of\n",
" Sargent Macroeconomic Theory, second edition, 1987.\n",
" \"\"\"\n",
"\n",
" fig, ax = plt.subplots(figsize=(10, 6))\n",
" ax.set_aspect(\"equal\")\n",
"\n",
" # Set axis\n",
" xmin, ymin = -3, -2\n",
" xmax, ymax = -xmin, -ymin\n",
" plt.axis([xmin, xmax, ymin, ymax])\n",
"\n",
" # Set axis labels\n",
" ax.set(xticks=[], yticks=[])\n",
" ax.set_xlabel(r\"$\\rho_2$\", fontsize=16)\n",
" ax.xaxis.set_label_position(\"top\")\n",
" ax.set_ylabel(r\"$\\rho_1$\", rotation=0, fontsize=16)\n",
" ax.yaxis.set_label_position(\"right\")\n",
"\n",
" # Draw (t1, t2) points\n",
" Ï1 = np.linspace(-2, 2, 100)\n",
" ax.plot(Ï1, -abs(Ï1) + 1, c=\"black\")\n",
" ax.plot(Ï1, np.full_like(Ï1, -1), c=\"black\")\n",
" ax.plot(Ï1, -(Ï1**2 / 4), c=\"black\")\n",
"\n",
" # Turn normal axes off\n",
" for spine in [\"left\", \"bottom\", \"top\", \"right\"]:\n",
" ax.spines[spine].set_visible(False)\n",
"\n",
" # Add arrows to represent axes\n",
" axes_arrows = {\"arrowstyle\": \"<|-|>\", \"lw\": 1.3}\n",
" ax.annotate(\"\", xy=(xmin, 0), xytext=(xmax, 0), arrowprops=axes_arrows)\n",
" ax.annotate(\"\", xy=(0, ymin), xytext=(0, ymax), arrowprops=axes_arrows)\n",
"\n",
" # Annotate the plot with equations\n",
" plot_arrowsl = {\"arrowstyle\": \"-|>\", \"connectionstyle\": \"arc3, rad=-0.2\"}\n",
" plot_arrowsr = {\"arrowstyle\": \"-|>\", \"connectionstyle\": \"arc3, rad=0.2\"}\n",
" ax.annotate(\n",
" r\"$\\rho_1 + \\rho_2 < 1$\",\n",
" xy=(0.5, 0.3),\n",
" xytext=(0.8, 0.6),\n",
" arrowprops=plot_arrowsr,\n",
" fontsize=\"12\",\n",
" )\n",
" ax.annotate(\n",
" r\"$\\rho_1 + \\rho_2 = 1$\",\n",
" xy=(0.38, 0.6),\n",
" xytext=(0.6, 0.8),\n",
" arrowprops=plot_arrowsr,\n",
" fontsize=\"12\",\n",
" )\n",
" ax.annotate(\n",
" r\"$\\rho_2 < 1 + \\rho_1$\",\n",
" xy=(-0.5, 0.3),\n",
" xytext=(-1.3, 0.6),\n",
" arrowprops=plot_arrowsl,\n",
" fontsize=\"12\",\n",
" )\n",
" ax.annotate(\n",
" r\"$\\rho_2 = 1 + \\rho_1$\",\n",
" xy=(-0.38, 0.6),\n",
" xytext=(-1, 0.8),\n",
" arrowprops=plot_arrowsl,\n",
" fontsize=\"12\",\n",
" )\n",
" ax.annotate(\n",
" r\"$\\rho_2 = -1$\",\n",
" xy=(1.5, -1),\n",
" xytext=(1.8, -1.3),\n",
" arrowprops=plot_arrowsl,\n",
" fontsize=\"12\",\n",
" )\n",
" ax.annotate(\n",
" r\"${\\rho_1}^2 + 4\\rho_2 = 0$\",\n",
" xy=(1.15, -0.35),\n",
" xytext=(1.5, -0.3),\n",
" arrowprops=plot_arrowsr,\n",
" fontsize=\"12\",\n",
" )\n",
" ax.annotate(\n",
" r\"${\\rho_1}^2 + 4\\rho_2 < 0$\",\n",
" xy=(1.4, -0.7),\n",
" xytext=(1.8, -0.6),\n",
" arrowprops=plot_arrowsr,\n",
" fontsize=\"12\",\n",
" )\n",
"\n",
" # Label categories of solutions\n",
" ax.text(1.5, 1, \"Explosive\\n growth\", ha=\"center\", fontsize=16)\n",
" ax.text(-1.5, 1, \"Explosive\\n oscillations\", ha=\"center\", fontsize=16)\n",
" ax.text(0.05, -1.5, \"Explosive oscillations\", ha=\"center\", fontsize=16)\n",
" ax.text(0.09, -0.5, \"Damped oscillations\", ha=\"center\", fontsize=16)\n",
"\n",
" # Add small marker to y-axis\n",
" ax.axhline(y=1.005, xmin=0.495, xmax=0.505, c=\"black\")\n",
" ax.text(-0.12, -1.12, \"-1\", fontsize=10)\n",
" ax.text(-0.12, 0.98, \"1\", fontsize=10)\n",
"\n",
" return fig\n",
"\n",
"\n",
"param_plot()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "bb8dd0f2",
"metadata": {},
"source": [
"The graph portrays regions in which the $ (\\lambda_1, \\lambda_2) $\n",
"root pairs implied by the $ (\\rho_1 = (\\alpha+\\beta), \\rho_2 = - \\beta) $\n",
"difference equation parameter pairs in the Samuelson model are such that:\n",
"\n",
"- $ (\\lambda_1, \\lambda_2) $ are complex with modulus less than\n",
" $ 1 $ - in this case, the $ \\{Y_t\\} $ sequence displays damped\n",
" oscillations. \n",
"- $ (\\lambda_1, \\lambda_2) $ are both real, but one is strictly\n",
" greater than $ 1 $ - this leads to explosive growth. \n",
"- $ (\\lambda_1, \\lambda_2) $ are both real, but one is strictly\n",
" less than $ -1 $ - this leads to explosive oscillations. \n",
"- $ (\\lambda_1, \\lambda_2) $ are both real and both are less than\n",
" $ 1 $ in absolute value - in this case, there is smooth\n",
" convergence to the steady state without damped cycles. \n",
"\n",
"\n",
"Later we’ll present the graph with a red mark showing the particular\n",
"point implied by the setting of $ (\\alpha,\\beta) $."
]
},
{
"cell_type": "markdown",
"id": "2ea8eeb2",
"metadata": {},
"source": [
"### Function to describe implications of characteristic polynomial"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e1085da2",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def categorize_solution(Ï1, Ï2):\n",
" \"\"\"\n",
" This function takes values of Ï1 and Ï2 and uses them\n",
" to classify the type of solution.\n",
" \"\"\"\n",
" discriminant = Ï1**2 + 4 * Ï2\n",
" if Ï2 > 1 + Ï1 or Ï2 < -1:\n",
" return \"Explosive oscillations\"\n",
" elif Ï1 + Ï2 > 1:\n",
" return \"Explosive growth\"\n",
" elif discriminant < 0:\n",
" return \"Damped oscillations\"\n",
" else:\n",
" return \"Steady state convergence\"\n",
"\n",
"def analyze_roots(α, β, verbose=True):\n",
" \"\"\"\n",
" Unified function to calculate roots and analyze their properties.\n",
" \"\"\"\n",
" Ï1 = α + β\n",
" Ï2 = -β\n",
" \n",
" # Compute characteristic polynomial roots\n",
" roots = np.roots([1, -Ï1, -Ï2])\n",
" \n",
" # Classify solution type\n",
" solution_type = categorize_solution(Ï1, Ï2)\n",
" \n",
" # Determine root properties\n",
" is_complex = all(isinstance(root, complex) for root in roots)\n",
" is_stable = all(abs(root) < 1 for root in roots)\n",
" \n",
" if verbose:\n",
" print(f\"Ï1 = {Ï1:.2f}, Ï2 = {Ï2:.2f}\")\n",
" print(f\"Roots: {[f'{root:.2f}' for root in roots]}\")\n",
" print(f\"Root type: {'Complex' if is_complex else 'Real'}\")\n",
" print(f\"Stability: {'Stable' if is_stable else 'Unstable'}\")\n",
" print(f\"Solution type: {solution_type}\")\n",
" \n",
" return {\n",
" 'roots': roots,\n",
" 'rho1': Ï1,\n",
" 'rho2': Ï2,\n",
" 'is_complex': is_complex,\n",
" 'is_stable': is_stable,\n",
" 'solution_type': solution_type\n",
" }"
]
},
{
"cell_type": "markdown",
"id": "cd7b80e5",
"metadata": {},
"source": [
"We also write a unified simulation function that can handle both\n",
"non-stochastic and stochastic versions of the model.\n",
"\n",
"It allows for government spending paths of a few simple forms which\n",
"we specify via a dictionary `g_params`"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "88bed20b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def simulate_samuelson(\n",
" y_0, y_1, α, β, γ=10, n=100, σ=0, g_params=None, seed=0\n",
"):\n",
" \"\"\"\n",
" Unified simulation function for Samuelson model.\n",
" \n",
" Parameters:\n",
" g_params: dict with keys 'g', 'g_t', 'duration' for government spending\n",
" seed: random seed for reproducible results\n",
" \"\"\"\n",
" analysis = analyze_roots(α, β, verbose=False)\n",
" Ï1, Ï2 = analysis['rho1'], analysis['rho2']\n",
" \n",
" # Initialize time series\n",
" y_t = [y_0, y_1]\n",
" \n",
" # Generate shocks if stochastic\n",
" if σ > 0:\n",
" np.random.seed(seed)\n",
" ϵ = np.random.normal(0, 1, n)\n",
" \n",
" # Simulate forward\n",
" for t in range(2, n):\n",
"\n",
" # Determine government spending\n",
" g = 0\n",
" if g_params:\n",
" g_val, g_t_val = g_params.get('g', 0), g_params.get('g_t', 0)\n",
" duration = g_params.get('duration', None)\n",
" if duration == 'permanent' and t >= g_t_val:\n",
" g = g_val\n",
" elif duration == 'one-off' and t == g_t_val:\n",
" g = g_val\n",
" elif duration is None:\n",
" g = g_val\n",
" \n",
" # Calculate next value\n",
" y_next = Ï1 * y_t[t-1] + Ï2 * y_t[t-2] + γ + g\n",
" if σ > 0:\n",
" y_next += σ * ϵ[t]\n",
" \n",
" y_t.append(y_next)\n",
" \n",
" return y_t, analysis"
]
},
{
"cell_type": "markdown",
"id": "71c1410c",
"metadata": {},
"source": [
"We will use this function to run simulations of the model.\n",
"\n",
"But before doing that, let’s test the analysis function"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2fb3e036",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"analysis = analyze_roots(α=1.3, β=0.4)"
]
},
{
"cell_type": "markdown",
"id": "1606b5d6",
"metadata": {},
"source": [
"### Function for plotting paths\n",
"\n",
"A useful function for our work below is"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c0e61f5d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def plot_y(function=None):\n",
" \"\"\"Function plots path of Y_t\"\"\"\n",
"\n",
" plt.subplots(figsize=(10, 6))\n",
" plt.plot(function)\n",
" plt.xlabel(\"$t$\")\n",
" plt.ylabel(\"$Y_t$\")\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"id": "343dfb03",
"metadata": {},
"source": [
"### Manual or “by hand†root calculations\n",
"\n",
"The following function calculates roots of the characteristic polynomial\n",
"using high school algebra.\n",
"\n",
"(We’ll calculate the roots in other ways later using `analyze_roots`.)\n",
"\n",
"The function also plots a $ Y_t $ starting from initial conditions\n",
"that we set"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "eac1edc8",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def y_nonstochastic(y_0=100, y_1=80, α=0.92, β=0.5, γ=10, n=80):\n",
" \"\"\"\n",
" This function calculates the roots of the characteristic polynomial\n",
" by hand and returns a path of y_t starting from initial conditions\n",
" \"\"\"\n",
" roots = []\n",
"\n",
" Ï1 = α + β\n",
" Ï2 = -β\n",
"\n",
" print(f\"Ï_1 is {Ï1:.2f}\")\n",
" print(f\"Ï_2 is {Ï2:.2f}\")\n",
"\n",
" discriminant = Ï1**2 + 4 * Ï2\n",
"\n",
" if discriminant == 0:\n",
" roots.append(-Ï1 / 2)\n",
" print(\"Single real root: \")\n",
" print(\"\".join(f\"{r:.2f}\" for r in roots))\n",
" elif discriminant > 0:\n",
" roots.append((-Ï1 + sqrt(discriminant).real) / 2)\n",
" roots.append((-Ï1 - sqrt(discriminant).real) / 2)\n",
" print(\"Two real roots: \")\n",
" print(\" \".join(f\"{r:.2f}\" for r in roots))\n",
" else:\n",
" roots.append((-Ï1 + sqrt(discriminant)) / 2)\n",
" roots.append((-Ï1 - sqrt(discriminant)) / 2)\n",
" print(\"Two complex roots: \")\n",
" print(\" \".join(f\"{r.real:.2f}{r.imag:+.2f}j\" for r in roots))\n",
"\n",
" if all(abs(root) < 1 for root in roots):\n",
" print(\"Absolute values of roots are less than one\")\n",
" else:\n",
" print(\"Absolute values of roots are not less than one\")\n",
"\n",
" def transition(x, t):\n",
" return Ï1 * x[t - 1] + Ï2 * x[t - 2] + γ\n",
"\n",
" y_t = [y_0, y_1]\n",
"\n",
" for t in range(2, n):\n",
" y_t.append(transition(y_t, t))\n",
"\n",
" return y_t\n",
"\n",
"\n",
"plot_y(y_nonstochastic())"
]
},
{
"cell_type": "markdown",
"id": "67482afc",
"metadata": {},
"source": [
"### Reverse-engineering parameters to generate damped cycles\n",
"\n",
"The next cell writes code that takes as inputs the modulus $ r $ and\n",
"phase $ \\phi $ of a conjugate pair of complex numbers in polar form\n",
"\n",
"$$\n",
"\\lambda_1 = r \\exp(i \\phi), \\quad \\lambda_2 = r \\exp(- i \\phi)\n",
"$$\n",
"\n",
"- The code assumes that these two complex numbers are the roots of the\n",
" characteristic polynomial \n",
"- It then reverse-engineers $ (\\alpha,\\beta) $ and $ (\\rho_1, \\rho_2) $,\n",
" pairs that would generate those roots "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "371f0807",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def f(r, Ï•):\n",
" \"\"\"\n",
" Takes modulus r and angle Ï• of complex number r exp(j Ï•)\n",
" and creates Ï1 and Ï2 of characteristic polynomial for which\n",
" r exp(j Ï•) and r exp(- j Ï•) are complex roots.\n",
"\n",
" Returns the multiplier coefficient $\\alpha$\n",
" and the accelerator coefficient $\\beta$\n",
" that verifies those roots.\n",
" \"\"\"\n",
" # Create complex conjugate pair from polar coordinates\n",
" g1 = cmath.rect(r, Ï•)\n",
" g2 = cmath.rect(r, -Ï•)\n",
"\n",
" # Calculate corresponding Ï1, Ï2 parameters\n",
" Ï1 = g1 + g2\n",
" Ï2 = -g1 * g2\n",
"\n",
" # Derive α and β coefficients from Ï parameters\n",
" β = -Ï2\n",
" α = Ï1 - β\n",
" return Ï1, Ï2, α, β"
]
},
{
"cell_type": "markdown",
"id": "96f39fab",
"metadata": {},
"source": [
"Now let’s use the function in an example.\n",
"\n",
"Here are the example parameters"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a52f62fe",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"r = 0.95\n",
"\n",
"# Cycle period in time units\n",
"period = 10\n",
"Ï• = 2 * math.pi / period\n",
"\n",
"# Apply the reverse-engineering function\n",
"Ï1, Ï2, α, β = f(r, Ï•)\n",
"\n",
"print(f\"α, β = {α:.2f}, {β:.2f}\")\n",
"print(f\"Ï1, Ï2 = {Ï1:.2f}, {Ï2:.2f}\")"
]
},
{
"cell_type": "markdown",
"id": "5f216eb9",
"metadata": {},
"source": [
"The real parts of the roots are"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "730a9975",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"print(f\"Ï1 = {Ï1.real:.2f}, Ï2 = {Ï2.real:.2f}\")"
]
},
{
"cell_type": "markdown",
"id": "f48fbf6d",
"metadata": {},
"source": [
"### Root finding using numpy\n",
"\n",
"Here we’ll use numpy to compute the roots of the characteristic\n",
"polynomial"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d07e9489",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"r1, r2 = np.roots([1, -Ï1, -Ï2])\n",
"\n",
"p1 = cmath.polar(r1)\n",
"p2 = cmath.polar(r2)\n",
"\n",
"print(f\"r, Ï• = {r:.2f}, {Ï•:.2f}\")\n",
"print(f\"p1, p2 = ({p1[0]:.2f}, {p1[1]:.2f}), ({p2[0]:.2f}, {p2[1]:.2f})\")\n",
"\n",
"print(f\"α, β = {α:.2f}, {β:.2f}\")\n",
"print(f\"Ï1, Ï2 = {Ï1:.2f}, {Ï2:.2f}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3d2b2d4d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def y_nonstochastic(y_0=100, y_1=80, α=0.9, β=0.8, γ=10, n=80):\n",
" \"\"\"\n",
" This function enlists numpy to calculate the roots of the characteristic\n",
" polynomial.\n",
" \"\"\"\n",
"\n",
" y_series, analysis = simulate_samuelson(y_0, y_1, α, β, γ, n, 0, None, 42)\n",
" \n",
" print(f\"Solution type: {analysis['solution_type']}\")\n",
" print(f\"Roots are {analysis['roots']}\")\n",
" print(f\"Root type: {'Complex' if analysis['is_complex'] else 'Real'}\")\n",
" print(f\"Stability: {'Stable' if analysis['is_stable'] else 'Unstable'}\")\n",
" \n",
" return y_series\n",
"\n",
"\n",
"plot_y(y_nonstochastic())"
]
},
{
"cell_type": "markdown",
"id": "c7cf748f",
"metadata": {},
"source": [
"### Reverse-engineered complex roots: Example\n",
"\n",
"The next cell studies the implications of reverse-engineered complex\n",
"roots.\n",
"\n",
"We’ll generate an *undamped* cycle of period 10"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f9f6c60",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"r = 1 # Generates undamped, nonexplosive cycles\n",
"\n",
"period = 10 # Length of cycle in units of time\n",
"Ï• = 2 * math.pi / period\n",
"\n",
"# Apply the reverse-engineering function f\n",
"Ï1, Ï2, α, β = f(r, Ï•)\n",
"\n",
"# Extract real parts for numerical computation\n",
"α = α.real\n",
"β = β.real\n",
"\n",
"print(f\"α, β = {α:.2f}, {β:.2f}\")\n",
"\n",
"ytemp = y_nonstochastic(α=α, β=β, y_0=20, y_1=30)\n",
"plot_y(ytemp)"
]
},
{
"cell_type": "markdown",
"id": "5d003fec",
"metadata": {},
"source": [
"### Digression: Using Sympy to find roots\n",
"\n",
"We can also use sympy to compute analytic formulas for the roots"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8876752d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"init_printing()\n",
"\n",
"r1 = Symbol(\"Ï_1\")\n",
"r2 = Symbol(\"Ï_2\")\n",
"z = Symbol(\"z\")\n",
"\n",
"sympy.solve(z**2 - r1 * z - r2, z)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9a227722",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"α = Symbol(\"α\")\n",
"β = Symbol(\"β\")\n",
"r1 = α + β\n",
"r2 = -β\n",
"\n",
"sympy.solve(z**2 - r1 * z - r2, z)"
]
},
{
"cell_type": "markdown",
"id": "806b253f",
"metadata": {},
"source": [
"## Stochastic shocks\n",
"\n",
"Now we’ll construct some code to simulate the stochastic version of the\n",
"model that emerges when we add a random shock process to aggregate\n",
"demand"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "910b5c9e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def y_stochastic(y_0=0, y_1=0, α=0.8, β=0.2, γ=10, n=100, σ=5):\n",
" \"\"\"\n",
" This function takes parameters of a stochastic version of\n",
" the model, analyzes the roots of the characteristic\n",
" polynomial and generates a simulation.\n",
" \"\"\"\n",
"\n",
" y_series, analysis = simulate_samuelson(y_0, y_1, α, β, γ, n, σ, None, 42)\n",
" \n",
" print(f\"Solution type: {analysis['solution_type']}\")\n",
" print(f\"Roots are {[f'{root:.2f}' for root in analysis['roots']]}\")\n",
" print(f\"Root type: {'Complex' if analysis['is_complex'] else 'Real'}\")\n",
" print(f\"Stability: {'Stable' if analysis['is_stable'] else 'Unstable'}\")\n",
" \n",
" return y_series\n",
"\n",
"\n",
"plot_y(y_stochastic())"
]
},
{
"cell_type": "markdown",
"id": "01a288c9",
"metadata": {},
"source": [
"Let’s do a simulation in which there are shocks and the characteristic\n",
"polynomial has complex roots"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3eba5df3",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"r = 0.97\n",
"\n",
"period = 10 # Length of cycle in units of time\n",
"Ï• = 2 * math.pi / period\n",
"\n",
"# Apply the reverse-engineering function f\n",
"Ï1, Ï2, α, β = f(r, Ï•)\n",
"\n",
"# Extract real parts for numerical computation\n",
"α = α.real\n",
"β = β.real\n",
"\n",
"print(f\"α, β = {α:.2f}, {β:.2f}\")\n",
"plot_y(y_stochastic(y_0=40, y_1=42, α=α, β=β, σ=2, n=100))"
]
},
{
"cell_type": "markdown",
"id": "6c2dc4cc",
"metadata": {},
"source": [
"## Government spending\n",
"\n",
"This function computes a response to either a permanent or one-off increase\n",
"in government expenditures"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57203d3d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def y_stochastic_g(\n",
" y_0=20, y_1=20, α=0.8, β=0.2, γ=10,\n",
" n=100, σ=2, g=0, g_t=0, duration=\"permanent\"\n",
"):\n",
" \"\"\"\n",
" This program computes a response to a permanent increase\n",
" in government expenditures that occurs at time 20\n",
" \"\"\"\n",
"\n",
" g_params = (\n",
" {'g': g, 'g_t': g_t, 'duration': duration} if g != 0 else None\n",
" )\n",
" y_series, analysis = simulate_samuelson(\n",
" y_0, y_1, α, β, γ, n, σ, g_params, 42\n",
" )\n",
" \n",
" print(f\"Solution type: {analysis['solution_type']}\")\n",
" print(f\"Roots: {analysis['roots']}\")\n",
" print(f\"Root type: {'Complex' if analysis['is_complex'] else 'Real'}\")\n",
" print(f\"Stability: {'Stable' if analysis['is_stable'] else 'Unstable'}\")\n",
" \n",
" return y_series"
]
},
{
"cell_type": "markdown",
"id": "c1abd047",
"metadata": {},
"source": [
"A permanent government spending shock can be simulated as follows"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7284c6e7",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_y(y_stochastic_g(g=10, g_t=20, duration=\"permanent\"))"
]
},
{
"cell_type": "markdown",
"id": "9e344fe0",
"metadata": {},
"source": [
"We can also see the response to a one time jump in government expenditures"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2b39356d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_y(y_stochastic_g(g=500, g_t=50, duration=\"one-off\"))"
]
},
{
"cell_type": "markdown",
"id": "644112e0",
"metadata": {},
"source": [
"## Wrapping everything into a class\n",
"\n",
"Up to now, we have written functions to do the work.\n",
"\n",
"Now we’ll roll up our sleeves and write a Python class called `Samuelson`\n",
"for the Samuelson model"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "39b41ed7",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"class Samuelson:\n",
" \"\"\"\n",
" This class represents the Samuelson model, otherwise known as the\n",
" multiplier-accelerator model. \n",
" The model combines the Keynesian multiplier\n",
" with the accelerator theory of investment.\n",
"\n",
" The path of output is governed by a linear \n",
" second-order difference equation\n",
"\n",
" .. math::\n",
"\n",
" Y_t = \\alpha (1 + \\beta) Y_{t-1} - \\alpha \\beta Y_{t-2}\n",
"\n",
" Parameters\n",
" ----------\n",
" y_0 : scalar\n",
" Initial condition for Y_0\n",
" y_1 : scalar\n",
" Initial condition for Y_1\n",
" α : scalar\n",
" Marginal propensity to consume\n",
" β : scalar\n",
" Accelerator coefficient\n",
" n : int\n",
" Number of iterations\n",
" σ : scalar\n",
" Volatility parameter. It must be greater than or equal to 0. Set\n",
" equal to 0 for a non-stochastic model.\n",
" g : scalar\n",
" Government spending shock\n",
" g_t : int\n",
" Time at which government spending shock occurs. Must be specified\n",
" when duration != None.\n",
" duration : {None, 'permanent', 'one-off'}\n",
" Specifies type of government spending shock. If none, government\n",
" spending equal to g for all t.\n",
"\n",
" \"\"\"\n",
"\n",
" def __init__(\n",
" self, y_0=100, y_1=50,\n",
" α=1.3, β=0.2, γ=10, n=100, σ=0, g=0, g_t=0, duration=None\n",
" ):\n",
"\n",
" self.y_0, self.y_1, self.α, self.β = y_0, y_1, α, β\n",
" self.n, self.g, self.g_t, self.duration = n, g, g_t, duration\n",
" self.γ, self.σ = γ, σ\n",
" \n",
" # Use unified analysis function\n",
" self.analysis = analyze_roots(α, β, verbose=False)\n",
" self.Ï1, self.Ï2 = self.analysis['rho1'], self.analysis['rho2']\n",
" self.roots = self.analysis['roots']\n",
"\n",
" def root_type(self):\n",
" return \"Complex conjugate\" if self.analysis['is_complex'] else \"Real\"\n",
"\n",
" def root_less_than_one(self):\n",
" return self.analysis['is_stable']\n",
"\n",
" def solution_type(self):\n",
" return self.analysis['solution_type']\n",
"\n",
" def generate_series(self, seed=0):\n",
" g_params = (\n",
" {'g': self.g, 'g_t': self.g_t, 'duration': self.duration} \n",
" if self.g != 0 else None\n",
" )\n",
" y_series, _ = simulate_samuelson(\n",
" self.y_0, self.y_1, self.α, self.β, self.γ, \n",
" self.n, self.σ, g_params, seed\n",
" )\n",
" return y_series\n",
"\n",
" def summary(self):\n",
" print(\"Summary\\n\" + \"-\" * 50)\n",
" print(f\"Root type: {self.root_type()}\")\n",
" print(f\"Solution type: {self.solution_type()}\")\n",
" print(f\"Roots: {str(self.roots)}\")\n",
"\n",
" if self.root_less_than_one() == True:\n",
" print(\"Absolute value of roots is less than one\")\n",
" else:\n",
" print(\"Absolute value of roots is not less than one\")\n",
"\n",
" if self.σ > 0:\n",
" print(\"Stochastic series with σ = \" + str(self.σ))\n",
" else:\n",
" print(\"Non-stochastic series\")\n",
"\n",
" if self.g != 0:\n",
" print(\"Government spending equal to \" + str(self.g))\n",
"\n",
" if self.duration != None:\n",
" print(\n",
" self.duration.capitalize()\n",
" + \" government spending shock at t = \"\n",
" + str(self.g_t)\n",
" )\n",
"\n",
" def plot(self, seed=0):\n",
" fig, ax = plt.subplots(figsize=(10, 6))\n",
" ax.plot(self.generate_series(seed))\n",
" ax.set(xlabel=\"iteration\", xlim=(0, self.n))\n",
" ax.set_ylabel(\"$Y_t$\", rotation=0)\n",
"\n",
" # Display model parameters on the plot\n",
" paramstr = (\n",
" f\"$α={self.α:.2f}$ \\n $β={self.β:.2f}$ \\n \"\n",
" f\"$\\\\gamma={self.γ:.2f}$ \\n $\\\\sigma={self.σ:.2f}$ \\n \"\n",
" f\"$\\\\rho_1={self.Ï1:.2f}$ \\n $\\\\rho_2={self.Ï2:.2f}$\"\n",
" )\n",
" props = dict(fc=\"white\", pad=10, alpha=0.5)\n",
" ax.text(\n",
" 0.87,\n",
" 0.05,\n",
" paramstr,\n",
" transform=ax.transAxes,\n",
" fontsize=12,\n",
" bbox=props,\n",
" va=\"bottom\",\n",
" )\n",
"\n",
" return fig\n",
"\n",
" def param_plot(self):\n",
"\n",
"\n",
" fig = param_plot()\n",
" ax = fig.gca()\n",
"\n",
" # Display eigenvalues in the legend\n",
" for i, root in enumerate(self.roots):\n",
" if isinstance(root, complex):\n",
"\n",
" # Handle sign formatting for complex number display\n",
" operator = [\"+\", \"\"]\n",
" root_real = self.roots[i].real\n",
" root_imag = self.roots[i].imag\n",
" label = (\n",
" rf\"$\\lambda_{i+1} = {root_real:.2f}\"\n",
" rf\"{operator[i]} {root_imag:.2f}i$\"\n",
" )\n",
" else:\n",
" label = rf\"$\\lambda_{i+1} = {self.roots[i].real:.2f}$\"\n",
"\n",
" # Add invisible point for legend entry\n",
" ax.scatter(\n",
" 0, 0, s=0, label=label\n",
" )\n",
"\n",
" # Mark current parameter values on the stability diagram\n",
" ax.scatter(\n",
" self.Ï1,\n",
" self.Ï2,\n",
" s=100,\n",
" c=\"red\",\n",
" marker=\"+\",\n",
" label=r\"$(\\rho_1, \\rho_2)$\",\n",
" zorder=5,\n",
" )\n",
"\n",
" plt.legend(fontsize=12, loc=3)\n",
"\n",
" return fig"
]
},
{
"cell_type": "markdown",
"id": "2601b10c",
"metadata": {},
"source": [
"### Illustration of Samuelson class\n",
"\n",
"Now we’ll put our Samuelson class to work on an example"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dcf00216",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"sam = Samuelson(α=0.8, β=0.5, σ=2, g=10, g_t=20, duration=\"permanent\")\n",
"sam.summary()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aaefe5bb",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"sam.plot()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "0c20f2ef",
"metadata": {},
"source": [
"### Using the graph\n",
"\n",
"We’ll use our graph to show where the roots lie and how their location\n",
"is consistent with the behavior of the path just graphed.\n",
"\n",
"The red $ + $ sign shows the location of the roots"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ec733939",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"sam.param_plot()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "752ced05",
"metadata": {},
"source": [
"## Using the LinearStateSpace class\n",
"\n",
"It turns out that we can use the [QuantEcon.py](https://quantecon.org/quantecon-py/)\n",
"[LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to do\n",
"much of the work that we have done from scratch above.\n",
"\n",
"Here is how we map the Samuelson model into an instance of a\n",
"`LinearStateSpace` class"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28cd0500",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"α = 0.8\n",
"β = 0.9\n",
"Ï1 = α + β\n",
"Ï2 = -β\n",
"γ = 10\n",
"σ = 1\n",
"g = 10\n",
"n = 100\n",
"\n",
"A = [[1, 0, 0], [γ + g, Ï1, Ï2], [0, 1, 0]]\n",
"\n",
"G = [\n",
" [γ + g, Ï1, Ï2], # Y_{t+1}\n",
" [γ, α, 0], # C_{t+1}\n",
" [0, β, -β], # I_{t+1}\n",
"]\n",
"\n",
"μ_0 = [1, 100, 50]\n",
"C = np.zeros((3, 1))\n",
"C[1] = σ # Shock variance\n",
"\n",
"sam_t = LinearStateSpace(A, C, G, mu_0=μ_0)\n",
"\n",
"x, y = sam_t.simulate(ts_length=n)\n",
"\n",
"fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8))\n",
"titles = [\"Output ($Y_t$)\", \"Consumption ($C_t$)\", \"Investment ($I_t$)\"]\n",
"colors = [\"darkblue\", \"red\", \"purple\"]\n",
"for ax, series, title, color in zip(axes, y, titles, colors):\n",
" ax.plot(series, color=color)\n",
" ax.set(title=title, xlim=(0, n))\n",
"\n",
"axes[-1].set_xlabel(\"iteration\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "b5ac59b4",
"metadata": {},
"source": [
"### Other methods in the `LinearStateSpace` class\n",
"\n",
"Let’s plot *impulse response functions* for the instance of the\n",
"Samuelson model using a method in the `LinearStateSpace` class"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0e0bae4e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"imres = sam_t.impulse_response()\n",
"imres = np.asarray(imres)\n",
"y1 = imres[:, :, 0]\n",
"y2 = imres[:, :, 1]\n",
"y1.shape"
]
},
{
"cell_type": "markdown",
"id": "981e7bc7",
"metadata": {},
"source": [
"Now let’s compute the zeros of the characteristic polynomial by simply\n",
"calculating the eigenvalues of $ A $"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "85ab24c2",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"A = np.asarray(A)\n",
"w, v = np.linalg.eig(A)\n",
"print(np.round(w, 2))"
]
},
{
"cell_type": "markdown",
"id": "cb7ebbad",
"metadata": {},
"source": [
"### Inheriting methods from `LinearStateSpace`\n",
"\n",
"We could also create a subclass of `LinearStateSpace` (inheriting all its\n",
"methods and attributes) to add more functions to use"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f205887",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"class SamuelsonLSS(LinearStateSpace):\n",
" \"\"\"\n",
" This subclass creates a Samuelson multiplier-accelerator model\n",
" as a linear state space system.\n",
" \"\"\"\n",
"\n",
" def __init__(self, y_0=100, y_1=50, α=0.8, β=0.9, γ=10, σ=1, g=10):\n",
"\n",
" self.α, self.β = α, β\n",
" self.y_0, self.y_1, self.g = y_0, y_1, g\n",
" self.γ, self.σ = γ, σ\n",
"\n",
" # Set initial state vector\n",
" self.initial_μ = [1, y_0, y_1]\n",
"\n",
" self.Ï1 = α + β\n",
" self.Ï2 = -β\n",
"\n",
" # Construct state transition matrix\n",
" self.A = [[1, 0, 0], [γ + g, self.Ï1, self.Ï2], [0, 1, 0]]\n",
"\n",
" # Construct observation matrix\n",
" self.G = [\n",
" [γ + g, self.Ï1, self.Ï2], # Y_{t+1}\n",
" [γ, α, 0], # C_{t+1}\n",
" [0, β, -β], # I_{t+1}\n",
" ] \n",
"\n",
" self.C = np.zeros((3, 1))\n",
" self.C[1] = σ # Shock variance\n",
"\n",
" # Initialize the LinearStateSpace instance\n",
" LinearStateSpace.__init__(\n",
" self, self.A, self.C, self.G, mu_0=self.initial_μ\n",
" )\n",
"\n",
" # Create unicode aliases for mu_0 and Sigma_0 in the parent class\n",
" @property\n",
" def μ_0(self):\n",
" return self.mu_0\n",
" \n",
" @μ_0.setter\n",
" def μ_0(self, value):\n",
" self.mu_0 = value\n",
" \n",
" @property\n",
" def Σ_0(self):\n",
" return self.Sigma_0\n",
" \n",
" @Σ_0.setter\n",
" def Σ_0(self, value):\n",
" self.Sigma_0 = value\n",
"\n",
" def plot_simulation(self, ts_length=100, stationary=True, seed=0):\n",
"\n",
" # Store original distribution parameters\n",
" temp_μ = self.μ_0\n",
" temp_Σ = self.Σ_0\n",
"\n",
" # Use stationary distribution for simulation\n",
" if stationary == True:\n",
" try:\n",
" (self.μ_x, self.μ_y, self.Σ_x, self.Σ_y, self.Σ_yx\n",
" ) = self.stationary_distributions()\n",
" self.μ_0 = self.μ_x\n",
" self.Σ_0 = self.Σ_x\n",
"\n",
" # Handle case where stationary distribution doesn't exist\n",
" except ValueError:\n",
" print(\"Stationary distribution does not exist\")\n",
"\n",
" np.random.seed(seed)\n",
" x, y = self.simulate(ts_length)\n",
"\n",
" fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8))\n",
" titles = [\"Output ($Y_t$)\", \n",
" \"Consumption ($C_t$)\", \n",
" \"Investment ($I_t$)\"]\n",
" colors = [\"darkblue\", \"red\", \"purple\"]\n",
" for ax, series, title, color in zip(axes, y, titles, colors):\n",
" ax.plot(series, color=color)\n",
" ax.set(title=title, xlim=(0, n))\n",
"\n",
" axes[-1].set_xlabel(\"iteration\")\n",
" plt.show()\n",
"\n",
" # Restore original distribution parameters\n",
" self.μ_0 = temp_μ\n",
" self.Σ_0 = temp_Σ\n",
"\n",
" def plot_irf(self, j=5):\n",
"\n",
" x, y = self.impulse_response(j)\n",
"\n",
" # Reshape impulse responses for plotting\n",
" yimf = np.array(y).flatten().reshape(j + 1, 3).T\n",
"\n",
" fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 8))\n",
" labels = [\"$Y_t$\", \"$C_t$\", \"$I_t$\"]\n",
" colors = [\"darkblue\", \"red\", \"purple\"]\n",
" for ax, series, label, color in zip(axes, yimf, labels, colors):\n",
" ax.plot(series, color=color)\n",
" ax.set(xlim=(0, j))\n",
" ax.set_ylabel(label, rotation=0, fontsize=14, labelpad=10)\n",
"\n",
" axes[0].set_title(\"Impulse response functions\")\n",
" axes[-1].set_xlabel(\"iteration\")\n",
" plt.show()\n",
"\n",
" def multipliers(self, j=5):\n",
" x, y = self.impulse_response(j)\n",
" return np.sum(np.array(y).flatten().reshape(j + 1, 3), axis=0)"
]
},
{
"cell_type": "markdown",
"id": "565e3edc",
"metadata": {},
"source": [
"### Illustrations\n",
"\n",
"Let’s show how we can use the `SamuelsonLSS`"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f532ce7f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"samlss = SamuelsonLSS()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ca79a34",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"samlss.plot_simulation(100, stationary=False)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e89eee8e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"samlss.plot_simulation(100, stationary=True)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e9442acd",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"samlss.plot_irf(100)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b1d63ff6",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"samlss.multipliers()"
]
},
{
"cell_type": "markdown",
"id": "80294a2d",
"metadata": {},
"source": [
"## Pure multiplier model\n",
"\n",
"Let’s shut down the accelerator by setting $ b=0 $ to get a pure\n",
"multiplier model\n",
"\n",
"- the absence of cycles gives an idea about why Samuelson included the\n",
" accelerator "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f1e342be",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"pure_multiplier = SamuelsonLSS(α=0.95, β=0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5540004f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"pure_multiplier.plot_simulation()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bb1fbbc2",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"pure_multiplier = SamuelsonLSS(α=0.8, β=0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "04d1615d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"pure_multiplier.plot_simulation()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7c33c505",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"pure_multiplier.plot_irf(100)"
]
},
{
"cell_type": "markdown",
"id": "10d447a4",
"metadata": {},
"source": [
"## Summary\n",
"\n",
"In this lecture, we wrote functions and classes to represent non-stochastic and\n",
"stochastic versions of the Samuelson (1939) multiplier-accelerator model, described\n",
"in [[Samuelson, 1939](https://python.quantecon.org/zreferences.html#id180)].\n",
"\n",
"We saw that different parameter values led to different output paths, which\n",
"could either be stationary, explosive, or oscillating.\n",
"\n",
"We also were able to represent the model using the [QuantEcon.py](https://quantecon.org/quantecon-py/)\n",
"[LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class."
]
}
],
"metadata": {
"date": 1782116504.004488,
"filename": "samuelson.md",
"kernelspec": {
"display_name": "Python",
"language": "python3",
"name": "python3"
},
"title": "Samuelson Multiplier-Accelerator"
},
"nbformat": 4,
"nbformat_minor": 5
}\n",
" \n",
"